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 TRUNCATEDCUBE

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 Silver Nanoparticles", 
C	   A. L. Gonzalez and C. Noguez, Physica Status Soildi (c) 4, 4118 - 4126 (2007).
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 SQUARE FACES ARE PARALLEL TO THE X,Y,Z AXIS.

      PARAMETER(mxnat=950000) 
      PARAMETER(Nver=24)      
      PARAMETER(nnor=24)                                   
      PARAMETER(m=2)
      REAL*8 AL,PI,X,Y,Z,modu,trun,trun_factor,truncation
      REAL*8 vdirx,vdiry,vdirz,vx,vy,vz
      REAL*8 dnx,dny,dnz,vv,norma,vvdir,dnn
      INTEGER IY,IX,IZ,IXYZ,NAT,mxnat,NZZ,nver,nnor,m,nlong,i
      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) 
   
c****************************************************************************


        write(*,*)'The value of the truncation factor needs to be'
        write(*,*)'a  number greater than 0 and'
        write(*,*)' less than or equal to 0.5' 
        Write(*,*)'For example:'
        Write(*,*)'If Truncation Factor  = 1/2 then write 0.5'
        Write(*,*)'If Truncation Factor = 1/4 then write 0.25 '
        Write(*,*)'If Truncation Factor = 1/8 then write 0.125'
        read(*,*) trun_factor


        truncation=1./trun_factor 
        If(trun_factor.le.0.or.trun_factor.gt.0.5)then
           write(*,*)'error in truncated value'
           stop'end of program'
        
           else if(trun_factor.gt.0.or.trun_factor.le.0.5)then
           write(*,*)'Cube with a truncated value= 1/',truncation  

          end if
        
        Write(*,*)'Give me a length to generate N entities to build up'
        write(*,*)' the truncated cube. For example:'
        write(*,*)'If you  choose 1/2 as the truncated value and:'
        write(*,*)'length = 10 ==> N = 1,051'
        write(*,*)'length = 20 ==> N = 7,501'
        write(*,*)'length = 30 ==> N =24,351'

        read(*,*)al
  

      PI = 4.d0*datan(1.d0)
      trun=al*trun_factor                


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

      OPEN(UNIT=2,FILE='truncube_normals.dat')
     
      
       if (al.le.0) then
           write(*,*)'error in length value'
           stop'end of program'
          end if
          Nlong=int(al) !number of atoms in the initial cube

c*********************************************************************
         !X,Y,Z coordinates of each vertex of the truncated cube. 

         vx(1)=0.0d0      
         vy(1)=trun     
         vz(1)=0.d0  
         vv(1)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2)

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

         vx(4)=AL  
         vy(4)=0.d0
         vz(4)=trun
         vv(4)=dsqrt(vx(4)**2+vy(4)**2+vz(4)**2)

         vx(5)=AL-TRUN  
         vy(5)=0.d0
         vz(5)=0.d0
         vv(5)=dsqrt(vx(5)**2+vy(5)**2+vz(5)**2)
        
         vx(6)=AL  
         vy(6)=trun
         vz(6)=0.d0
         vv(6)=dsqrt(vx(6)**2+vy(6)**2+vz(6)**2)

         vx(7)=AL-TRUN  
         vy(7)=0.d0
         vz(7)=AL  
         vv(7)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(8)=AL   
         vy(8)=0.d0
         vz(8)=AL-TRUN 
         vv(8)=dsqrt(vx(8)**2+vy(8)**2+vz(8)**2)
         
         vx(9)=AL    
         vy(9)=trun
         vz(9)=AL     
         vv(9)=dsqrt(vx(9)**2+vy(9)**2+vz(9)**2)

         vx(10)=0.d0
         vy(10)=trun
         vz(10)=AL     
         vv(10)=dsqrt(vx(10)**2+vy(10)**2+vz(10)**2)

         vx(11)=0.d0
         vy(11)=0.d0
         vz(11)=AL-TRUN   
         vv(11)=dsqrt(vx(11)**2+vy(11)**2+vz(11)**2)

         vx(12)=trun
         vy(12)=0.d0
         vz(12)=AL 
         vv(12)=dsqrt(vx(12)**2+vy(12)**2+vz(12)**2)

         vx(13)=AL      
         vy(13)=AL-TRUN  
         vz(13)=al   
         vv(13)=dsqrt(vx(13)**2+vy(13)**2+vz(13)**2)

         vx(14)=AL  
         vy(14)=AL    
         vz(14)=AL-TRUN 
         vv(14)=dsqrt(vx(14)**2+vy(14)**2+vz(14)**2)

         vx(15)=AL-TRUN
         vy(15)=AL  
         vz(15)= AL 
         vv(15)=dsqrt(vx(15)**2+vy(15)**2+vz(15)**2)

         vx(16)=trun
         vy(16)=AL     
         vz(16)=AL    
         vv(16)=dsqrt(vx(16)**2+vy(16)**2+vz(16)**2)

         vx(17)=0.d0
         vy(17)=AL    
         vz(17)=AL-TRUN  
         vv(17)=dsqrt(vx(17)**2+vy(17)**2+vz(17)**2)

         vx(18)=0.d0
         vy(18)=AL-TRUN 
         vz(18)=al 
         vv(18)=dsqrt(vx(18)**2+vy(18)**2+vz(18)**2)

         vx(19)=0.d0
         vy(19)=al  
         vz(19)=trun  
         vv(19)=dsqrt(vx(19)**2+vy(19)**2+vz(19)**2)

         vx(20)=trun
         vy(20)=al   
         vz(20)=0.d0
         vv(20)=dsqrt(vx(20)**2+vy(20)**2+vz(20)**2)

         vx(21)=0.d0
         vy(21)=AL-TRUN  
         vz(21)=0.d0
         vv(21)=dsqrt(vx(21)**2+vy(21)**2+vz(21)**2)

         vx(22)=al    
         vy(22)=AL-TRUN 
         vz(22)=0.d0
         vv(22)=dsqrt(vx(22)**2+vy(22)**2+vz(22)**2)

         vx(23)=AL-TRUN
         vy(23)=al   
         vz(23)=0.d0
         vv(23)=dsqrt(vx(23)**2+vy(23)**2+vz(23)**2)

         vx(24)=al    
         vy(24)=al    
         vz(24)=trun
         vv(24)=dsqrt(vx(24)**2+vy(24)**2+vz(24)**2)


