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 OCTAHEDRON


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  3. "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: Each axis that pass through two oposite vertex
      ! is parallel to the x,y, and z axis, respectively.

      PARAMETER(mxnat=12500000) 
      PARAMETER(nnor=8)       
      INTEGER IX,IY,IZ,IXYZ,NAT,mxnat,Nlong   
      REAL*8 AL,X,Y,Z,DNX,DNY,DNZ,Sn
      DIMENSION IXYZ(mxnat,3)         
      DIMENSION DNX(nnor),DNY(nnor),DNZ(Nnor) 
   
        Write(*,*)'Give me a length to generate N entities to build up'
        write(*,*)' the octahedron. For example:'
        write(*,*)'Length = 10 ==> N =   231'
        write(*,*)'Length = 20 ==> N = 1,561'
        write(*,*)'Length = 30 ==> N = 4,991'

        read(*,*)al
           

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

    
       If (AL.le.0) goto 50
        
          
          Nlong=INT(AL) 
          write(*,*)Nlong

c****************************************************************************
          ! Normal vector components to each face of the octahedron

          DNX(1)=1.d0/sqrt(3.)
          DNY(1)=1.d0/sqrt(3.)
          DNZ(1)=1.d0/sqrt(3.)
        
          DNX(2)=-1.d0/sqrt(3.)
          DNY(2)=1.d0/sqrt(3.)
          DNZ(2)=1.d0/sqrt(3.)
         
          DNX(3)=-1.d0/sqrt(3.)
          DNY(3)=-1.d0/sqrt(3.)
          DNZ(3)=1.d0/sqrt(3.)
        
          DNX(4)=1.d0/sqrt(3.)
          DNY(4)=-1.d0/sqrt(3.)
          DNZ(4)=1.d0/sqrt(3.)
        
          DNX(5)=1.d0/sqrt(3.)
          DNY(5)=1.d0/sqrt(3.)
          DNZ(5)=-1.d0/sqrt(3.)
        
          DNX(6)=-1.d0/sqrt(3.)
          DNY(6)=1.d0/sqrt(3.)
          DNZ(6)=-1.d0/sqrt(3.)
        
          DNX(7)=-1.d0/sqrt(3.)
          DNY(7)=-1.d0/sqrt(3.)
          DNZ(7)=-1.d0/sqrt(3.)
     
          DNX(8)=1.d0/sqrt(3.)
          DNY(8)=-1.d0/sqrt(3.)
          DNZ(8)=-1.d0/sqrt(3.)
        
      NAT=0      

C************************************************************************      
       !Building up the octahedron

      !Sn vector in the direction (111)/sqrt(3),distance from the CM
      ! to a plane

      Sn=(Nlong+1)/6.d0   

      DO 80 IZ=0,Nlong

         Z=IZ
      DO 10 IY=0,Nlong
         Y=IY

         DO  IX=0,Nlong
            
        NAT=NAT+1
               X=IX

      If (((x-Sn-Nlong/2.)*DNX(1)+(y-Sn-Nlong/2.)*DNY(1)+
     *(z-Sn-Nlong/2.)*DNZ(1).le.0).and.((x+Sn-Nlong/2.)*DNX(2)
     * +(y-Sn-Nlong/2.)*DNY(2)
     * +(z-Sn-Nlong/2.)*DNZ(2).le.0).and.
     * ((x+Sn-Nlong/2.)*DNX(3)+(y+Sn-Nlong/2.)*DNY(3)+
     * (z-Sn-Nlong/2.)*DNZ(3).le.0).and.
     *((x-Sn-Nlong/2.)*DNX(4)+(y+Sn-Nlong/2.)*DNY(4)+
     * (z-Sn-Nlong/2.)*DNZ(4).le.0)
     * .and.
     *((x-Sn-Nlong/2.)*DNX(5)+(y-Sn-Nlong/2.)*DNY(5)+
     * (z+Sn-Nlong/2.)*DNZ(5).le.0)
     *.and.
     *((x+Sn-Nlong/2.)*DNX(6)+(y-Sn-Nlong/2.)*DNY(6)+
     * (z+Sn-Nlong/2.)*DNZ(6).le.0)
     * .and.
     *((x+Sn-Nlong/2.)*DNX(7)+(y+Sn-Nlong/2.)*DNY(7)+
     * (z+Sn-Nlong/2.)*DNZ(7).le.0)
     * .and.
     *((x-Sn-Nlong/2.)*DNX(8)+(y+Sn-Nlong/2.)*DNY(8)
     * +(z+Sn-Nlong/2.)*DNZ(8).le.0))
     *   IXYZ(NAT,1)=INT(X) 
               
        If (((x-Sn-Nlong/2.)*DNX(1)+(y-Sn-Nlong/2.)*DNY(1)+
     * (z-Sn-Nlong/2.)*DNZ(1).gt.0).or.((x+Sn-Nlong/2.)*DNX(2)
     *   +(y-Sn-Nlong/2.)*DNY(2)
     * +(z-Sn-Nlong/2.)*DNZ(2).gt.0).or.
     * ((x+Sn-Nlong/2.)*DNX(3)+(y+Sn-Nlong/2.)*DNY(3)
     *   +(z-Sn-Nlong/2.)*DNZ(3).gt.0).or.
     * ((x-Sn-Nlong/2.)*DNX(4)+(y+Sn-Nlong/2.)*DNY(4)
     *   +(z-Sn-Nlong/2.)*DNZ(4).gt.0)
     * .or.
     *((x-Sn-Nlong/2.)*DNX(5)+(y-Sn-Nlong/2.)*DNY(5)+
     *   (z+Sn-Nlong/2.)*DNZ(5).gt.0)
     * .or.
     *((x+Sn-Nlong/2.)*DNX(6)+(y-Sn-Nlong/2.)*DNY(6)+
     *   (z+Sn-Nlong/2.)*DNZ(6).gt.0)
     * .or.
     *((x+Sn-Nlong/2.)*DNX(7)+(y+Sn-Nlong/2.)*DNY(7)+
     *   (z+Sn-Nlong/2.)*DNZ(7).gt.0)
     * .or.
     * ((x-Sn-Nlong/2.)*DNX(8)+(y+Sn-Nlong/2.)*DNY(8)
     *   +(z+Sn-Nlong/2.)*DNZ(8).gt.0))
     *     goto 20
      
        If (IX.ne.IXYZ(nat,1)) write(*,*)'IX',IX,'IXYZ',IXYZ(nat,1)
          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

     
     
    
         
      

    
