C This program builds up a solid object using a cubic array of N entities to form a
C truncated cube to be employed in DDSCAT. 

        PROGRAM PENTAGONAL_PRISM


C TERMS AND CONDITIONS
C This program has been devised for its use in research within the academic commnunity. 
C The authors are: Ana Lilia Gonzalez and Cecilia Noguez. There is no warranty for any 
C part of the program, which is provided by the authors without warranty of any 
C kind, either express or implied.
C The use of this program in any form, subject to the restriction that in any paper or 
C other academic publication containing results wholly or partially derived from the 
C results of use of this program, the following papers must be cited in the normal manner:
C  1. "Influence of morphology on the optical properties of metal nanoparticles" 
C      A. L. Gonzalez and C. Noguez, J. Comput. Theor. Nanosci. 4, 231 - 238 (2007).
C  2. "Optical properties of elongated noble metal nanoparticles",
C	   A. L. Gonzalez, J. A. Reyes-Esqueda, and C. Noguez, 
C          Journal of Physical Chemistry C 112 (19), 7356 -736 (2008). 
C Other papers relevant to this program and method should be cited as best 
C scientific practice dictates.
C The authors do not supply support services of any kind, including, but not 
C limited to, installation, compilation and use of the software.
C Please report errors and modifications to cecilia@fisica.unam.mx


      !TARGET ORIENTATION: The prism is along  the Z axis
      !One edge of the pentagonal cross section is parallel to the x axis.
      !The opposite vertex to that edge is parallel the y axis


      parameter(mxnat=950000)
      PARAMETER(Nver=12)  
      parameter(nnor=17)     
      parameter(m=2)
      IMPLICIT REAL*8 (A-H,O-Q,T-Z) 
      REAL*8 AL,A,B,PI,X,Y,Z,alphamed,betha,modu,apot
      REAL*8 D,H,hprin,GAMMA,vdirx,vdiry,vdirz,vx,vy,vz
      real*8 dnx,dny,dnz,vv,norma,vvdir,dnn,razon
      INTEGER NH,NA,NB,IY,IX,IZ,IXYZ,NAT,mxnat,NZZ,nver,nnor,m,i2
      DIMENSION IXYZ(mxNAT,3)  
      DIMENSION vdirx(m),vdiry(m),vdirz(m),dnn(Nnor)
      dimension vv(nver),norma(nnor),vvdir(m),modu(nnor)       
      DIMENSION VX(nver),VY(Nver),VZ(Nver),DNX(nnor),DNY(nnor),DNZ(Nnor) 
   

      WRITE(*,*)'Give the aspect ratio of the prism'
      READ(*,*)RAZON

      If(razon.le.1)then
         write(*,*)'error, value not permited.'
         stop'end of program'
      else if (razon.gt.1) then 
         write(*,*)'aspect ratio of the prism =',razon
      end if

      write(*,*)'Give a length to generate N dipole entities.'
      write(*,*)'For example, if aspect ratio=2'
      write(*,*)'length = 10 then N = '
      write(*,*)'length = 15 then N = '
      write(*,*)'length = 20 then N = ' 
      read(*,*)al

      If(al.le.0)then
         stop'error, value not permited'
      end if

      PI = 4.d0*datan(1.d0)
      alphamed=  3.d0*PI/10.d0 
      betha=4.d0*PI/10.d0  
      GAMMA=PI/5.d0      
      b=Al*Dcos(alphamed)  
      d=Al*Dcos(betha)
      h=Al*Dsin(betha)
      apot=Al/(2.d0*Dtan(gamma)) 
      hprin=Dsqrt(3.d0*(al/2.d0)**2-apot**2)

      OPEN(UNIT=1,FILE='coordinates.dat') 
      OPEN(UNIT=2,FILE='pentaprism_normales.dat')


      !ojo la razon que doy aqui no incluye la altura de las puntas,
      !es decir, estoy midiendo el ancho del pentagono entre la longitud
      !de las caras rectangulares del prisma.

