!!! Program for fast athermal flipping of the 2D Penrose Tiling
!!! Usage: (echo MCS ; cat TILING.IN ) | ./ptmc > TILING.OUT
!!! MCS is number of sweeps per tiling node
!!! NOTE: array it(nodes,10) uses : basis e for first 5 entries, and inverted basis in exactly same order for the other 5 entries
!!! NOTE: tiling must include connectivity table, that can be created like this:
!!! cat tiling.in|./pttable > tiling_with_ctable

!!! Flipping is done exlusively by _computing_ indices of the connectivity table (without creating temporary neighbor lists) : this is the fastest way.

function iwrap(i,j) result(a)
  integer i,j,a
  a=i
  if(a>j)a=a-j
  if(a<1)a=a+j
end function iwrap
function ix(ii) result (a)
  implicit none
  integer, dimension(5) :: ii
  integer :: i,is,a
  a=0
  do i=1,5
     a=a+ii(i)
  enddo
end function
function ixa(ii) result (a)
  implicit none
  integer, dimension(5) :: ii
  integer :: i,is,a
  a=0
  do i=1,5
     a=a+abs(ii(i))
  enddo
end function
function ipo(v) result(a)
  integer v(5),a,s,i
  a=0
  s=0
  do i=1,5
     s=s+abs(v(i))
  enddo
  if(s==1)then
     do i=1,5
        if(v(i)==1)a=i
        if(v(i)==-1)a=i+5
     enddo
  endif
end function ipo
function bound5(v,ib,b,e) result(a)
  IMPLICIT double precision (A-H,O-Z)
  integer v(5),a(5),iw(2),ib(2,5)
  dimension b(2,2),r(2),s(2),e(2,5)
  r(:)=0
  do io=1,5
     r(:)=r(:)+v(io)*e(:,io)   ! cartesian coordinates
  enddo
  xb=b(1,1)*b(2,2)-b(1,2)*b(2,1)
  s(1)=(r(1)*b(2,2)-r(2)*b(2,1))/xb
  s(2)=(r(2)*b(1,1)-r(1)*b(1,2))/xb
  iw(:)=floor(s(:))
  a(:)=0
  do io=1,5
     a(io)=v(io)-iw(1)*ib(1,io)-iw(2)*ib(2,io)
  enddo
!  print "('bound5 ',5i4,'  x  ',5i4,'  sub ',2f7.3,2i3)",v,a,s,iw
end function bound5
recursive function ibound(ii) result(v)
 integer, dimension(5) :: ii,v
  integer :: is=0
  is=ix(ii)
  if(is>2)then
     ii(1:5)=ii(1:5)-1
     ii=ibound(ii)
  else if (is<-2) then
     ii(1:5)=ii(1:5)+1
     ii=ibound(ii)
  endif
  v=ii
end function ibound
function position(iv,e) result(a)
  IMPLICIT double precision (A-H,O-Z)
  integer, dimension(5) :: iv
  dimension a(2)
  dimension e(2,5)
  integer :: i,j,io
  a(1:2)=0
  do io=1,5
     a(1)=a(1)+iv(io)*e(1,io)
     a(2)=a(2)+iv(io)*e(2,io)
  enddo
end function position
function iboundapp(ii,ax,ay) result(v)
  integer, dimension(5) :: ii,v,ax,ay,iw,vv
  integer :: i,j,s,sa
  interface
     recursive function ibound(ii) result(v)
       integer, dimension(5) :: ii,v
       integer :: is=0
     end function ibound
  end interface
  sa=1000
  do i=-1,1
     do j=-1,1
        if(sa>1)then
           iw=ii+i*ax+j*ay
!           print *,i,j,"  xx",iw," yy",ii
           vv=ibound(iw)
           s=0
           do k=1,5
              s=s+abs(vv(k))
           enddo
           if(s<sa)then
              sa=s
              v=vv
           endif
        endif
     enddo
  enddo
end function iboundapp
        
