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 ROUNDED_DECAHEDRON

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 pentagonal base is parallel to the xy plane
       ! and  one edge of the base is parallel to the x axis with
       ! the opposite vertex  parallel to the y axis.

       !TRUNCATION: Only done on the 5 pentagonal base vertex.

      PARAMETER(mxnat=9999999)
      PARAMETER(Nver=7)        
      PARAMETER(NPUNT=5)       
      PARAMETER(nnor=10)       
      PARAMETER(m=2)
      REAL*8 AL,A,B,PI,X,Y,Z,alphamed,betha,modu,apot,TRUN
      REAL*8 D,H,hprin,GAMMA,vdirx,vdiry,vdirz,vx,vy,vz
      REAL*8 P_X,P_Y,P_Z,VC_X,VC_Y,VC_Z,trun_factor,truncation
      REAL*8 dnx,dny,dnz,vv,norma,vvdir,dnn,VE_X,VE_Y,VE_Z
      REAL*8 dnx_E,dny_E,dnz_E
      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) 
      DIMENSION VE_X(NPUNT),VE_Y(NPUNT),VE_Z(NPUNT)
      DIMENSION DNX_E(nnor),DNY_E(nnor),DNZ_E(Nnor)
      DIMENSION VC_X(nnor),VC_Y(nnor),VC_Z(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(*,*)'The maximum value 0.5 gives'
        Write(*,*)' a shape similar to the regular decahedron.'
        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(*,*)'Decahedron with a truncated value= 1/',truncation  

          end if
        
        Write(*,*)'Give me a length to generate N entities to build up'
        write(*,*)' the rounded decahedron. For example:'
        write(*,*)'If you  choose 1/8 as the truncated value and:'
        write(*,*)'length = 10 ==> N = 606'
        write(*,*)'length = 20 ==> N = 4,829'
        write(*,*)'length = 30 ==> N =16,302'
        read(*,*)Al  


       if (al.le.0) then
           write(*,*)'error in length value'
           stop'end of program'
          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='rounded_normals.dat')

c****************************************************************************
         ! Coordinates of each  initial regular decahedron vertex
       
         vx(1)=-al/2.d0-d    
         vy(1)=dsqrt(apot**2-d**2-d*al)   
         vz(1)=0.d0  
         vv(1)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2)

         vx(2)=0.d0  
         vy(2)=b+h-apot 
         vz(2)=0.d0  
         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)=0.d0  
         vv(3)=dsqrt(vx(3)**2+vy(3)**2+vz(3)**2)

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

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

         vx(6)=0.d0  
         vy(6)=0.d0
         vz(6)=hprin   
         vv(6)=dsqrt(vx(6)**2+vy(6)**2+vz(6)**2)

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


C******************************************************************
C  10 vertex coordinates to generate the truncated decahedron


         TRUN=AL*trun_factor  

         VC_X(1)=-TRUN*DSIN(ALPHAMED)    
         VC_Y(1)=VY(2)-TRUN*DCOS(ALPHAMED)
         VC_Z(1)=0.d0 

         VC_X(2)=TRUN*DSIN(ALPHAMED)        
         VC_Y(2)=VY(2)-TRUN*DCOS(ALPHAMED)
         VC_Z(2)=0.d0 

         VC_X(3)=(AL-TRUN)*DSIN(ALPHAMED)      
         VC_Y(3)=VY(2)-(AL-TRUN)*DCOS(ALPHAMED) 
         VC_Z(3)=0.d0 

         VC_X(4)=VX(3)-TRUN*DCOS(BETHA)    
         VC_Y(4)=VY(3)-TRUN*DSIN(BETHA)
         VC_Z(4)=0.d0 

         VC_X(5)=VX(4)+TRUN*DCOS(BETHA)      
         VC_Y(5)=VY(4)+TRUN*DSIN(BETHA)
         VC_Z(5)=0.d0 

         VC_X(6)= VX(4)-TRUN    
         VC_Y(6)=VY(4)
         VC_Z(6)=0.d0 

         VC_X(7)=VX(5)+TRUN      
         VC_Y(7)=VY(4)
         VC_Z(7)=0.d0 

         VC_X(8)=VX(5)-TRUN*DCOS(BETHA)   
         VC_Y(8)=VY(4)+TRUN*DSIN(BETHA)
         VC_Z(8)=0.d0 

         vC_X(9)= -(VX(3)-TRUN*DCOS(BETHA))     
         VC_Y(9)=VY(3)-TRUN*DSIN(BETHA)
         VC_Z(9)=0.d0 

          VC_X(10)=-(AL-TRUN)*DSIN(ALPHAMED)       
         VC_Y(10)=VY(2)-(AL-TRUN)*DCOS(ALPHAMED) 
         VC_Z(10)=0.d0 