c****************************************************************************
          !Vertex coordinates of the pentagonal prism

          
         vx(1)=-al/2.d0-d  
         vy(1)=dsqrt(apot**2-d**2-d*al)   
         vz(1)=(al+2.d0*d)*razon/2.d0 -hprin 
         vv(1)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2)

         vx(2)=0.d0  
         vy(2)=b+h-apot 
         vz(2)=(al+2.d0*d)*razon/2.d0-hprin 
         vv(2)=dsqrt(vx(2)**2+vy(2)**2+vz(2)**2)

         vx(3)=al/2.d0+d  
         vy(3)=dsqrt(apot**2-d**2-d*al)  
         vz(3)=(al+2.d0*d)*razon/2.d0-hprin  
         vv(3)=dsqrt(vx(3)**2+vy(3)**2+vz(3)**2)

         vx(4)=al/2.d0  
         vy(4)=-apot
         vz(4)=(al+2.d0*d)*razon/2.d0-hprin   
         vv(4)=dsqrt(vx(4)**2+vy(4)**2+vz(4)**2)

         vx(5)=-al/2.d0  
         vy(5)=-apot
         vz(5)=(al+2.d0*d)*razon/2.d0-hprin 
         vv(5)=dsqrt(vx(5)**2+vy(5)**2+vz(5)**2)

         vx(6)=0.d0  
         vy(6)=0.d0
         vz(6)=(al+2.d0*d)*razon/2.d0   
         vv(6)=dsqrt(vx(6)**2+vy(6)**2+vz(6)**2)

         vx(7)=0.d0  
         vy(7)=0.d0 
         vz(7)=-(al+2.d0*d)*razon/2.d0
         vv(7)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2 )

         vx(8)=-al/2.d0-d    
         vy(8)=dsqrt(apot**2-d**2-d*al)   
         vz(8)=-(al+2.d0*d)*razon/2.d0 +hprin 
         vv(8)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2) 

         vx(9)=0.d0  
         vy(9)=b+h-apot 
         vz(9)=-(al+2.d0*d)*razon/2.d0+hprin 
         vv(9)=dsqrt(vx(2)**2+vy(2)**2+vz(2)**2)

         vx(10)=al/2.d0+d  
         vy(10)=dsqrt(apot**2-d**2-d*al)  
         vz(10)=-(al+2.d0*d)*razon/2.d0+hprin 
         vv(10)=dsqrt(vx(3)**2+vy(3)**2+vz(3)**2)

         vx(11)=al/2.d0  
         vy(11)=-apot
         vz(11)=-(al+2.d0*d)*razon/2.d0 +hprin 
         vv(11)=dsqrt(vx(4)**2+vy(4)**2+vz(4)**2)

         vx(12)=-al/2.d0  
         vy(12)=-apot
         vz(12)=-(al+2.d0*d)*razon/2.d0+hprin 
         vv(12)=dsqrt(vx(5)**2+vy(5)**2+vz(5)**2)