program Main
  !CCCC plain athermal PT flipping
  IMPLICIT double precision (A-H,O-Z)
  integer, dimension(5) :: iy,iv,iu,iw    ! 5D integer vectors
  integer li(40)   ! list of nodes needed for flip update
  integer p1,p2,ll
  integer,parameter :: seed = 86456, tnodes=155000
  dimension e(2,5),iba(2,5),nit(tnodes)
  integer(4) it(tnodes,10),itt(100,10),itd(10)
  integer(4) ii,jj,mm,n,nn
  integer q1,q2,q3,q4,q5,q6,qq(3),pp(3),rr(3),ss(3),pt(3),pu(3),pq(3)
  integer, dimension(tnodes,5) :: pv
  dimension bar(2,2)
  interface
     function iwrap(i,j) result(a)
       integer i,j,a
     end function iwrap
     function position(iv,e) result(a)
       IMPLICIT double precision (A-H,O-Z)
       integer, dimension(5) :: iv
       dimension a(2)
       dimension e(2,5)
       integer :: i,j,io
     end function position
     function bound5(v,ib,b,e) result(a)
       IMPLICIT double precision (A-H,O-Z)
       integer v(5),a(5),iw(2),ib(2,5)
       dimension b(2,2),r(2),s(2),e(2,5)
     end function bound5
     function ix(ii) result (a)
       implicit none
       integer, dimension(5) :: ii
       integer :: i,is,a
     end function ix
     function ixa(ii) result (a)
       implicit none
       integer, dimension(5) :: ii
       integer :: i,is,a
     end function ixa
     recursive function ibound(ii) result(v)
       integer, dimension(5) :: ii,v
       integer :: is=0
     end function ibound
     function iboundapp(ii,ax,ay) result(v)
       integer, dimension(5) :: ii,v,ax,ay,iw,vv
       integer :: i,j,s,sa
     end function iboundapp
     function ipo(v) result(a)
       integer v(5),a,s,i
     end function ipo
  end interface
  read(*,*) mcs
  it(1:tnodes,1:10)=0
  pi=4.0*atan(1.0)
  do i=1,5
     e(1,i)=cos(2.0*pi*(i-1)/5.0)
     e(2,i)=sin(2.0*pi*(i-1)/5.0)
  enddo
  !      do i=1,5
  !         write(*,*)e(1,i),e(2,i)
  !     enddo
  read(*,*)iba(1,1:5)
  read(*,*)iba(2,1:5)
  bar(1,:)=position(iba(1,:),e)
  bar(2,:)=position(iba(2,:),e)
  read(*,*)nn
  do i=1,nn
     read(*,*)iy
     iv=ibound(iy)
     iv=bound5(iy,iba,bar,e)
     pv(i,:)=iv(:)
     p1=ix(iy)
!     print *,i,p1," xx  ",pv(i,:)
  enddo
  !  print *,"dadada"
  read(*,*)
  do i=1,nn
     read(*,*)it(i,:)
  enddo
  do i=1,nn
     nit(i)=0
     do j=1,10
        if(it(i,j)>0)then
           nit(i)=nit(i)+1
        endif
     enddo
!     print *,i,nit(i)
  enddo
!  do io=1,nn
!     print "(i5,'  ini-nit=',i3,'  it  ',10i3)",io,nit(io),it(io,:)
!  enddo
  mm=0   ! num. of flips
  do jj=1,mcs
     do ii=1,nn
        call random_number(x)
        n=x*nn+1
        if(nit(n)==3)then
           mm=mm+1
!           print *,jj,mm
!           print "('debug it b ',i5,' x  ',10i5)",n,it(n,:)
           ! we want to figure out directly the other three nodes needed for update.
           ! first we must decide whether angles are 4-2-4 (skihex) OR 3-4-3 (fathex)
           is1=0   ! neighbors from first 5 directions
           is2=0   ! neighbors from first 5 directions
           n3=0
           do io=1,10
              if(it(n,io)>0.and.io<6)is1=is1+1
              if(it(n,io)>0.and.io>5)is2=is2+1
              ! vertex position could be updated here, but is postponed later (for debugging)
!              if(io<6)pv(n,io)=pv(n,io)+1 ! add +1 along dir j
!              if(io>5)pv(n,io-5)=pv(n,io-5)-1 ! add -1 along dir j-5
           enddo
           ! is1 and is2 must be 3-0 or 2-1
           icase=0
           if(is1==3)icase=1   !skihex, first 5
           if(is2==3)icase=2   !skihex, second 5
           if(is1==2)icase=3   !fathex, angle "4" is first 5
           if(is2==2)icase=4   !fathex, angle "4" is second 5
           if(icase==0)then
              print *,"Connectivity update error: unknown case!"
              print *,"is1=",is1,"  is2=",is2
              print *,n," it ",it(n,:)
              goto 4444
           endif
           k0=0 
           if(icase==2.or.icase==4)then
              k0=5
           endif
           pt=(/2,1,3/) ! swap first and second index
           pu=(/1,2,1/) ! needed for triplet of new vertices
           ka=0
           if(icase<3)then  ! SKIHEX
              do k=1,5