C************************************************************************
          ! NORMAL VECTOR COMPONENTS TO EACH FACE OF THE TRUNCATED CUBE


         
         Do 45 i=1,nver,3
 
         Vdirx(1)=(vx(i+2)-vx(i+1))/dsqrt((vx(i+2)-vx(i+1))**2+
     *    (vy(i+2)-vy(i+1))**2+(vz(i+2)-vz(i+1))**2)
         Vdiry(1)=(vy(i+2)-vy(i+1))/dsqrt((vx(i+2)-vx(i+1))**2+
     *    (vy(i+2)-vy(i+1))**2+(vz(i+2)-vz(i+1))**2)
         vdirz(1)=(vz(i+2)-vz(i+1))/dsqrt((vx(i+2)-vx(i+1))**2+
     *    (vy(i+2)-vy(i+1))**2+(vz(i+2)-vz(i+1))**2)



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


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

         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

 
C************************************************************************      
       ! BUILDING UP  THE TRUNCATED CUBE 

      NAT=0      

      DO 80 IZ=0,nlong

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

         DO  IX=0,nlong
            
        NAT=NAT+1
               X=REAL(IX)
        !Testing if the point (x,y,z) is inside or on one of the  planes that gives 
        !the shape of the truncated cube

      If (((x-vx(1))*dnx(1)+(y-vy(1))*DNy(1)+
     * (z-vz(1))*dnz(1).le.0).and.((x-vx(4))*dnx(4)+(y-vy(4))*DNy(4)
     * +(z-vz(4))*dnz(4).le.0).and.
     * ((x-vx(7))*dnx(7)+(y-vy(7))*DNy(7)+(z-vz(7))*dnz(7).le.0).and.
     * ((x-vx(10))*dnx(10)+(y-vy(10))*DNy(10)+(z-vz(10))*dnz(10).le.0)
     * .and.
     * ((x-vx(13))*dnx(13)+(y-vy(13))*DNy(13)+(z-vz(13))*dnz(13).le.0)
     * .and.
     * ((x-vx(16))*dnx(16)+(y-vy(16))*DNy(16)+(z-vz(16))*dnz(16).le.0)
     * .and.
     * ((x-vx(19))*dnx(19)+(y-vy(19))*DNy(19)+(z-vz(19))*dnz(19).le.0)
     * .and.
     * ((x-vx(22))*dnx(22)+(y-vy(22))*DNy(22)+(z-vz(22))*dnz(22).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(4))*dnx(4)+(y-vy(4))*DNy(4)+(z-vz(4))*dnz(4).gt.0)
     *   .or.
     *   ((x-vx(7))*dnx(7)+(y-vy(7))*DNy(7)+(z-vz(7))*dnz(7).gt.0)
     *    .or.
     * ((x-vx(10))*dnx(10)+(y-vy(10))*DNy(10)+(z-vz(10))*dnz(10).gt.0)
     *    .or.
     * ((x-vx(13))*dnx(13)+(y-vy(13))*DNy(13)+(z-vz(13))*dnz(13).gt.0)
     *    .or.
     * ((x-vx(16))*dnx(16)+(y-vy(16))*DNy(16)+(z-vz(16))*dnz(16).gt.0)
     *    .or.
     * ((x-vx(19))*dnx(19)+(y-vy(19))*DNy(19)+(z-vz(19))*dnz(19).gt.0)
     *    .or.
     * ((x-vx(22))*dnx(22)+(y-vy(22))*DNy(22)+(z-vz(22))*dnz(22).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

 
