!!! Program flipping of the 2D Penrose Tiling
!!! Hamiltonian is "fat-fat" side by side
!!! 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
function v5v2p(iv,ep,bap) result(v)
  IMPLICIT double precision (A-H,O-Z)
  dimension v(2),bap(2,5),x(2)
  integer io,iv(5)
  x(:)=0.0
  do io=1,5
     x(1)=iv(io)*ep(1,io)
     x(2)=iv(io)*ep(2,io)
  enddo
  v(1)=sqrt(5.)*2*(x(1)*bap(2,2)-x(2)*bap(2,1))/(bap(1,1)*bap(2,2)-bap(1,2)*bap(2,1))
  v(2)=sqrt(5.)*2*(x(2)*bap(1,1)-x(1)*bap(1,2))/(bap(1,1)*bap(2,2)-bap(1,2)*bap(2,1))
end function v5v2p
subroutine perpsum(iv,is,ep,ps,ps2) ! add/subtract mean perp and mean square perp
  IMPLICIT double precision (A-H,O-Z)
  integer is,js,io,iv(5)
  dimension r(2),ps(3),ps2(3),a(3),ep(2,5)
  r(:)=0.
  do io=1,5
     r(:)=r(:)+iv(io)*ep(:,io)
  enddo
  a(1:2)=r(:)
  js=0
  do io=1,5
     js=js+iv(io)
  enddo
  js=mod(js,5)
  a(3)=js
  ps(:)=ps(:)+is*a(:)
  ps2(:)=ps2(:)+is*a(:)**2
  return
end subroutine perpsum
function wperp(ps,ps2,n) result(wp)
  IMPLICIT double precision (A-H,O-Z)
  integer(8) n
  dimension wp(2),ps(3),ps2(3)
  wp(1)=ps2(1)/n+ps2(2)/n-(ps(1)/n)**2-(ps(2)/n)**2
  wp(2)=ps2(3)/n-(ps(3)/n)**2    ! discrete perp space
end function wperp
subroutine rround(x,i,dx)
  IMPLICIT double precision (A-H,O-Z)
  integer i
  i=x
  dx=x-float(i)
  if(dx.gt.0.5)then 
     i=i+1
     dx=x-i
  endif
  if(dx.lt.-0.5)then
     i=i-1
     dx=x-i
  endif
  return
end subroutine rround
recursive function pathwalk(node,dir,turns,po,mxpo,td,it) result(lastnode)
  integer,parameter :: tnodes=155000
  integer ia(2),td(-5:5,10),it(tnodes,10),k,j,turns(10),lastnode,node,po,dir,mxpo,nextnode
  k=td(turns(po),dir)  ! new direction
  nextnode=it(node,k)     ! new node
!  print *,"in-pathwalk node dir : ",lastnode,k
  if(nextnode>0.and.po+1<mxpo+1)then
     lastnode=pathwalk(nextnode,k,turns,po+1,mxpo,td,it)
  else
     lastnode=nextnode
  endif
end function pathwalk