!                 print *,"skihex ",k0+k,it(n,k0+k),iwrap(k0+k+1,5),it(n,iwrap(k0+k+1,5))
                 if(it(n,k0+k)>0.and.it(n,k0+iwrap(k+1,5))>0)then  ! xx0x0 0xx0x x0xx0 0x0xx x0x0x
                    ka=k
                    rr(1)=k
                    rr(2)=iwrap(k+1,5)
                    rr(3)=iwrap(k+3,5)
                    ss(1)=rr(3)
                    ss(2)=rr(3)
                    ss(3)=rr(2)
                 endif
              enddo
           endif
           if(icase>2)then  ! FATHEX
              do k=1,5
!                 print *,"fathex ",k0+k,it(n,k0+k),k0+iwrap(k+2,5),it(n,k0+iwrap(k+2,5))
                 if(it(n,k0+k)>0.and.it(n,k0+iwrap(k+2,5))>0)then !x0x00 0x0x0 00x0x x00x0 0x00x
                    ka=k
                    rr(1)=k
                    rr(2)=iwrap(k+2,5)
                    rr(3)=iwrap(k+1,5)+5
                    ss(1)=rr(3)
                    ss(2)=rr(3)
                    ss(3)=rr(2)
                 endif
              enddo
           endif
           do io=1,3
              pp(io)=it(n,iwrap(k0+rr(io),10))   ! three neighbors of the flippable vertex N
              qq(io)=it(pp(pu(io)),iwrap(k0+ss(io),10))  ! three HEX vertices that will become neighbors
           enddo
!test            li(1)=n
!test            li(2:4)=pp(:)
!test            li(5:7)=qq(:)
!test            nli=7
!test !!!!!! TEST before
!test            do i=2,7   ! add all first neighbors of the three nodes to the list
!test               ip=li(i)
!test               do j=1,10
!test                  iq=it(ip,j)
!test                  if(iq/=n.and.iq>0) then
!test                     ie=0
!test                     do k=1,m ! list through existing nodes and add to the list only when needed
!test                        if(ie==0) then
!test                           if(iq==li(k))ie=1
!test                        endif
!test                     enddo
!test                     if(ie==0)then
!test                        nli=nli+1
!test                        li(nli)=iq
!test                     endif
!test                  endif
!test               enddo
!test            enddo  ! finished constructing nodelist for update
!test            print *,"Testing list : nli = ",nli
!test            itt(1:nli,:)=0   ! this will be for current node
!test            m=0
!test            ier=0
!test            do i=1,nli-1
!test               do j=i+1,nli
!test                  iu(:)=pv(li(i),:)
!test                  iv(:)=pv(li(j),:)
!test                  iw=iboundapp(iv-iu,iba(1,:),iba(2,:))
!test                  k=ipo(iw)  ! k>0 only when difference between two nodes is ONE
!test                  if(k>0)then
!test                     itt(i,k)=li(j)
!test                     if(itt(i,k)/=it(li(i),k))ier=1
!test                     if(k<6)m=k+5
!test                     if(k>5)m=k-5
!test                     itt(j,m)=li(i)
!test                     if(itt(j,m)/=it(li(j),m))ier=1
!test                  endif
!test               enddo
!test            enddo
!test !           if(ier==1)then
!test               print *,"ERROR test-before - connectivity tables n, ier : ",n,ier
!test               print *,n," IT",it(n,:)
!test               print *,n,"ITT",itt(1,:)
!test               do j=2,7
!test                  print *,li(j)," it",it(li(j),:)
!test                  print *,li(j),"itt",itt(j,:)
!test               enddo
!test !              goto 4444
!test !           endif
!test !!!!!!! end TEST before           
           ier=0
           do io=1,3
              if(pp(io)==0)ier=1
              if(qq(io)==0)ier=1
           enddo
           if(ier==1)then
              do io=1,3
                 print "('debug it b ',i5,' pp ',10i5)",pp(io),it(pp(io),:)
              enddo
              do io=1,3
