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 ICOSAHEDRON

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 Absorbance of Col loidal Suspensions of Silver Polyhedral Nanoparticles", 
C      A. L. Gonzalez, C. Noguez, G. P. Ortiz, G. Rodriguez-Gattorno, J. Phys. Chem. B. 
C      109 17512 - 17517 (2005).
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

      PARAMETER(mxnat=950000) 
      PARAMETER(Nver=12)
      PARAMETER(nnor=20)
      PARAMETER(m=2)
      REAL*8 vdirx,vdiry,vdirz,vx,vy,vz,modu,x,y,z
      REAL*8 dnx,dny,dnz,vv,norma,vvdir,dnn,factor
      INTEGER N,NA,NB,IY,IX,IZ,IXYZ,NAT,mxnat,NZZ,nver,nnor,m
      INTEGER Natom,i
      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 me a factor to generate N entities to build up'
        write(*,*)'the icosahedron. For example:'
	write(*,*)'     factor = 10  ==> N =   4093'
	write(*,*)'     factor = 20  ==> N =  32953'
        write(*,*)'     factor = 30  ==> N = 111231'
        write(*,*)'     factor = 40  ==> N = 263675'

      read(*,*)factor  
 
      OPEN(UNIT=1,FILE='coordinates.dat') 
      OPEN(UNIT=2,FILE='ico_normals.dat')   
      OPEN(UNIT=3,FILE='ico_vertex.dat')

          !X,Y,Z COORDINATES OF EACH VERTEX OF THE ICOSAHEDRON 
       
         vx(1)=0.0d0*factor   
         vy(1)=0.0d0*factor  
         vz(1)=1.1755705045849463d0*factor  
         vv(1)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2)

         vx(2)=1.0514622242382672d0*factor
         vy(2)=0.d0*factor 
         vz(2)=0.5257311121191333d0*factor  
         vv(2)=dsqrt(vx(2)**2+vy(2)**2+vz(2)**2)

         vx(3)=0.3249196962329063d0*factor 
         vy(3)=-1.0d0*factor  
         vz(3)= 0.5257311121191333d0*factor 
         vv(3)=dsqrt(vx(3)**2+vy(3)**2+vz(3)**2)

         vx(4)=-0.8506508083520399d0*factor  
         vy(4)=-0.6180339887498949d0*factor
         vz(4)=0.5257311121191333d0*factor
         vv(4)=dsqrt(vx(4)**2+vy(4)**2+vz(4)**2)

         vx(5)=-0.8506508083520399d0*factor 
         vy(5)=0.6180339887498949d0*factor
         vz(5)=0.5257311121191333d0*factor 
         vv(5)=dsqrt(vx(5)**2+vy(5)**2+vz(5)**2)

         vx(6)=0.3249196962329063d0*factor  
         vy(6)=1.0d0*factor
         vz(6)=0.5257311121191333d0*factor 
         vv(6)=dsqrt(vx(6)**2+vy(6)**2+vz(6)**2)

         vx(7)=0.8506508083520399d0*factor 
         vy(7)=0.6180339887498949d0*factor
         vz(7)=-0.5257311121191333d0*factor
         vv(7)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(8)=0.8506508083520399d0*factor  
         vy(8)=-0.6180339887498949d0*factor 
         vz(8)=-0.5257311121191333d0*factor
         vv(8)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(9)=-0.3249196962329063d0*factor  
         vy(9)=-1.0d0*factor 
         vz(9)=-0.5257311121191333d0*factor
         vv(9)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(10)=-1.0514622242382672d0*factor  
         vy(10)= 0.0d0*factor
         vz(10)=-0.5257311121191333d0*factor
         vv(10)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(11)=-0.3249196962329063d0*factor 
         vy(11)=1.0d0*factor 
         vz(11)=-0.5257311121191333d0*factor
         vv(11)=dsqrt(vx(7)**2+vy(7)**2+vz(7)**2)

         vx(12)=0.0D0*factor    
         vy(12)=0.0D0*factor
         vz(12)=-1.1755705045849463*factor  
         vv(12)=dsqrt(vx(1)**2+vy(1)**2+vz(1)**2) 

 4       FORMAT(I3,2x,'vx=',F14.8,2x,'vy=',F14.8,2x,'vz=',F14.8)
        Do 5 i=1,12
 
         write(3,4)i,vx(i),vy(i),vz(i)

 5       continue


         ! NORMAL VECTOR COMPONENTS TO EACH FACE OF THE ICOSAHEDRON

          Do 45 i=1,4 

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


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

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

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

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

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

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

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

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

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

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


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

         Dnn(6)=Dsqrt(dnx(6)**2+dny(6)**2+dnz(6)**2)
         write(2,3)6,dnx(6),dny(6),dnz(6)
       
           Do 55 i=2,5 

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

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


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

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

 55      continue
       
           Do 65 i=2,5

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

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

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

         Dnn(2*i+3)=Dsqrt(dnx(2*i+3)**2+dny(2*i+3)**2+dnz(2*i+3)**2)
         write(2,3)2*i+3,dnx(2*i+3),dny(2*i+3),dnz(2*i+3)
        
 65      continue

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

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


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

         Dnn(15)=Dsqrt(dnx(15)**2+dny(15)**2+dnz(15)**2)
         
         write(2,3)15,dnx(15),dny(15),dnz(15)
 
          Do 75 i=7,10 

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

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

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

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

 75      continue

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

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


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

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

        ! BUILDING UP  THE ICOSAHEDRON

         NA=INT(vx(2))  
         NB=INT(vy(11))
         NZZ=INT(vz(1))  

          NAT=0

      DO 80 IZ=-NZZ,NZZ

         z=real(iz)
      DO 10 IY=-NB,NB
         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(1))*dnx(2)+(y-vy(1))*DNy(2)+(z-vz(1))*dnz(2).le.0).and.
     *    ((x-vx(1))*dnx(3)+(y-vy(1))*DNy(3)+(z-vz(1))*dnz(3).le.0).and.
     *    ((x-vx(1))*dnx(4)+(y-vy(1))*DNy(4)+(z-vz(1))*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(6))*dnx(6)+(y-vy(6))*DNy(6)+(z-vz(6))*dnz(6).le.0).and. 
     *    ((x-vx(7))*dnx(7)+(y-vy(7))*DNy(7)+(z-vz(7))*dnz(7).le.0).and.     
     *    ((x-vx(8))*dnx(8)+(y-vy(8))*DNy(8)+(z-vz(8))*dnz(8).le.0).and.
     *    ((x-vx(9))*dnx(9)+(y-vy(9))*DNy(9)+(z-vz(9))*dnz(9).le.0).and.     
     *    ((x-vx(9))*dnx(10)+(y-vy(9))*DNy(10)+(z-vz(9))*dnz(10).le.0)
     *  .and.     
     *  ((x-vx(10))*dnx(11)+(y-vy(10))*DNy(11)+(z-vz(10))*dnz(11).le.0)
     *  .and.     
     *  ((x-vx(10))*dnx(12)+(y-vy(10))*DNy(12)+(z-vz(10))*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(11))*dnx(15)+(y-vy(11))*DNy(15)+(z-vz(11))*dnz(15).le.0)
     *  .and.     
     *  ((x-vx(12))*dnx(16)+(y-vy(12))*DNy(16)+(z-vz(12))*dnz(16).le.0)
     *  .and.
     *  ((x-vx(12))*dnx(17)+(y-vy(12))*DNy(17)+(z-vz(12))*dnz(17).le.0)
     *  .and.
     *  ((x-vx(12))*dnx(18)+(y-vy(12))*DNy(18)+(z-vz(12))*dnz(18).le.0)
     *  .and.
     *  ((x-vx(12))*dnx(19)+(y-vy(12))*DNy(19)+(z-vz(12))*dnz(19).le.0)
     *  .and.
     *  ((x-vx(12))*dnx(20)+(y-vy(12))*DNy(20)+(z-vz(12))*dnz(20).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(1))*dnx(2)+(y-vy(1))*DNy(2)+(z-vz(1))*dnz(2).gt.0).or.
     *    ((x-vx(1))*dnx(3)+(y-vy(1))*DNy(3)+(z-vz(1))*dnz(3).gt.0).or.
     *    ((x-vx(1))*dnx(4)+(y-vy(1))*DNy(4)+(z-vz(1))*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(6))*dnx(6)+(y-vy(6))*DNy(6)+(z-vz(6))*dnz(6).gt.0).or.
     *    ((x-vx(7))*dnx(7)+(y-vy(7))*DNy(7)+(z-vz(7))*dnz(7).gt.0).or.     
     *    ((x-vx(8))*dnx(8)+(y-vy(8))*DNy(8)+(z-vz(8))*dnz(8).gt.0).or.
     *    ((x-vx(9))*dnx(9)+(y-vy(9))*DNy(9)+(z-vz(9))*dnz(9).gt.0).or.     
     *    ((x-vx(9))*dnx(10)+(y-vy(9))*DNy(10)+(z-vz(9))*dnz(10).gt.0)
     *  .or.     
     *  ((x-vx(10))*dnx(11)+(y-vy(10))*DNy(11)+(z-vz(10))*dnz(11).gt.0)
     *  .or.     
     *  ((x-vx(10))*dnx(12)+(y-vy(10))*DNy(12)+(z-vz(10))*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(11))*dnx(15)+(y-vy(11))*DNy(15)+(z-vz(11))*dnz(15).gt.0)
     *  .or.     
     *  ((x-vx(12))*dnx(16)+(y-vy(12))*DNy(16)+(z-vz(12))*dnz(16).gt.0)
     *  .or.
     *  ((x-vx(12))*dnx(17)+(y-vy(12))*DNy(17)+(z-vz(12))*dnz(17).gt.0)
     *  .or.
     *  ((x-vx(12))*dnx(18)+(y-vy(12))*DNy(18)+(z-vz(12))*dnz(18).gt.0)
     *  .or.
     *  ((x-vx(12))*dnx(19)+(y-vy(12))*DNy(19)+(z-vz(12))*dnz(19).gt.0)
     *  .or.
     *  ((x-vx(12))*dnx(20)+(y-vy(12))*DNy(20)+(z-vz(12))*dnz(20).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

     
     
    
         
      

    