!!!!!  MAIN

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, pdim=600
  dimension e(2,5),ep(2,5),iba(2,5),nit(tnodes)
  integer(4) it(tnodes,10),itt(100,10),itd(10),pth(20)
  dimension d(-pdim:pdim,-pdim:pdim),cn(4,2)
  integer(8) ii,jj,mm,nn,mcs,n,nnff,nnss
  integer q1,q2,q3,q4,q5,q6,qq(3),pp(3),rr(3),ss(3),pt(3),pu(3),rp(3)
  integer tr1(10),tl1(10),tr2(10),tl2(10),tr3(10),tl3(10),tr4(10),tl4(10),tr5(10),tl5(10)
  integer turns(10)
  integer td(-5:5,10)
  integer, dimension(tnodes,5) :: pv
  dimension bar(2,2),bap(2,2),abap(2),ps(3),ps2(3),wp(2),xh(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
     function v5v2p(iv,ep,bap) result(v)
       IMPLICIT double precision (A-H,O-Z)
       dimension v(2),bap(2,5),x(2)
       integer io,iv(5)
     end function v5v2p
     subroutine perpsum(iv,is,ep,ps,ps2) ! add/subtract mean perp and mean square perp
       IMPLICIT double precision (A-H,O-Z)
       integer is,js,io,iv(5)
       dimension r(2),ps(3),ps2(3),a(3),ep(2,5)
     end subroutine perpsum
     function wperp(ps,ps2,n) result(wp)
       IMPLICIT double precision (A-H,O-Z)
       integer(8) n
       dimension wp(2),ps(3),ps2(3)
     end function wperp
     subroutine rround(x,i,dx)
       IMPLICIT double precision (A-H,O-Z)
       integer i
     end subroutine rround
     recursive function pathwalk(node,dir,turns,po,mxpo,td,it) result(lastnode)
       integer,parameter :: tnodes=155000
       integer ia(2),td(-5:5,10),it(tnodes,10),k,j,turns(10),lastnode,node,po,dir,mxpo,nextnode
     end function pathwalk
  end interface
  read(*,*) mcs,efafa,esksk,eve,tem,nsap,npix0,perpsize   ! MCS energy-zz temperature
  npix=2*npix0+1
  erule1=0.0
  ifafa=0
  isksk=0
  irule1=0
  ive=0
  if(abs(efafa)>0.00001)ifafa=1
  if(abs(esksk)>0.00001)isksk=1
  if(abs(eve)>0.00001)ive=1
  print *,"fat-fat ",ifafa,efafa
  print *,"ski-ski ",isksk,esksk
  print *,"v-rule ",ive,eve
  print *,"TEM ",tem
  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)
     ep(1,i)=cos(4.0*pi*(i-1)/5.0)
     ep(2,i)=sin(4.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)
  bap(1,:)=position(iba(1,:),ep)
  bap(2,:)=position(iba(2,:),ep)
  abap(1)=sqrt(bap(1,1)**2+bap(1,2)**2)
  abap(2)=sqrt(bap(2,1)**2+bap(2,2)**2)
  latctr=1 ! lattice centering
  iw(:)=iba(1,:)+iba(2,:)
  do io=1,5
     if(mod(iw(io),2)==1)latctr=0
  enddo
  print *,"Lattice centering : ",latctr
  ps(:)=0.
  ps2(:)=0.
  read(*,*)nn
  do i=1,nn
     read(*,*)iy
     iv=ibound(iy)
     iv=bound5(iy,iba,bar,e)
     pv(i,:)=iv(:)
     p1=ix(iy)
     call perpsum(iv,1,ep,ps,ps2)
