!!! creates tiling+connectivity table for a plain tailing input.
!!! Usage: cat TILING.IN | ./pttable > tiling.ctable

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
  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"
  do i=1,nn-1
     do j=i+1,nn
        iu(:)=pv(j,:)
        iv(:)=pv(i,:)
        iw=iu-iv
        iw=iboundapp(iw,iba(1,:),iba(2,:))
        k=ipo(iw)  ! k>0 only when difference between two nodes is ONE
        if(k>0)then
           it(i,k)=j
           if(k<6)m=k+5
           if(k>5)m=k-5
           it(j,m)=i
!           print *,i,j,k,m," yy",iw
        endif
     enddo
  enddo

3333 continue
  print *, iba(1,:)
  print *, iba(2,:)
  print *, nn," 1"
  do i=1,nn
     print *, pv(i,:),"     0"
  enddo
  print *,"ctable"
  do i=1,nn
     print *, it(i,:)
  enddo
4444 continue
end program Main
      