C******************************************************************
         !  NORMAL VECTOR COMPONENTS TO EACH FACE OF THE 
         !  REGULAR DECAHEDRON 

          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+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(7)-vx(i))/dsqrt((vx(7)-vx(i))**2+
     *    (vy(7)-vy(i))**2+(vz(7)-vz(i))**2)
         Vdiry(2)=(vy(7)-vy(i))/dsqrt((vx(7)-vx(i))**2+
     *    (vy(7)-vy(i))**2+(vz(7)-vz(i))**2)
         vdirz(2)=(vz(7)-vz(i))/dsqrt((vx(7)-vx(i))**2+
     *    (vy(7)-vy(i))**2+(vz(7)-vz(i))**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

C***********************************************************************
         !Here comes the truncation
         !The 5 truncation planes are parallel to the Z axis
         
         !Components of 5 Normal vectors to build up the truncated decahedron.
         

         Do 65 i=1,5    
            
            P_X=VC_X(2*i-1)
            P_Y=VC_Y(2*i-1)
            P_Z=2.d0 

         Vdirx(1)=(P_X-VC_X(2*i-1))/dsqrt((P_X-VC_X(2*i-1))**2+
     *    (P_Y-VC_Y(2*i-1))**2+(P_Z-VC_Z(2*i-1))**2)
         Vdiry(1)=(P_Y-VC_Y(2*i-1))/dsqrt((P_X-VC_X(2*i-1))**2+
     *    (P_Y-VC_Y(2*i-1))**2+(P_Z-VC_Z(2*i-1))**2)
         vdirz(1)=(P_Z-VC_Z(2*i-1))/dsqrt((P_X-VC_X(2*i-1))**2+
     *    (P_Y-VC_Y(2*i-1))**2+(P_Z-VC_Z(2*i-1))**2)


       Vdirx(2)=(VC_X(2*i)-VC_X(2*i-1))/dsqrt((VC_X(2*i)-VC_X(2*i-1))**2
     *   + (VC_Y(2*i)-VC_Y(2*i-1))**2+(VC_Z(2*i)-VC_Z(2*i-1))**2)
       Vdiry(2)=(VC_Y(2*i)-VC_Y(2*i-1))/dsqrt((VC_X(2*i)-VC_X(2*i-1))**2
     *   + (VC_Y(2*i)-VC_Y(2*i-1))**2+(VC_Z(2*i)-VC_Z(2*i-1))**2)
       vdirz(2)=(VC_Z(2*i)-VC_Z(2*i-1))/dsqrt((VC_X(2*i)-VC_X(2*i-1))**2
     *   + (VC_Y(2*i)-VC_Y(2*i-1))**2+(VC_Z(2*i)-VC_Z(2*i-1))**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_E(i)=(vdiry(2)*vdirz(1)-vdirz(2)*vdiry(1))/modu(i)
         DNy_E(i)=(vdirz(2)*vdirx(1)-vdirx(2)*vdirz(1))/modu(i)
         Dnz_E(i)=(vdirx(2)*vdiry(1)-vdiry(2)*vdirx(1))/modu(i)

         Dnn(i)=Dsqrt(dnx_E(i)**2+dny_E(i)**2+dnz_E(i)**2)

         write(2,3)dnx_E(i),dny_E(i),dnz_E(i)
         
  65      continue
                     
C**********************************************************************
         ! Building up the rounded decahedron.

         NB=INT(apot)  
         NA=INT(vx(3))
         NZZ=INT(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(5))*dnx(6)+(y-vy(5))*DNy(6)+(z-vz(5))*dnz(6).le.0).and.
     * ((x-vx(1))*dnx(7)+(y-vy(1))*DNy(7)+(z-vz(1))*dnz(7).le.0).and.
     * ((x-vx(2))*dnx(8)+(y-vy(2))*DNy(8)+(z-vz(2))*dnz(8).le.0).and.
     * ((x-vx(3))*dnx(9)+(y-vy(3))*DNy(9)+(z-vz(3))*dnz(9).le.0).and.
     * ((x-vx(4))*dnx(10)+(y-vy(4))*DNy(10)+(z-vz(4))*dnz(10).le.0).and.
     * ((x-VC_X(1))*dnx_E(1)+(y-VC_Y(1))*DNy_E(1)
     * +(z-VC_Z(1))*dnz_E(1).gt.0).and.
     * ((x-VC_X(3))*dnx_E(2)+(y-VC_Y(3))*DNy_E(2)
     * +(z-VC_Z(3))*dnz_E(2).gt.0).and.
     * ((x-VC_X(5))*dnx_E(3)+(y-VC_Y(5))*DNy_E(3)
     * +(z-VC_Z(5))*dnz_E(3).gt.0).and.
     * ((x-VC_X(7))*dnx_E(4)+(y-VC_Y(7))*DNy_E(4)
     * +(z-VC_Z(7))*dnz_E(4).gt.0).and.
     * ((x-VC_X(9))*dnx_E(5)+(y-VC_Y(9))*DNy_E(5)
     * +(z-VC_Z(9))*dnz_E(5).gt.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(5))*dnx(6)+(y-vy(5))*DNy(6)+(z-vz(5))*dnz(6).gt.0).or.
     * ((x-vx(1))*dnx(7)+(y-vy(1))*DNy(7)+(z-vz(1))*dnz(7).gt.0).or.
     * ((x-vx(2))*dnx(8)+(y-vy(2))*DNy(8)+(z-vz(2))*dnz(8).gt.0).or.
     * ((x-vx(3))*dnx(9)+(y-vy(3))*DNy(9)+(z-vz(3))*dnz(9).gt.0).or.
     * ((x-vx(4))*dnx(10)+(y-vy(4))*DNy(10)+(z-vz(4))*dnz(10).gt.0).or.
     * ((x-VC_X(1))*dnx_E(1)+(y-VC_Y(1))*DNy_E(1)
     *    +(z-VC_Z(1))*dnz_E(1).le.0).or.
     * ((x-VC_X(3))*dnx_E(2)+(y-VC_Y(3))*DNy_E(2)
     * +(z-VC_Z(3))*dnz_E(2).le.0).or.  
     * ((x-VC_X(5))*dnx_E(3)+(y-VC_Y(5))*DNy_E(3)
     * +(z-VC_Z(5))*dnz_E(3).le.0).or.
     * ((x-VC_X(7))*dnx_E(4)+(y-VC_Y(7))*DNy_E(4)
     * +(z-VC_Z(7))*dnz_E(4).le.0).or.
     * ((x-VC_X(9))*dnx_E(5)+(y-VC_Y(9))*DNy_E(5)
     * +(z-VC_Z(9))*dnz_E(5).le.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

     
     
    
         
      

    