!     print *,i,p1," xx  ",pv(i,:)
  enddo
  do i=1,nn
  enddo
  wp=wperp(ps,ps2,nn)  ! initial perp variance
  print "('Initial perp variance ',f8.4)",wp(1)
  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
  enddo
  half=0.499999 ! these are for perp image
  cn(1,1)=-half
  cn(1,2)=-half
  cn(2,1)=-half
  cn(2,2)=half
  cn(3,1)=half
  cn(3,2)=-half
  cn(4,1)=half
  cn(4,2)=half
  dmx=0d0
   ! old way
  tl1=(/9,10,6,7,8,4,5,1,2,3/)  ! turn left (by 1)
  tr1=(/8,9,10,6,7,3,4,5,1,2/)  ! turn right (by 1)
  tl2=(/2,3,4,5,1,7,8,9,10,6/) ! turn left by 2
  tr2=(/5,1,2,3,4,10,6,7,8,9/) ! turn right by 2
  tl3=(/10,6,7,8,9,5,1,2,3,4/)  ! turn left by 3
  tr3=(/7,8,9,10,6,2,3,4,5,1/)  ! turn right by 3
  tl4=(/3,4,5,1,2,8,9,10,6,7/)  ! turn left by 4
  tr4=(/4,5,1,2,3,9,10,6,7,8/)  ! turn right by 4
  tl5=(/6,7,8,9,10,1,2,3,4,5/)  ! turn left by 5  just for completness
  tr5=(/6,7,8,9,10,1,2,3,4,5/)  ! turn right by 5   just for completness
  ! new way
  td(0,:)=(/1,2,3,4,5,6,7,8,9,10/)  ! no turn
  td(-1,:)=(/9,10,6,7,8,4,5,1,2,3/)  ! turn left (by 1)
  td(-2,:)=(/2,3,4,5,1,7,8,9,10,6/) ! turn left by 2
  td(-3,:)=(/10,6,7,8,9,5,1,2,3,4/)  ! turn left by 3
  td(-4,:)=(/3,4,5,1,2,8,9,10,6,7/)  ! turn left by 4
  td(-5,:)=(/6,7,8,9,10,1,2,3,4,5/)  ! turn left by 5  just for completness
  td(1,:)=(/8,9,10,6,7,3,4,5,1,2/)  ! turn right (by 1)
  td(2,:)=(/5,1,2,3,4,10,6,7,8,9/) ! turn right by 2
  td(3,:)=(/7,8,9,10,6,2,3,4,5,1/)  ! turn right by 3
  td(4,:)=(/4,5,1,2,3,9,10,6,7,8/)  ! turn right by 4
  td(5,:)=(/6,7,8,9,10,1,2,3,4,5/)  ! turn right by 5   just for completness
  nnff=0  ! total fat-fat
  nnss=0  ! total ski-ski
  nsmp=0    ! number of samples for perp images
  do i=1,nn
     do j=1,5
        if(it(i,j)>0)then
           k=iwrap(j+5,10)
           if(it(i,k)>0) then  ! two straight edges
              nei=0
              kk=j
              do io=1,4
                 kk=tr1(kk)
                 ii=it(i,kk)
                 if(ii>0)then
                    nei=nei+1
                    mio=io
                 endif
              enddo
              if(nei==1.and.mio==1)nnss=nnss+1
              if(nei==1.and.mio==4)nnss=nnss+1
              if(nei==1.and.mio==2)nnff=nnff+1
              if(nei==1.and.mio==3)nnff=nnff+1
           endif
        endif
     enddo
  enddo
!  nnss=nnss/4
!  nnff=nnff/4
  print *," ski-ski fat-fat init : ",nnss,nnff
  egy=efafa*nnff+esksk*nnss ! initial energy
  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
!           print *,"VERTEX ",n
!           print *,n
!           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
              rp(io)=iwrap(k0+rr(io),10)
           enddo
           do io=1,3
              pp(io)=it(n,rp(io))  ! 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
           nffo=0
           nffn=0
           nsso=0
           nssn=0
           nvo=0
           nvn=0
           nr1o=0
           nr1n=0
           ir1=iwrap(rp(1)+5,10)
           ir2=iwrap(rp(2)+5,10)
           ir3=iwrap(rp(3)+5,10)
           ! HAMILTONIAN : adjacent fat rhombi
!           print *,"ICASE ",icase
           if(icase<3) then ! skinnyhex
              if(irule1==1)then
              endif
              if(ive==1) then ! V-rule
                 k1=tr2(rp(3))   ! turn by +2
                 j1=it(pp(3),k1)
                 if(j1>0) then
                    k2=tr2(k1)  ! turn by +2
                    j2=it(j1,k2)
                    if(j2>0) then
                       j3=it(j2,ir3)
                       if(j3>0)nvo=nvo+1
                    endif
                 endif
                 k1=tl2(rp(3))   ! turn by -2
                 j1=it(pp(3),k1)
                 if(j1>0) then 
                    k2=tl2(k1)  ! turn by -2
                    j2=it(j1,k2)
                    if(j2>0) then
                       j3=it(j2,ir3)
                       if(j3>0)nvo=nvo+1
                    endif
                 endif
                 ! V-rule, AFTER
                 k1=tl2(ir3)   ! turn by -2
                 j1=it(qq(3),k1)
                 if(j1>0) then
                    k2=tl2(k1)  ! turn by -2
                    j2=it(j1,k2)
                    if(j2>0) then
                       j3=it(j2,rp(3))
                       if(j3>0)nvn=nvn+1
                    endif
                 endif
                 k1=tr2(ir3) 
                 j1=it(qq(3),k1)
                 if(j1>0) then
                    k2=tr2(k1)
                    j2=it(j1,k2)
                    if(j2>0) then
                       j3=it(j2,rp(3))
                       if(j3>0)nvn=nvn+1
                    endif
                 endif
              endif ! ive==1
              if(ifafa==1) then ! ADJACENT FAT's
                 j1=it(pp(1),rp(1))
                 j2=it(qq(3),rp(1))
                 if(j1>0.and.j2>0)nffo=nffo+1
                 j3=it(pp(2),rp(2))
                 j4=it(qq(3),rp(2))
                 if(j3>0.and.j4>0)nffo=nffo+1
