C******************************************************************
C  This program calculate the vertex positions of 2D quasiperiodic tilings
C  using the generalized dual method. It is just a simple program
C  from were you can start other applications like the 3D version (it is just
C  very easy to do,from the present code).
C  First, you give the star vectors e_{i},
C  i.e. a set of unitary vectors that have the required symmetry. For example,
C  in the Penrose it is required 5 fold symmetry, thus, the star is a set
C  of 5 vectors that point to the vertex of a pentagon. Then, all intersections
C  of a family of lines perpendicular to each vector is calculated. Following
C  the formulas of G.G. Naumis, J.L. Aragón,"Analytic expressions for the 
C  vertex coordinates of quasiperiodic lattices",Z.  Kristallogr. 218, 397 (2003).
C  This program has been made by Gerardo G. Naumis and José Luis Aragon,
C  in 2002, at the Instituto de Fisica, National University of Mexico.
C  To draw the tiling, just make a program that joins all points separated
C  by a distance "1". See our mathematica version of the program in the
C  same webpage:  www.fisica.unam.mx/naumis
C  If you have comments, contact: G. Naumis naumis@fisica.unam.mx
C*********************************************************************
C  Variable declaration
      REAL EX(5),EY(5),UX(5),UY(5),GAMMA(5),rx,ry,arx0,ary0,aux,AL,PI
      REAL Aij,arx,ary
      COMMON H,NTP,U,V,PI,N0,NTG,N,NFR
      integer iseed,irand,i,ID(200),KD(200),KA
      double precision rand
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C    Parameters
       
       PI=4.0*ATAN(1.0)
       AL=2.0*PI/5.0
C*******************************************************************
C  Following the notation of Steinhardt and Levine, on needs to define
C  phases that shifts the family of lines on each vector. 
C  Different isomorphisms (different types of tilings) are obtained 
C  if these phases are changed. If the sum of these phases is 0, 
C  one obtain the Penrose isomorphism class, which is the most famous.
C  see: J Socolar, Steinhardt, Phys. Rev. B 34, 617 (1986).
C  However, many people forget to tell you that you must be careful 
C  with singular tilings,i.e., if these shifts are so that
C  one has points where three of more lines cross, the tiles are no longer
C  rhombuses, one can have hexagons (duals of three ´lines that meet in a point)
C  octagons (four lines), decagon (5 lines). This last one ocuurs at the origin
C  if all phases are set to zero.
C  A FAQ is "why I don´t get the right Penrose?. First, you will always
C  get different portions of the Penrose tiling unless you know which are the
C  exact parameters of the tile that you want to compare. Usually,instead one
C  needs to check that both tilings are in the same local isomorphism class,
C  that means that each patch of a given radius, is observed in both tilings.
C  To check the Penrose, as first check, verify the proportions of each 
C  vertex configuration and second neighbour. A table was published in:
C   G.G. Naumis, "Density of states and first spectral moments of a 
C  quasiperiodic lattice", J. Phys: Condens. Matter 11 (1999) 7 143-7 154.
C  I an exotic symmetry is required, change the number of star vectors
C************************************************************************  
       GAMMA(1)=0.0
       GAMMA(2)=-0.50
       GAMMA(3)=0.25
       GAMMA(4)=0.25
       GAMMA(5)=0.44
C       write(*,*) AINT(0.1135)
C       write(*,*) AINT(2.5760)
C       write(*,*) AINT(-1.2563)
C************************************************************************
C      Definition of star (basis vectors), and generation of auxiliary
C      vectors that allow to get the intersection of each family of lines
C************************************************************************
       DO 10 i=1,5
       ri=REAL(i)-1.0
       EX(i)=cos(AL*ri)
       EY(i)=sin(AL*ri)
       UX(i)=-sin(AL*ri)
       UY(i)=cos(AL*ri)
       write(*,*) EX(i),EY(i)
10     CONTINUE
C**************************************************************************
C  Just for the record, two files are generated as output, one with the
C  intersection points of each family ("pensorseint.dat"). This is known 
C  as the skeleton. The other file, "penrsoe.dat" contains the vertex positions.
C  The Average lattice is also calculated and store in "penroseav.dat"
C**************************************************************************
      OPEN(2,FILE='penroseav.dat',STATUS='UNKNOWN')
      OPEN(4,FILE='penrose.dat',STATUS='UNKNOWN')
      OPEN(6,FILE='penroseint.dat',STATUS='UNKNOWN')