!                 print "('debug it b ',i5,' qq ',10i5)",qq(io),it(qq(io),:)
              enddo
              print *,"------------------------------------------------------"
              print *,"error node-list pp qq : ",pp(:),qq(:)
              print *,"error it : ",it(n,:)
              print *,"error k0 ka :  ",icase,k0,ka
              print *,"error rr ss : ",rr(:),ss(:)
!              goto 4444
           endif
           do io=1,10! UPDATE POSITION HERE!
              if(it(n,io)>0) then
                 if(io<6)pv(n,io)=pv(n,io)+1 ! add +1 along dir j
                 if(io>5)pv(n,io-5)=pv(n,io-5)-1 ! add -1 along dir j-5
              endif
           enddo
           do io=1,3   ! UPDATE TABLES
              ia=0
              if(io<3)ia=5
              ip=k0+rr(io)
              it(n,iwrap(ip,10))=0   ! delete old neighbors of node N
              it(n,iwrap(ip+5,10))=qq(pt(io))  ! pt - because we need to swap index 1 and 2
              it(pp(io),iwrap(ip+5,10))=0   ! delete node N from entry of q1..q3
              it(qq(pt(io)),iwrap(ip,10))=n
          enddo
          nit(pp(:))=nit(pp(:))-1
          nit(qq(:))=nit(qq(:))+1
          q1=iwrap(k0+rr(1),10)
          q2=iwrap(k0+rr(2),10)
          q3=iwrap(k0+rr(3),10)
!debug          print *,"XXXXXXXXX ",n,rr(:),it(qq(1),q2),it(qq(2),q1),it(qq(3),q3)
!debug          print *,"debug node-list pp qq : ",pp(:),qq(:)
!debug          print *,"debug k0 ka :  ",icase,k0,ka
!debug          print *,"debug rr ss : ",rr(:),ss(:)
!debug          print "('debug it a ',i5,' x  ',10i5)",n,it(n,:)
!debug          do io=1,3
!debug             print "('debug it a ',i5,' pp ',10i5)",pp(io),it(pp(io),:)
!debug          enddo
!debug          do io=1,3
!debug             print "('debug it a ',i5,' qq ',10i5)",qq(io),it(qq(io),:)
!debug          enddo
!debug          print *,"------------------------------------------------------"
!TEST   after
!test           itt(1:nli,:)=0   ! this will be for current node
!test           m=0
!test           ier=0
!test           do i=1,nli-1
!test              do j=i+1,nli
!test                 iu(:)=pv(li(i),:)
!test                 iv(:)=pv(li(j),:)
!test                 iw=iboundapp(iv-iu,iba(1,:),iba(2,:))
!test                 k=ipo(iw)  ! k>0 only when difference between two nodes is ONE
!test                 if(k>0)then
!test                    itt(i,k)=li(j)
!test                    if(itt(i,k)/=it(li(i),k))ier=1
!test                    if(k<6)m=k+5
!test                    if(k>5)m=k-5
!test                    itt(j,m)=li(i)
!test                    if(itt(j,m)/=it(li(j),m))ier=1
!test                 endif
!test              enddo
!test           enddo
!test           if(ier==1)then
!test              print *,"ERROR test-after - connectivity tables : n, ier ",n,ier
!test              print *,n," IT",it(n,:)
!test              print *,n,"ITT",itt(1,:)
!test              do j=2,7
!test                 print *,li(j)," it",it(li(j),:)
!test                 print *,li(j),"itt",itt(j,:)
!test              enddo
!test              goto 4444
!test           endif
       endif  ! end nit(n)==3
    enddo  ! end ii loop over nodes
     if(mod(jj,100)==0) then  ! test node for wrap-around every 100 MCS
        do ii=1,nn
           pv(ii,:)=bound5(pv(ii,:),iba,bar,e)
        enddo
     endif
  enddo  ! end jj MCS loop nn
  x1=mm
  x2=mcs*nn
3333 continue
  print *, iba(1,:)
  print *, iba(2,:)
  print *, nn," 1"
  do i=1,nn
     print *, pv(i,:),"     0"
  enddo
!  print *,"flipstat",mm,mcs*nn,x1/x2
  print *,"ctable"
  do i=1,nn
     print *, it(i,:)
  enddo
4444 continue
end program Main
      
