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 HEXAGONAL_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. 
      !Two edges of the hexagonal cross section are parallel to the the x axis.
   

      parameter(mxnat=950000)                             
      PARAMETER(Nver=14)     
      parameter(nnor=18)
      parameter(m=2) 
      REAL*8 AL,AL_2,PI,X,Y,Z,modu,apot,r,ancho_hexa
      REAL*8 hprin,vdirx,vdiry,vdirz,vx,vy,vz
      real*8 dnx,dny,dnz,vv,norma,vvdir,dnn,razon
      INTEGER NX,NY,NZ,IY,IX,IZ,IXYZ,NAT,mxnat,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 = 6,871 '
      write(*,*)'length = 15 then N = 22,961'
      write(*,*)'length = 20 then N = 55,701' 
      read(*,*)r

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

      PI = 4.d0*datan(1.d0)


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

        
C*************************************************************************
         !Vertex of the prism
      
         vx(1)=-r    
         vy(1)=0.D0   
         vz(1)=0.D0  
         vv(1)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2) 


         Do 10 i=1,5 

         vx(i+1)=vx(i)*Dcos(pi/3.d0)+vy(i)*Dsin(pi/3.d0)    
         vy(i+1)=-vx(i)*Dsin(pi/3.d0)+vy(i)*Dcos(pi/3.d0) 
         vz(i+1)=0.d0


  10      continue
 
         al=dsqrt((vx(2)-vx(1))**2+(vy(2)-vy(1))**2)

         
         hprin=r   !height of the ends

        
         vx(7)=0.D0  
         vy(7)=0.D0   
         vz(7)=-hprin  


         ANCHO_HEXA=VX(4)-VX(1)

         AL_2=RAZON*ANCHO_HEXA-2.D0*hprin


         Do 20 i=1,6

         vx(i+7)=vx(i)    
         vy(i+7)=vy(i)  
         vz(i+7)=al_2  
 

 20      continue

         vx(14)=0.D0 
         vy(14)=0.D0   
         vz(14)=al_2+hprin  
         vv(14)=dsqrt(vx(14)**2+vy(14)**2+vz(14)**2)


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

          Do 30 i=8,12

         Vdirx(1)=(vx(14)-vx(i))/dsqrt((vx(14)-vx(i))**2+
     *    (vy(14)-vy(i))**2+(vz(14)-vz(i))**2)
         Vdiry(1)=(vy(14)-vy(i))/dsqrt((vx(14)-vx(i))**2+
     *    (vy(14)-vy(i))**2+(vz(14)-vz(i))**2)
         vdirz(1)=(vz(14)-vz(i))/dsqrt((vx(14)-vx(i))**2+
     *    (vy(14)-vy(i))**2+(vz(14)-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-1)=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-1)=(vdiry(1)*vdirz(2)-vdirz(1)*vdiry(2))/modu(i-1)
         DNy(i-1)=(vdirz(1)*vdirx(2)-vdirx(1)*vdirz(2))/modu(i-1)
         Dnz(i-1)=(vdirx(1)*vdiry(2)-vdiry(1)*vdirx(2))/modu(i-1)

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

 3       FORMAT('Nx=',F13.9,4x,'Ny=',F13.9,4x,'Nz=',F13.9)
         write(2,3)dnx(i-1),dny(i-1),dnz(i-1)
         
 30   continue

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


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


         modu(12)=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(12)=(vdiry(1)*vdirz(2)-vdirz(1)*vdiry(2))/modu(12)
         DNy(12)=(vdirz(1)*vdirx(2)-vdirx(1)*vdirz(2))/modu(12)
         Dnz(12)=(vdirx(1)*vdiry(2)-vdiry(1)*vdirx(2))/modu(12)

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

         write(2,3)dnx(12),dny(12),dnz(12)
 

         Do 40 i=1,6 

       DNx(i)= DNx(i+6)   
       DNY(i)= DNY(i+6) 
       DNZ(i)= -DNZ(i+6) 
        write(2,3)dnx(i),dny(i),dnz(i)

 40     continue
       

       DNx(13)= 0.D0   
       DNY(13)= 1.D0 
       DNZ(13)= 0.D0 
         write(2,3)dnx(13),dny(13),dnz(13)

       DNx(16)= 0.D0   
       DNY(16)= -1.D0 
       DNZ(16)= 0.D0
         write(2,3)dnx(16),dny(16),dnz(16)

 

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


         Vdirx(2)=(vx(3)-vx(10))/dsqrt((vx(3)-vx(10))**2+
     *    (vy(3)-vy(10))**2+(vz(3)-vz(10))**2)
         Vdiry(2)=(vy(3)-vy(10))/dsqrt((vx(3)-vx(10))**2+
     *    (vy(3)-vy(10))**2+(vz(3)-vz(10))**2)
         vdirz(2)=(vz(3)-vz(10))/dsqrt((vx(3)-vx(10))**2+
     *    (vy(3)-vy(10))**2+(vz(3)-vz(10))**2)


         modu(14)=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(14)=(vdiry(1)*vdirz(2)-vdirz(1)*vdiry(2))/modu(14)
         DNy(14)=(vdirz(1)*vdirx(2)-vdirx(1)*vdirz(2))/modu(14)
         Dnz(14)=(vdirx(1)*vdiry(2)-vdiry(1)*vdirx(2))/modu(14)

         Dnn(14)=Dsqrt(dnx(14)**2+dny(14)**2+dnz(14)**2)
         write(2,3)dnx(14),dny(14),dnz(14)
                  

       DNx(15)= DNx(14)   
       DNY(15)= -DNY(14) 
       DNZ(15)= DNZ(14) 

         Dnn(15)=Dsqrt(dnx(15)**2+dny(15)**2+dnz(15)**2)
         write(2,3)dnx(15),dny(15),dnz(15)
                  
       DNx(17)= -DNx(14)   
       DNY(17)= -DNY(14) 
       DNZ(17)=  DNZ(14) 

         Dnn(17)=Dsqrt(dnx(17)**2+dny(17)**2+dnz(17)**2)
         write(2,3)dnx(17),dny(17),dnz(17)
                  
       DNx(18)= -DNx(14)   
       DNY(18)= DNY(14) 
       DNZ(18)= DNZ(14) 
          Dnn(18)=Dsqrt(dnx(18)**2+dny(18)**2+dnz(18)**2)
         write(2,3)dnx(18),dny(18),dnz(18)
                  
        