!!                 print *,' skihex old ',j1,j2,j3,j4
                 ! skihex AFTER
                 j1=it(qq(1),ir2)
                 j2=it(pp(3),ir2)
                 if(j1>0.and.j2>0)nffn=nffn+1
                 j3=it(qq(2),ir1)
                 j4=it(pp(3),ir1)
                 if(j3>0.and.j4>0)nffn=nffn+1
!                 print *,' skihex new ',j1,j2,j3,j4
              endif ! ifafa==1
              if(isksk==1) then
                 nsso1=0
                 nsso2=0
                 nsso3=0
                 nsso4=0
                 nssn1=0
                 nssn2=0
                 nssn3=0
                 nssn4=0
                 if(it(qq(1),rp(1))>0)nsso1=nsso1+1
                 if(it(qq(2),rp(2))>0)nsso2=nsso2+1
                 if(it(qq(1),rp(3))>0)nsso3=nsso3+1
                 if(it(qq(2),rp(3))>0)nsso4=nsso4+1
                 if(it(pp(1),ir2)>0)nssn1=nssn1+1
                 if(it(pp(2),ir1)>0)nssn2=nssn2+1
                 if(it(pp(1),ir3)>0)nssn3=nssn3+1
                 if(it(pp(2),ir3)>0)nssn4=nssn4+1
                 nsso=nsso1+nsso2+nsso3+nsso4
                 nssn=nssn1+nssn2+nssn3+nssn4
!                 print "('qq1 rp1 it',4i4)",qq(1),rp(1),it(qq(1),rp(1)),nsso1
!                 print "('qq2 rp2 it',4i4)",qq(2),rp(2),it(qq(2),rp(2)),nsso2
!                 print "('qq1 rp3 it',4i4)",qq(1),rp(3),it(qq(1),rp(3)),nsso3
!                 print "('qq2 rp3 it',4i4)",qq(2),rp(3),it(qq(2),rp(3)),nsso4
!                 print "('pp1 ir2 it',4i4)",pp(1),ir2,it(pp(1),ir2),nssn1
!                 print "('pp2 ir1 it',4i4)",pp(2),ir1,it(pp(2),ir1),nssn2
!                 print "('pp1 ir3 it',4i4)",pp(1),ir3,it(pp(1),ir3),nssn3
!                 print "('pp2 ir3 it',4i4)",pp(2),ir3,it(pp(2),ir3),nssn4
!                 print *,"nssn nsso ",nssn,nsso
              endif
           endif ! icase<3
           if(icase>2) then ! fathex
              if(irule1==1)then
              endif
              if(ive==1) then ! V-rule from fathex: quite tricky
              endif
              if(ifafa==1) then
                 j1=it(qq(1),rp(3))
                 j2=it(qq(2),rp(3))
                 j3=it(pp(3),rp(3))
                 if(j1>0.and.j3>0)nffo=nffo+1
                 if(j2>0.and.j3>0)nffo=nffo+1
                 j4=it(qq(1),rp(1))
                 j5=it(pp(1),rp(1))
                 if(j4>0.and.j5>0)nffo=nffo+1
                 j6=it(qq(2),rp(2))
                 j7=it(pp(2),rp(2))
                 if(j6>0.and.j7>0)nffo=nffo+1
