C****************************************************************** C This program calculate the vertex positions of 2D quasiperiodic tilings C APPROXIMANTS using the generalized dual method. C 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 C APPROXIMANTS are computed by using rational approximants for the star vectors C For a Penrose, this is done by C using an appropiate basis star. Here we present two options C as detailed below C If you have comments, contact: G. Naumis naumis@fisica.unam.mx C********************************************************************* C Variable declaration DOUBLE PRECISION EX(5),EY(5),UX(5),UY(5),GAMMA(5),rx,ry,SNORME(5) DOUBLE PRECISION arx0,ary0,aux,AL,PI,SGAMMA,NORMA,arxd,aryd DOUBLE PRECISION Aij,arx,ary, TAU,SIGMA,BETA,ai,aj DOUBLE PRECISION VX(5),VY(5),DOT1,DOT2,DOT3,DOT4,VXD(5),VYD(5) DOUBLE PRECISION PHD(5),QHD(5),OMEGA(5,5),MU(5,5),ALFA(5) DOUBLE PRECISION PHDP(5),QHDP(5),RM,RN,EXN(5),EYN(5) DOUBLE PRECISION MUNORM(5),NORMOME(5) DOUBLE PRECISION PHDNORM,QHDNORM,PQNORM DOUBLE PRECISION TAUPLUS,BETAPLUS,SIGMAPLUS COMMON H,NTP,U,V,PI,N0,NTG,N,NFR integer iseed,irand,i,ID(200),KD(200),KA,P,Q,FIBONUM(100),FMAX integer AMUX,NK(5),SUMK,SI,NDF,CELL integer NP,NQ,PPLUS,QPLUS 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.231040 GAMMA(2)=-0.33976 GAMMA(3)=-0.144633 GAMMA(4)=-0.110439 GAMMA(5)=-1.0*(GAMMA(1)+GAMMA(2)+GAMMA(3)+GAMMA(4)) SGAMMA=GAMMA(1)+GAMMA(2)+GAMMA(3)+GAMMA(4)+GAMMA(5) write(*,*) 'Gammas',GAMMA(1),GAMMA(2),GAMMA(3),GAMMA(4),GAMMA(5) C========================================================================== C TYPE OF UNIT CELL CHOOSE TWO VALUES C CELL=1 RECTANGLE, CELL=0 THIN RHOMBUS C C========================================================================= CELL=1 C************************************************************************ C Define approximant of TAU for Penrose tiling approximant C Fibonacci Numbers using FMAX C TAU=P/Q where P is the FMAX th Fibonacci number C C For CELL=1, works for FMAX>=4 C For CELL=0 works for FMAX>=3 C************************************************************************ FMAX=4 FIBONUM(1)=1 FIBONUM(2)=1 DO 2 i=3,FMAX FIBONUM(i)=FIBONUM(i-1)+FIBONUM(i-2) 2 CONTINUE 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************************************************************************ PPLUS=FIBONUM(FMAX) QPLUS=FIBONUM(FMAX-1) P=FIBONUM(FMAX-1) Q=FIBONUM(FMAX-2) TAU=DBLE(P)/DBLE(Q) SIGMA=TAU-1.0 BETA=DBLE(3.0-TAU) TAUPLUS=DBLE(P)/DBLE(Q) SIGMAPLUS=TAU-1.0 BETAPLUS=DBLE(3.0-TAU) C write(*,*) 'Computing using TAU,P,Q',TAU,P,Q write(*,*)'Using APPROXIMANT P and Q:',PPLUS,QPLUS write(*,*) 'Isomorphism class given by summing GAMMA phases=',SGAMMA C======================================================================== C Define P and Q vectors in plane P_{pq} C We modified the method of getting rational approximants by C Entin-Wohlman,Klemant, J. Phys. France, Vol. 49, 567 (1988) C THIS PRODUCES FATS AND THIN RHOMBUS UNIT CELL C======================================================================== if (CELL.EQ.0) then PHD(1)=0.0 PHD(2)=DBLE(PPLUS) PHD(3)=DBLE(P) PHD(4)=-1.0*DBLE(P) PHD(5)=-1.0*DBLE(PPLUS) QHD(1)=-1.0*DBLE(QPLUS) QHD(2)=0.0 QHD(3)=1.0*DBLE(QPLUS) QHD(4)=1.0*DBLE(Q) QHD(5)=-1.0*DBLE(Q) C write(*,*) 'P VECTOR',PHD(1),',',PHD(2),PHD(3),PHD(4),PHD(5) C write(*,*) 'Q VECTOR',QHD(1),',',QHD(2),QHD(3),QHD(4),QHD(5) C====================================================================== C Projection of P and Q vectors over irrational plane P C====================================================================== DO J=1,5 PHDP(J)=0.0 QHDP(J)=0.0 DO I=1,5 ai=DBLE(I) aj=DBLE(J) PHDP(J)=PHDP(J)+PHD(I)*cos((ai-aj)*AL) QHDP(J)=QHDP(J)+QHD(I)*cos((ai-aj)*AL) ENDDO ENDDO C======================================================================= write(*,*) 'P VECTOR',PHD(1),',',PHD(2),PHD(3),PHD(4),PHD(5) write(*,*) 'P in PLANE P',PHDP(1),',',PHDP(2),PHDP(3),PHDP(4),PHDP(5) DOT3=0.0 PHDNORM=0.0 QHDNORM=0.0 C====================================================================== C Calculate GRID lines as vector omega C====================================================================== DO 4 I=1,5 DO 5 J=1,5 OMEGA(I,J)=QHD(I)*PHD(J)-PHD(I)*QHD(J) 5 CONTINUE DOT3=DOT3+PHDP(I)*QHDP(I) PHDNORM=PHDNORM+PHDP(I)*PHDP(I) QHDNORM=QHDNORM+QHDP(I)*QHDP(I) 4 CONTINUE C write(*,*) 'OMEGA(1,1 a 3',OMEGA(1,1),OMEGA(1,2),OMEGA(1,3) C write(*,*) 'OMEGA(1,1 a 4)',OMEGA(1,4),OMEGA(1,5) C write(*,*) C write(*,*) 'OMEGA(2,1 a 3',OMEGA(2,1),OMEGA(2,2),OMEGA(2,3) C write(*,*) 'OMEGA(2,1 a 4)',OMEGA(2,4),OMEGA(2,5) C======================================================================= C Calculate angles between grids C======================================================================= C Norm of OMEGA DO I=1,5 DOT1=0.0 DO j=1,5 DOT1=DOT1+OMEGA(I,J)*OMEGA(I,J) ENDDO NORMOME(I)=DSQRT(DOT1) ENDDO C Dot product and angle between omegas ALFA(1)=0.0 DO i=2,5 DOT3=0.0 DO j=1,5 DOT3=DOT3+OMEGA(i,j)*OMEGA(1,j) ENDDO DOT4=DOT3/(NORMOME(i)*NORMOME(1)) C DOT4=DOT3 ALFA(i)=ACOS(DOT4) write(*,*) 'Vector 1 with i, DOT3,DOT4',i,DOT3,NORMOME(i) write(*,*) 'DOT4',DOT4 write(*,*) 'Alfa ANGLE in DEGREES',ALFA(i)*360.0/(2.0*3.1416) ENDDO DO 13 i=1,3 EX(i)=1.0*NORMOME(i)*cos(ALFA(i))/NORMOME(1) EY(i)=1.0*NORMOME(i)*sin(ALFA(i))/NORMOME(1) write(*,*) 'APPROXIMANT STAR:',i,EX(i),EY(i) 13 CONTINUE DO i=4,5 EX(i)=1.0*NORMOME(i)*cos(ALFA(i))/NORMOME(1) EY(i)=-1.0*NORMOME(i)*sin(ALFA(i))/NORMOME(1) ENDDO write(*,*) 'APPROXIMANT STAR:',i,EX(i),EY(i) write(*,*) 'SUM OF E IN X',EX(1)+EX(2)+EX(3)+EX(4)+EX(5) write(*,*) 'SUM OF E IN Y',EY(1)+EY(2)+EY(3)+EY(4)+EY(5) endif C HERE ENDS FIRST OPTION of CELL C HERE STARTS SECOND OPTION RECTANGULAR CELL if (CELL.EQ.1) then C========================================================================= C Alternative rational star obtained by considering a rationalization C of a rotation matrix in 5D, following Lord, Ramaksrihna C PRODUCES ORTHOROMBIC (RECTANGLE) UNIT CELL if succesive rationals C of tau are taken C Obtained by writing the angles in a pentagon in terms of tau C and using approximants for it. It can be used two different C approximants for EX and EY. If they are succesive of C tau, a rectangular unit cell is found C========================================================================= EX(1)=1.0 EY(1)=0.0 EX(2)=SIGMAPLUS/2.0 EY(2)=(TAU*BETA)/2.0 EX(3)=(-1.0*TAUPLUS)/2.0 EY(3)=BETA/2.0 EX(4)=(-1.0*TAUPLUS)/2.0 EY(4)=(-1.0*BETA)/2.0 EX(5)=SIGMAPLUS/2.0 EY(5)=(-1.0*TAU*BETA)/2.0 C DOT4=REAL(P*P+Q*Q) C EX(1)=REAL(P*P)/DOT4 C EY(1)=-REAL(Q)/SQRT(DOT4) C EX(2)=REAL(P*P)/DOT4 C EY(2)=REAL(Q)/SQRT(DOT4) C EX(3)=-REAL(Q*Q)/(DOT4) C EY(3)=REAL(P)/SQRT(DOT4) C EX(4)=-2.0*REAL(P*Q)/DOT4 C EY(4)=0.0 C EX(5)=-REAL(Q*Q)/(DOT4) C EY(5)=REAL(P)/SQRT(DOT4) endif C====================================================================== C END SECTION OPTION OF UNIT CELL C====================================================================== C====================================================================== C STARS NOT NORMALIZED C For approximant, the basis is not normalized. In principle, we can trick C the program by using a new basis which replaces the base EN C This basis is normalized to the SQUARE. However, one can show C that the result is the same as EX and EY. So I keept this base C in case ones want to check this c====================================================================== DO i=1,5 SNORME(i)=(EX(i)*EX(i)+EY(i)*EY(i)) EXN(i)=EX(i)/SNORME(i) EYN(i)=EY(i)/SNORME(i) C EXN(i)=EX(i) C EYN(i)=EY(i) write(*,*) 'Final Basis Used E',i,EX(i),EY(i) ENDDO C*********************************************************************** C When approximants are used, the dual has two options. One is using C the same vectors as the grid, given by vectors (EX(i),EY(i)) C or have the same symmetry of the limiting case using a set C of vectors pointing to the star in dual space (VX(i),VY(j)) C By default, the program uses the base (VX,VY) C*********************************************************************** DO 15 i=1,5 C========================================================================= C Base in PQ-PLANE where the approximant lives C It produces 10 different tiles for Penrose C========================================================================= ai=DBLE(i)-1.0 VX(i)=EX(i) VY(i)=EY(i) C========================================================================= C Base in the irrational plane where a pentagonal star is used C It produces 2 kinds of tilings for Penrose C Usually, this is the required tiling C========================================================================== VXD(i)=cos(ai*AL) VYD(i)=sin(ai*AL) 15 CONTINUE C************************************************************************* C Compute orthonal vectors to star C************************************************************************* DO 19 i=1,5 UX(i)=-EY(i) UY(i)=EX(i) 19 CONTINUE write(*,*)'Loop' 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 "Penrose.xyz" is a file in XYZ format, to be read in Xcrysden C************************************************************************** OPEN(2,FILE='penroseav.dat',STATUS='UNKNOWN') OPEN(4,FILE='penrosed.dat',STATUS='UNKNOWN') OPEN(6,FILE='penroseint.dat',STATUS='UNKNOWN') OPEN(10,FILE='penrose.xyz',STATUS='UNKNOWN') OPEN(12,FILE='indices.dat',STATUS='UNKNOWN') OPEN(14,FILE='penrose.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,4 DO 20 j=i+1,5 C write(*,*) 'i,j',i,j Aij=((EX(i)*EY(j))-(EY(i)*EX(j))) NP=-1*PPLUS-5 NQ=-1*QPLUS-5 C write(*,*) 'n,m',NP,NQ DO 30 n=NP,(1*PPLUS)+5 DO 30 m=NQ,(1*QPLUS)+5 C============================================================================== C Generate intersection of grids using perpendicular vectors C arx0 and ary0 are coordinates of intersection C For approximants, the star is not normalized corrections for taken C into account this are needed C============================================================================== RM=(DBLE(m)+GAMMA(j)) RN=(DBLE(n)+GAMMA(i)) arx0=(-1.0*RM*EY(i)+RN*EY(J))/Aij ary0=(RM*EX(i)-RN*EX(j))/Aij co1=(arx0*EX(i))+(ary0*EY(i))-DBLE(n)-GAMMA(i) co2=(arx0*EX(j))+(ary0*EY(j))-DBLE(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 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**************************************************************************** C============================================================================== C Compute dual with n and m C This can be done using two basis C In the aproximant basis (VX,VY) or in the irrational C plane using star irrational basis (VXD,VYD) which is a pentagonal star C For approximants, E is not normalized. Corrections are needed into C obtainig integer coordinates. Thus we use a basis E/(1E1*1E1) C============================================================================== rx=(DBLE(n)*VX(i))+((DBLE(m))*VX(j)) ry=(DBLE(n)*VY(i))+((DBLE(m))*VY(j)) rxd=(DBLE(n)*VXD(i))+(DBLE(m)*VXD(j)) ryd=(DBLE(n)*VYD(i))+(DBLE(m)*VYD(j)) arxp=rx aryp=ry NK(i)=n NK(j)=m DO 40 L=1,5 if ((L.NE.i).AND.(L.NE.j)) then aux=((arx0*EX(L))+(ary0*EY(L))) aux=aux-GAMMA(L) if (aux.GE.0.0) then n1=AINT(aux) else n1=AINT(aux)-1 endif NK(L)=n1 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+(DBLE(n1)*VX(L)) ary=ry+(DBLE(n1)*VY(L)) arxp=arxp+((aux-0.5)*VX(L)) aryp=aryp+((aux-0.5)*VY(L)) rx=arx ry=ary arxd=rxd+(DBLE(n1)*VXD(L)) aryd=ryd+(DBLE(n1)*VYD(L)) rxd=arxd ryd=aryd endif 40 CONTINUE C n1=-(NK(1)+NK(2)+NK(3)+NK(4)) C NK(5)=n1 C arx=rx+(REAL(n1)*VX(5)) C ary=ry+(REAL(n1)*VY(5)) C arxp=arxp+((aux-0.5)*VX(5)) C aryp=aryp+((aux-0.5)*VY(5)) C rx=arx C ry=ary SUMK=(NK(1)+NK(2)+NK(3)+NK(4)+NK(5)) 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 write(*,*) 'SUMK',SUMK if ((ABS(rx).LT.12.0).AND.(ABS(ry).LT.12.0)) then if ((rx.GE.-5.2).AND.(ry.GE.-5.2)) then write(4,*) rx,ry,SUMK write(4,*) rx-VX(i),ry-VY(i),SUMK write(4,*) rx-VX(j),ry-VY(j),SUMK write(4,*) rx-VX(i)-VX(j),ry-VY(i)-VY(j),SUMK write(14,*) rxd,ryd,SUMK write(14,*) rxd-VXD(i),ryd-VYD(i),SUMK write(14,*) rxd-VXD(j),ryd-VYD(j),SUMK write(14,*) rxd-VXD(i)-VXD(j),ryd-VYD(i)-VYD(j),SUMK C if (ABS(SUMK).EQ.0) then C write(14,*) rx,ry,SUMK C endif C write(14,*) rx-VXD(i),ry-VYD(i),SUMK-1 C write(14,*) rx-VXD(j),ry-VYD(j),SUMK-1 C write(14,*) rx-VXD(i)-VXD(j),ry-VYD(i)-VYD(j),SUMK-2 C endif 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 write(10,*) 'C',rx,ry,1.0 write(10,*) 'C',rx-EX(i),ry-EY(i),1.0 write(10,*) 'C',rx-EX(j),ry-EY(j),1.0 write(10,*) 'C',rx-EX(i)-EX(j),ry-EY(i)-EY(j),1.0 write(12,*) i,j,rxd,ryd,NK(1),NK(2),NK(3),NK(4),NK(5),SUMK endif endif 30 CONTINUE 20 CONTINUE C CLOSE(2) C CLOSE(4) C CLOSE(6) C CLOSE(10) C CLOSE(12) C CLOSE(14) C************************************************************** C This feature makes a nice drawing of the family of lines C************************************************************** OPEN(8,FILE='lineas.dat',STATUS='UNKNOWN') AMUX=8 DO 120 i=1,5 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)-(REAL(AMUX)/2.0) y=(REAL(n)-(REAL(AMUX)/2.0)+GAMMA(i)-x*EX(i))/EY(i) else y=(REAL(naux)*(0.01)-(REAL(AMUX)/2.0)+GAMMA(i)) x=REAL(n)-(REAL(AMUX)/2.0)+GAMMA(i) endif write(8,*) x,y,1.0 119 CONTINUE 120 CONTINUE CLOSE(8) STOP END