**********************************************************************
        !Buiding up the prism

         NX=INT(r)  
         NY=INT(vy(2))
        apot=r*Dcos(pi/6.d0) 
         NZ_inf=INT(hprin)  
         NZ_sup=int(al_2+hprin)
          NAT=0
      
      DO 80 IZ=-NZ_inf,NZ_sup

         z=real(iz)
      DO 90 IY=-NY,NY
         Y=REAL(IY)

         DO  IX=-NX,NX
            
        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(6))*dnx(6)+(y-vy(6))*DNy(6)+(z-vz(6))*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(12))*dnx(11)+(y-vy(12))*DNy(11)+(z-vz(12))*dnz(11).le.0)
     *.and.
     * ((x-vx(13))*dnx(12)+(y-vy(13))*DNy(12)+(z-vz(13))*dnz(12).le.0)
     * .and.
     * ((x-vx(10))*dnx(13)+(y-vy(10))*DNy(13)+(z-vz(10))*dnz(13).le.0)
     * .and.
     * ((x-vx(11))*dnx(14)+(y-vy(11))*DNy(14)+(z-vz(11))*dnz(14).le.0)
     * .and.
     * ((x-vx(12))*dnx(15)+(y-vy(12))*DNy(15)+(z-vz(12))*dnz(15).le.0)
     * .and.
     * ((x-vx(13))*dnx(16)+(y-vy(13))*DNy(16)+(z-vz(13))*dnz(16).le.0)
     * .and.
     * ((x-vx(8))*dnx(17)+(y-vy(8))*DNy(17)+(z-vz(8))*dnz(17).le.0)
     * .and.
     * ((x-vx(9))*dnx(18)+(y-vy(9))*DNy(18)+(z-vz(9))*dnz(18).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(6))*dnx(6)+(y-vy(6))*DNy(6)+(z-vz(6))*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(12))*dnx(11)+(y-vy(12))*DNy(11)+(z-vz(12))*dnz(11).gt.0)
     * .or.
     * ((x-vx(13))*dnx(12)+(y-vy(13))*DNy(12)+(z-vz(13))*dnz(12).gt.0)
     *  .or.
     * ((x-vx(10))*dnx(13)+(y-vy(10))*DNy(13)+(z-vz(10))*dnz(13).gt.0)
     * .or.
     * ((x-vx(11))*dnx(14)+(y-vy(11))*DNy(14)+(z-vz(11))*dnz(14).gt.0)
     * .or.
     * ((x-vx(12))*dnx(15)+(y-vy(12))*DNy(15)+(z-vz(12))*dnz(15).gt.0)
     *.or.
     * ((x-vx(13))*dnx(16)+(y-vy(13))*DNy(16)+(z-vz(13))*dnz(16).gt.0)
     * .or.
     * ((x-vx(8))*dnx(17)+(y-vy(8))*DNy(17)+(z-vz(8))*dnz(17).gt.0).or.
     * ((x-vx(9))*dnx(18)+(y-vy(9))*DNy(18)+(z-vz(9))*dnz(18).gt.0))
     *   goto 95
      
          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

 95      end do
         
 90       CONTINUE
 80       CONTINUE
         

 50       STOP
    
        END

     
     
    
         
      

    