C**************************************************************************
C A loop is carried over I and J which are couples of vectors that generate
C all the intersections. For each copule I and J, we have a family of lines
C indexed by m and n, which are the "numeration" of each line over the directions
C I and J.
C**************************************************************************
      DO 20 i=1,5
       DO 20 j=i+1,5
        if (i.NE.j) then 
         Aij=((EX(i)*EY(j))-(EY(i)*EX(j)))

         DO 30 n=-10,10
         DO 30 m=-10,10
          rx=(REAL(n)*EX(i))+(REAL(m)*EX(j))
          ry=(REAL(n)*EY(i))+(REAL(m)*EY(j))
         arx0=(-(REAL(m)+GAMMA(j))*EY(i)+(REAL(n)+GAMMA(i))*EY(j))/Aij
         ary0=(((REAL(m)+GAMMA(j))*EX(i)-(REAL(n)+GAMMA(i))*EX(j)))/Aij
         co1=(arx0*EX(i))+(ary0*EY(i))-REAL(n)-GAMMA(i)
         co2=(arx0*EX(j))+(ary0*EY(j))-REAL(m)-GAMMA(j)
C***************************************************************************
C Write the family and interseection coordinates of two lines,
C which in fact correspond to the center of a rhombus.
C***************************************************************************
C         write(*,*) 'N',N,' M',M,' INT=',arx0,ary0
          if ((ABS(co1).GT.0.0001).OR.(ABS(co2).GT.0.0001)) then
           write(*,*) 'error',I,J,n,m,co1,co2
          endif
C         write(*,*) 'I',I,' J',J,' N',N,' M',M,co1
C****************************************************************************
C Dual transformation 
C Calculate ordinal coordinates with respect to each family of lines
C of the generated intersection. Of course, coordinates I and J are known
C already, since they were used to generate the point. Thus, such step
C is not performed for lines of family I and J.
C****************************************************************************
          arxp=rx
          aryp=ry
         DO 40 L=1,5
          if ((L.NE.i).AND.(L.NE.j)) then
          aux=(arx0*EX(L))+(ary0*EY(L))-GAMMA(L)
           if (aux.GE.0.0) then
           n1=AINT(aux)
           else
            n1=AINT(aux)-1
           endif
C******************************************************************
C Perform in a recursive way the dual transformation by making a linear
C combination with the right coefficients of the star vectors.
C******************************************************************
           arx=rx+(REAL(n1)*EX(L))
           ary=ry+(REAL(n1)*EY(L))
           arxp=arxp+((aux-0.5)*EX(L))
           aryp=aryp+((aux-0.5)*EY(L))
           rx=arx
           ry=ary
          endif
40    CONTINUE
C  Once a point is found in the tiling, the corresponding other 3 points of
C  a rhombus are generated by adding the proper vectors
        if ((ABS(rx).LT.9.0).AND.(ABS(ry).LT.9.0)) then
        write(6,*) arx0,ary0
        write(4,*) rx,ry,1.0
        write(4,*) rx-EX(i),ry-EY(i),1.0
        write(4,*) rx-EX(j),ry-EY(j),1.0
        write(4,*) rx-EX(i)-EX(j),ry-EY(i)-EY(j),1.0
        write(2,*) arxp,aryp,1.0 
        write(2,*) arxp-EX(i),aryp-EY(i),1.0  
        write(2,*) arxp-EX(j),aryp-EY(j),1.0  
        write(2,*) arxp-EX(i)-EX(j),aryp-EY(i)-EY(j),1.0     
        endif
30    CONTINUE
        endif
20    CONTINUE
      CLOSE(2) 
      CLOSE(4) 
      CLOSE(6)  
C**************************************************************
C This feature makes a nice drawing of the family of lines
C**************************************************************	   
      OPEN(8,FILE='lineas.dat',STATUS='UNKNOWN')
       AMUX=8.0
       DO 120 i=1,5
C         write(8,*) EX(i),EY(i)
        Do 119 n=1,AMUX
         DO 119 naux=1,1000
         if (i.NE.1) then
         x=REAL(naux)*(0.01)-(AMUX/2.0)
         y=(REAL(n)-(AMUX/2.0)+GAMMA(i)-x*EX(i))/EY(i)
         else
          y=(REAL(naux)*(0.01)-(AMUX/2.0)+GAMMA(i))
          x=REAL(n)-(AMUX/2.0)+GAMMA(i)
         endif
        write(8,*) x,y,1.0
119    CONTINUE
120    CONTINUE
       CLOSE(8)
      STOP   
      END