C***************************************************************************
          !Normal vectors to each face of the prism

          Do 45 i=1,4 

         Vdirx(1)=(vx(i+1)-vx(i))/dsqrt((vx(i+1)-vx(i))**2+
     *    (vy(i+1)-vy(i))**2+(vz(i+1)-vz(i))**2)
         Vdiry(1)=(vy(i+1)-vy(i))/dsqrt((vx(i+1)-vx(i))**2+
     *    (vy(i+1)-vy(i))**2+(vz(i+1)-vz(i))**2)
         vdirz(1)=(vz(i+1)-vz(i))/dsqrt((vx(i+1)-vx(i))**2+
     *    (vy(i+1)-vy(i))**2+(vz(i+1)-vz(i))**2)


         Vdirx(2)=(vx(6)-vx(i))/dsqrt((vx(6)-vx(i))**2+
     *    (vy(6)-vy(i))**2+(vz(6)-vz(i))**2)
         Vdiry(2)=(vy(6)-vy(i))/dsqrt((vx(6)-vx(i))**2+
     *    (vy(6)-vy(i))**2+(vz(6)-vz(i))**2)
         vdirz(2)=(vz(6)-vz(i))/dsqrt((vx(6)-vx(i))**2+
     *    (vy(6)-vy(i))**2+(vz(6)-vz(i))**2)


         modu(i)=dsqrt((vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))**2+
     * (vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))**2+
     * (vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))**2)
      
         DNx(i)=(vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))/modu(i)
         DNy(i)=(vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))/modu(i)
         Dnz(i)=(vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))/modu(i)

         Dnn(i)=Dsqrt(dnx(i)**2+dny(i)**2+dnz(i)**2)

 3       FORMAT('Nx=',F13.9,4x,'Ny=',F13.9,4x,'Nz=',F13.9)

         write(2,3)dnx(i),dny(i),dnz(i)
         
         
 45      continue

            
         Vdirx(1)=(vx(1)-vx(5))/dsqrt((vx(1)-vx(5))**2+
     *    (vy(1)-vy(5))**2+(vz(1)-vz(5))**2)
         Vdiry(1)=(vy(1)-vy(5))/dsqrt((vx(1)-vx(5))**2+
     *    (vy(1)-vy(5))**2+(vz(1)-vz(5))**2)
         vdirz(1)=(vz(1)-vz(5))/dsqrt((vx(1)-vx(5))**2+
     *    (vy(1)-vy(5))**2+(vz(1)-vz(5))**2)


         Vdirx(2)=(vx(6)-vx(5))/dsqrt((vx(6)-vx(5))**2+
     *    (vy(6)-vy(5))**2+(vz(6)-vz(5))**2)
         Vdiry(2)=(vy(6)-vy(5))/dsqrt((vx(6)-vx(5))**2+
     *    (vy(6)-vy(5))**2+(vz(6)-vz(5))**2)
         vdirz(2)=(vz(6)-vz(5))/dsqrt((vx(6)-vx(5))**2+
     *    (vy(6)-vy(5))**2+(vz(6)-vz(5))**2)


         modu(5)=dsqrt((vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))**2+
     * (vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))**2+
     * (vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))**2)
      
         DNx(5)=(vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))/modu(5)
         DNy(5)=(vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))/modu(5)
         Dnz(5)=(vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))/modu(5)

         Dnn(5)=Dsqrt(dnx(5)**2+dny(5)**2+dnz(5)**2)
         
         write(2,3)dnx(5),dny(5),dnz(5)
         
         DNx(6)=dnx(5)
         dny(6)=dny(5)
         dnz(6)=-1.*dnz(5)
         dnn(6)=dnn(5)
         write(2,3)dnx(6),dny(6),dnz(6)
 
           Do 55 i=1,4

         Vdirx(1)=(vx(i+8)-vx(i+7))/dsqrt((vx(i+8)-vx(i+7))**2+
     *    (vy(i+8)-vy(i+7))**2+(vz(i+8)-vz(i+7))**2)
         Vdiry(1)=(vy(i+8)-vy(i+7))/dsqrt((vx(i+8)-vx(i+7))**2+
     *    (vy(i+8)-vy(i+7))**2+(vz(i+8)-vz(i+7))**2)
         vdirz(1)=(vz(i+8)-vz(i+7))/dsqrt((vx(i+8)-vx(i+7))**2+
     *    (vy(i+8)-vy(i+7))**2+(vz(i+8)-vz(i+7))**2)


         Vdirx(2)=(vx(7)-vx(i+7))/dsqrt((vx(7)-vx(i+7))**2+
     *    (vy(7)-vy(i+7))**2+(vz(7)-vz(i+7))**2)
         Vdiry(2)=(vy(7)-vy(i+7))/dsqrt((vx(7)-vx(i+7))**2+
     *    (vy(7)-vy(i+7))**2+(vz(7)-vz(i+7))**2)
         vdirz(2)=(vz(7)-vz(i+7))/dsqrt((vx(7)-vx(i+7))**2+
     *    (vy(7)-vy(i+7))**2+(vz(7)-vz(i+7))**2)


         modu(i)=dsqrt((vdiry(1)*vdirz(2)-vdirz(1)*vdiry(2))**2+
     * (vdirz(1)*vdirx(2)-vdirx(1)*vdirz(2))**2+
     * (vdirx(1)*vdiry(2)-vdiry(1)*vdirx(2))**2)
      
         DNx(i+6)=(vdiry(1)*vdirz(2)-vdirz(1)*vdiry(2))/modu(i)
         DNy(i+6)=(vdirz(1)*vdirx(2)-vdirx(1)*vdirz(2))/modu(i)
         Dnz(i+6)=(vdirx(1)*vdiry(2)-vdiry(1)*vdirx(2))/modu(i)

         Dnn(i+6)=Dsqrt(dnx(i+6)**2+dny(i+6)**2+dnz(i+6)**2)

         write(2,3)dnx(i+6),dny(i+6),dnz(i+6)     

 55      continue


          Do 66 i=1,4

         Vdirx(1)=(vx(i+7)-vx(i))/dsqrt((vx(i+7)-vx(i))**2+
     *    (vy(i+7)-vy(i))**2+(vz(i+7)-vz(i))**2)
         Vdiry(1)=(vy(i+7)-vy(i))/dsqrt((vx(i+7)-vx(i))**2+
     *    (vy(i+7)-vy(i))**2+(vz(i+7)-vz(i))**2)
         vdirz(1)=(vz(i+7)-vz(i))/dsqrt((vx(i+7)-vx(i))**2+
     *    (vy(i+7)-vy(i))**2+(vz(i+7)-vz(i))**2)

         Vdirx(2)=(vx(i+1)-vx(i))/dsqrt((vx(i+1)-vx(i))**2+
     *    (vy(i+1)-vy(i))**2+(vz(i+1)-vz(i))**2)
         Vdiry(2)=(vy(i+1)-vy(i))/dsqrt((vx(i+1)-vx(i))**2+
     *    (vy(i+1)-vy(i))**2+(vz(i+1)-vz(i))**2)
         vdirz(2)=(vz(i+1)-vz(i))/dsqrt((vx(i+1)-vx(i))**2+
     *    (vy(i+1)-vy(i))**2+(vz(i+1)-vz(i))**2)


         modu(i)=dsqrt((vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))**2+
     * (vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))**2+
     * (vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))**2)
      
         DNx(i+10)=(vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))/modu(i)
         DNy(i+10)=(vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))/modu(i)
         Dnz(i+10)=(vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))/modu(i)

         Dnn(i+10)=Dsqrt(dnx(i+10)**2+dny(i+10)**2+dnz(i+10)**2)

         write(2,3)dnx(i+10),dny(i+10),dnz(i+10)
         
         
 66   continue


         Vdirx(1)=(vx(12)-vx(5))/dsqrt((vx(12)-vx(5))**2+
     *    (vy(12)-vy(5))**2+(vz(12)-vz(5))**2)
         Vdiry(1)=(vy(12)-vy(5))/dsqrt((vx(12)-vx(5))**2+
     *    (vy(12)-vy(5))**2+(vz(12)-vz(5))**2)
         vdirz(1)=(vz(12)-vz(5))/dsqrt((vx(12)-vx(5))**2+
     *    (vy(12)-vy(5))**2+(vz(12)-vz(5))**2)


         Vdirx(2)=(vx(1)-vx(5))/dsqrt((vx(1)-vx(5))**2+
     *    (vy(1)-vy(5))**2+(vz(1)-vz(5))**2)
         Vdiry(2)=(vy(1)-vy(5))/dsqrt((vx(1)-vx(5))**2+
     *    (vy(1)-vy(5))**2+(vz(1)-vz(5))**2)
         vdirz(2)=(vz(1)-vz(5))/dsqrt((vx(1)-vx(5))**2+
     *    (vy(1)-vy(5))**2+(vz(1)-vz(5))**2)


         modu(i)=dsqrt((vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))**2+
     * (vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))**2+
     * (vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))**2)
      
         DNx(15)=(vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))/modu(i)
         DNy(15)=(vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))/modu(i)
         Dnz(15)=(vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))/modu(i)

         Dnn(15)=Dsqrt(dnx(15)**2+dny(15)**2+dnz(15)**2)

         write(2,3)dnx(15),dny(15),dnz(15)
         