!              print *,' fathex old ',j1,j2,j3,j4,j5,j6,j7,'  pq ',pp(:),qq(:)
!              print *,qq(1),rp(1),it(qq(1),rp(1))
!              print *,qq(2),rp(2),it(qq(2),rp(2))
              ! fathex AFTER
                 j1=it(pp(1),ir3)
                 j2=it(pp(2),ir3)
                 j3=it(qq(3),ir3)
                 if(j1>0.and.j3>0)nffn=nffn+1
                 if(j2>0.and.j3>0)nffn=nffn+1
                 j4=it(pp(1),ir2)
                 j5=it(qq(1),ir2)
                 if(j4>0.and.j5>0)nffn=nffn+1
                 j6=it(pp(2),ir1)
                 j7=it(qq(2),ir1)
                 if(j6>0.and.j7>0)nffn=nffn+1
!                 print *,'fafa fathex new ',j1,j2,j3,j4,j5,j6,j7
              endif ! ifafa==1
              if(isksk==1) then
                 if(it(qq(3),rp(1))>0)nsso=nsso+1
                 if(it(qq(3),rp(2))>0)nsso=nsso+1
                 if(it(pp(3),ir1)>0)nssn=nssn+1
                 if(it(pp(3),ir2)>0)nssn=nssn+1
!                 print "('qq3 rp1 it',i4,i3,i4)",qq(3),rp(1),it(qq(3),rp(1))
!                 print "('qq3 rp2 it',i4,i3,i4)",qq(3),rp(2),it(qq(3),rp(2))
!                 print "('pp3 ir1 it',i4,i3,i4)",pp(3),ir1,it(pp(3),ir1)
!                 print "('pp3 ir2 it',i4,i3,i4)",pp(3),ir2,it(pp(3),ir2)
!                 print *,"nssn nsso ",nssn,nsso
              endif
           endif  ! icase>2
           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 *,"debug node-list pp qq : ",pp(:),qq(:)
              print *,"debug it : ",it(n,:)
!              print *,"debug k0 ka :  ",icase,k0,ka
!              print *,"debug rr ss : ",rr(:),ss(:)
!              goto 4444
           endif
           de=efafa*(nffn-nffo)+eve*(nvn-nvo)+esksk*(nssn-nsso)+erule1*(nr1n-nr1o)
           call random_number(x)
           if(de>0.00001) then
              expo=exp(-de/tem)
              if(x>expo)goto 5555 ! REJECT !!!
           endif
           mm=mm+1 ! accepted
           nnff=nnff+nffn-nffo
           nnss=nnss+nssn-nsso
           egy=efafa*nnff+esksk*nnss
!           if(modulo(mm,1000)==0)write(98,"(f15.2,2i15)")egy,nnss,nnff
!           if(modulo(mm,1000)==0)write(98,*)e,nnss,nnff
!           if(modulo(mm,1)==0) then
! write(101,"(f15.2,i12,2i5,i12,3i5,'  ',8i3)")de,nnff,nffn,nffo,nnss,nssn,nsso,icase,nssn1,nssn2,nssn3,nssn4,nsso1,nsso2,nsso3,nsso4
!              write(101,"(10i9)")it(pp(2),:)
           !           endif
!           call perpsum(pv(n,:),-1,ep,ps,ps2)
           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
!           call perpsum(pv(n,:),1,ep,ps,ps2)
           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
!           write(101,"(10i9)")it(pp(2),:)
!           if(icase<3)goto 3333
!           goto 3333
5555       continue
        endif  ! end nit(n)==3
     enddo  ! end ii loop over nn nodes
     if(mod(jj,nsap)==0) then  ! test node for wrap-around every nsap MCS
        nsmp=nsmp+1
        ps(:)=0.
        ps2(:)=0.
        do ii=1,nn
           pv(ii,:)=bound5(pv(ii,:),iba,bar,e)
           call perpsum(pv(ii,:),1,ep,ps,ps2)
        enddo
        wp=wperp(ps,ps2,nn)
        write(97,"(2f9.4)")wp(1)
        write(98,"(f15.2)")egy
        ps(:)=0.
        ps2(:)=0.
        do ii=1,nn
           pv(ii,:)=bound5(pv(ii,:),iba,bar,e)
           call perpsum(pv(ii,:),1,ep,ps,ps2)
        enddo
