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 TETRAHEDRAL_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 
      ! The axis that join two opposite vertex of the square cross section
      ! is parallel to the x, y axis, respectively.

      parameter(mxnat=950000) 
      PARAMETER(Nver=10)
      parameter(nnor=12)
      parameter(m=2)
      IMPLICIT REAL*8 (A-H,O-Q,T-Z) 
      REAL*8 AL,PI,X,Y,modu,dRAZ,dl_2
      REAL*8 vdirx,vdiry,vdirz,vx,vy,vz
      real*8 dnx,dny,dnz,vv,vvdir,dnn
      INTEGER IY,IX,IZ,IXYZ,NAT,mxnat,NZZ,nver,nnor,m,nlong,nlongz
      DIMENSION IXYZ(mxNAT,3)  
      DIMENSION vdirx(m),vdiry(m),vdirz(m),dnn(Nnor)
      dimension vv(nver),norma(nnor),vvdir(m)       
      DIMENSION VX(nver),VY(Nver),VZ(Nver),DNX(nnor),DNY(nnor),DNZ(Nnor) 
   

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

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

      write(*,*)'Give a length to generate N dipole entities.'
      write(*,*)'For example, if aspect ratio=2'
      write(*,*)'length = 10 then N = 841 '
      write(*,*)'length = 20 then N = 5,981'
      write(*,*)'length = 30 then N = 19,421' 
      read(*,*)al

      If(al.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='tetraprism_normales.dat')
 

          
          dl_2=al*draz-al     !length of the rectangular faces
        
          Nlong=int(al)
 
          Nlongz=int(al+dl_2) 
       
C******************************************************************
         !Vertex coordinates of the prism

         vx(1)=al       
         vy(1)=al/2.d0     
         vz(1)= al/2.d0      
         vv(1)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2) 

         vx(2)=al/2.d0       
         vy(2)=al
         vz(2)=al/2.d0     
         vv(2)=dsqrt(vx(2)**2+vy(2)**2+vz(2)**2)

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

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

         vx(5)=al      !
         vy(5)=al/2.d0     
         vz(5)=al/2.d0 +dl_2    
         vv(5)=dsqrt(vx(5)**2+vy(5)**2+vz(5)**2)

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

         vx(7)=0.d0!
         vy(7)=al/2.d0     
         vz(7)=al/2.d0+dl_2     
         vv(7)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(8)=al/2.d0     
         vy(8)=0.d0
         vz(8)=al/2.d0+dl_2       !
         vv(8)=dsqrt(vx(8)**2+vy(8)**2+vz(8)**2)
         
         vx(9)=al/2.d0     
         vy(9)=al/2.d0     
         vz(9)=al+dl_2
         vv(9)=dsqrt(vx(9)**2+vy(9)**2+vz(9)**2)

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

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

         dnx(1)=1.d0
         dny(1)=1.d0
         dnz(1)=1.d0

          dnx(2)=-1.d0
         dny(2)=1.d0
         dnz(2)=1.d0

          dnx(3)=-1.d0
         dny(3)=-1.d0
         dnz(3)=1.d0

          dnx(4)=1.d0
         dny(4)=-1.d0
         dnz(4)=1.d0

          dnx(5)=1.d0
         dny(5)=1.d0
         dnz(5)=-1.d0

          dnx(6)=-1.d0
         dny(6)=1.d0
         dnz(6)=-1.d0

          dnx(7)=-1.d0
         dny(7)=-1.d0
         dnz(7)=-1.d0

          dnx(8)=1.d0
         dny(8)=-1.d0
         dnz(8)=-1.d0

          dnx(9)=1.d0
         dny(9)=1.d0
         dnz(9)=0.d0

          dnx(10)=-1.d0
         dny(10)=1.d0
         dnz(10)=0.d0

          dnx(11)=-1.d0
         dny(11)=-1.d0
         dnz(11)=0.d0


          dnx(12)=1.d0
         dny(12)=-1.d0
         dnz(12)=0.d0        

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

       Do 45 i=1,12

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


C****************************************************************
        !Building up the tetrahedral prism

      NAT=0      

      DO 80 IZ=0,nlongz

         z=real(iz)
      DO 10 IY=0,nlong
         Y=REAL(IY)

         DO  IX=0,nlong
            
        NAT=NAT+1
               X=REAL(IX)

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

          IXYZ(NAT,3)=int(z)
          IXYZ(NAT,2)=INT(y)  

          
          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