C**********************************************************************
         !Building up the pentagonal prism       
     
         NB=INT(apot)  
         NA=INT(vx(3))
         NZZ=INT((al+2.d0*d)*razon/2.d0+hprin) 
         Nh=INT(vy(2))
          NAT=0
      
      DO 80 IZ=-NZZ,NZZ

         z=real(iz)
      DO 10 IY=-NB,Nh
         Y=REAL(IY)

         DO  IX=-NA,NA
            
        NAT=NAT+1
               X=REAL(IX)

      If (((x-vx(1))*dnx(1)+(y-vy(1))*DNy(1)+
     * (z-vz(1))*dnz(1).le.0).and.
     * ((x-vx(2))*dnx(2)+(y-vy(2))*DNy(2)+(z-vz(2))*dnz(2).le.0).and.
     * ((x-vx(3))*dnx(3)+(y-vy(3))*DNy(3)+(z-vz(3))*dnz(3).le.0).and.
     * ((x-vx(4))*dnx(4)+(y-vy(4))*DNy(4)+(z-vz(4))*dnz(4).le.0).and.
     * ((x-vx(5))*dnx(5)+(y-vy(5))*DNy(5)+(z-vz(5))*dnz(5).le.0).and.
     * ((x-vx(12))*dnx(6)+(y-vy(12))*DNy(6)+(z-vz(12))*dnz(6).le.0).and.
     * ((x-vx(8))*dnx(7)+(y-vy(8))*DNy(7)+(z-vz(8))*dnz(7).le.0).and.
     * ((x-vx(9))*dnx(8)+(y-vy(9))*DNy(8)+(z-vz(9))*dnz(8).le.0).and.
     * ((x-vx(10))*dnx(9)+(y-vy(10))*DNy(9)+(z-vz(10))*dnz(9).le.0).and.
     * ((x-vx(11))*dnx(10)+(y-vy(11))*DNy(10)+
     *   (z-vz(11))*dnz(10).le.0).and.
     * ((x-vx(1))*dnx(11)+(y-vy(1))*DNy(11)+(z-vz(1))*dnz(11).le.0).and.
     * ((x-vx(2))*dnx(12)+(y-vy(2))*DNy(12)+(z-vz(2))*dnz(12).le.0).and.
     * ((x-vx(3))*dnx(13)+(y-vy(3))*DNy(13)+(z-vz(3))*dnz(13).le.0).and.
     * ((x-vx(4))*dnx(14)+(y-vy(4))*DNy(14)+(z-vz(4))*dnz(14).le.0).and.
     * ((x-vx(5))*dnx(15)+(y-vy(5))*DNy(15)+(z-vz(5))*dnz(15).le.0))
     *   IXYZ(NAT,1)=INT(X) 
        
         If (((x-vx(1))*dnx(1)+(y-vy(1))*DNy(1)
     *      +(z-vz(1))*dnz(1).gt.0).or.
     * ((x-vx(2))*dnx(2)+(y-vy(2))*DNy(2)+(z-vz(2))*dnz(2).gt.0).or.
     * ((x-vx(3))*dnx(3)+(y-vy(3))*DNy(3)+(z-vz(3))*dnz(3).gt.0).or.
     * ((x-vx(4))*dnx(4)+(y-vy(4))*DNy(4)+(z-vz(4))*dnz(4).gt.0).or.
     * ((x-vx(5))*dnx(5)+(y-vy(5))*DNy(5)+(z-vz(5))*dnz(5).gt.0).or.
     * ((x-vx(12))*dnx(6)+(y-vy(12))*DNy(6)+(z-vz(12))*dnz(6).gt.0).or.
     * ((x-vx(8))*dnx(7)+(y-vy(8))*DNy(7)+(z-vz(8))*dnz(7).gt.0).or.
     * ((x-vx(9))*dnx(8)+(y-vy(9))*DNy(8)+(z-vz(9))*dnz(8).gt.0).or.
     * ((x-vx(10))*dnx(9)+(y-vy(10))*DNy(9)+(z-vz(10))*dnz(9).gt.0).or.
     * ((x-vx(11))*dnx(10)+(y-vy(11))*DNy(10)+
     * (z-vz(11))*dnz(10).gt.0).or.
     * ((x-vx(1))*dnx(11)+(y-vy(1))*DNy(11)+(z-vz(1))*dnz(11).gt.0).or.
     * ((x-vx(2))*dnx(12)+(y-vy(2))*DNy(12)+(z-vz(2))*dnz(12).gt.0).or.
     * ((x-vx(3))*dnx(13)+(y-vy(3))*DNy(13)+(z-vz(3))*dnz(13).gt.0).or.
     * ((x-vx(4))*dnx(14)+(y-vy(4))*DNy(14)+(z-vz(4))*dnz(14).gt.0).or.
     * ((x-vx(5))*dnx(15)+(y-vy(5))*DNy(15)+(z-vz(5))*dnz(15).gt.0))
     *   goto 20
      
          IXYZ(NAT,2)=INT(y)
          IXYZ(NAT,3)=int(z)
          
          WRITE(1,*) IXYZ(NAT,1),IXYZ(NAT,2),IXYZ(nat,3) !x,y,z coordinates

 20          end do
         
 10       CONTINUE
 80       CONTINUE

 50       STOP
    
        END

     
     
    
         
      

    