!        do ii=1,nn
!           xh(:)=0.
!           do io=1,5
!              xh(:)=xh(:)+pv(ii,io)*ep(:,io)
!           enddo
!           xh(:)=xh(:)-ps(1:2)/nn
!           write(96,"(2f9.4)")xh
!        enddo
        do ii=1,nn
           xh(:)=0.
           do io=1,5
              xh(:)=xh(:)+pv(ii,io)*ep(:,io)
           enddo
           xh(:)=xh(:)-ps(1:2)/nn
           w1=xh(1)/perpsize*npix
           w2=xh(2)/perpsize*npix
           call rround(w1,m1,dx1)
           call rround(w2,m2,dx2)
           tarea=0d0
           do jo=1,4
              w1=xh(1)/perpsize*npix+cn(jo,1)
              w2=xh(2)/perpsize*npix+cn(jo,2)
              call rround(w1,k1,dk1)
              call rround(w2,k2,dk2)
              if(k1>-pdim-1.and.k1<pdim+1.and.k2>-pdim-1.and.k2<pdim+1) then
                 dm1=(k1-m1+dk1)
                 dm2=(k2-m2+dk2)
                 if(dm1.lt.0)dm1=-dm1
                 if(dm2.lt.0)dm2=-dm2
                 area=(1-dm1)*(1-dm2)
                 d(k1,k2)=d(k1,k2)+area
                 tarea=tarea+area
                 if(d(k1,k2).gt.dmx)dmx=d(k1,k2)
              endif
           enddo
        enddo
     endif
  enddo  ! end jj MCS loop
3333 continue
  x1=mm
  x2=mcs*nn
  if(nsap>0) then
     af=abap(1)*abap(2)/(perpsize/npix)**2
     open(95,file='perp.info')
     write(95,"(i6,'   pixels')")npix0
     write(95,"(f6.2,'   perpsize')")perpsize
     write(95,"(f7.4,'   maxden')")dmx/nsmp*af
     close(95)
     do i1=-npix0,npix0
        do i2=-npix0,npix0
           write(96)d(i1,i2)/nsmp*af
        enddo
     enddo
  endif
  print "('tem acc : ',f8.2,f8.4)",tem,x1/x2
  print *," ski-ski fat-fat BEFORE CHECK : ",nnss,nnff
  nnff=0  ! total fat-fat
  nnss=0  ! total ski-ski
  do i=1,nn
     do j=1,5   ! HALF is enough to have each once!
        if(it(i,j)>0)then
           k=iwrap(j+5,10)
           if(it(i,k)>0) then  ! two straight edges
              nei=0
              kk=j
              do io=1,4   ! we only check 180 degrees to the right - to have each eactly once!
                 kk=tr1(kk)
                 ii=it(i,kk)
                 if(ii>0)then
                    nei=nei+1
                    mio=io
                 endif
              enddo
              if(nei==1.and.mio==1)nnss=nnss+1
              if(nei==1.and.mio==4)nnss=nnss+1
              if(nei==1.and.mio==2)nnff=nnff+1
              if(nei==1.and.mio==3)nnff=nnff+1
           endif
        endif
     enddo
  enddo
!  nnss=nnss/4
!  nnff=nnff/4
  print *," ski-ski fat-fat CHECK : ",nnss,nnff
  write(99,"(5i8)")iba(1,:)
  write(99,"(5i8)")iba(2,:)
  write(99,"(i10,'     1')")nn
  do i=1,nn
     write(99,"(5i8,'         0')")pv(i,:)
  enddo
  write(99,"('ctable')")
  do i=1,nn
     write(99,"(10i10)") it(i,:)
  enddo
4444 continue
end program Main
      
