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
function fraco(r,b) result(a)
  IMPLICIT double precision (A-H,O-Z)
  dimension b(2,2),r(2),a(2),e(2,5)
  xb=b(1,1)*b(2,2)-b(1,2)*b(2,1)
  a(1)=(r(1)*b(2,2)-r(2)*b(2,1))/xb
  a(2)=(r(2)*b(1,1)-r(1)*b(1,2))/xb
end function fraco
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
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
        
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(25)   ! list of nodes needed for flip update
  integer p1,p2,q1,q2,ll,ss
  integer it(130000,10)
  integer,parameter :: seed = 86456, tnodes=155000
  dimension e(2,5),iba(2,5)
  dimension bar(2,2),cell(4,2),bbox(4)  ! basis in real space, and cell corners position, BoundingBox
  integer, dimension(tnodes,5) :: pv
  dimension r2(2),s1(2),s2(2),s3(2),u2(2),bco(2)
  interface
     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 ix(ii) result (a)
       implicit none
       integer, dimension(5) :: ii
       integer :: i,is,a
     end function ix
     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 fraco(r,b) result(a)
       IMPLICIT double precision (A-H,O-Z)
       dimension b(2,2),r(2),a(2),e(2,5)
     end function fraco
     recursive function ibound(ii) result(v)
       integer, dimension(5) :: ii,v
       integer :: is=0
     end function ibound
  end interface
  mcs=100
  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
  read(*,*)iba(1,1:5)
  read(*,*)iba(2,1:5)
  bar(1,:)=position(iba(1,:),e)
  bar(2,:)=position(iba(2,:),e)
  cell(1,:)=0
  cell(2,:)=bar(1,:)
  cell(3,:)=bar(1,:)+bar(2,:)
  cell(4,:)=bar(2,:)
  xmin=10000
  xmax=-10000
  ymin=10000
  ymax=-10000
  do io=1,4
     if(cell(io,1)<xmin)xmin=cell(io,1)
     if(cell(io,2)<ymin)ymin=cell(io,2)
     if(cell(io,1)>xmax)xmax=cell(io,1)
     if(cell(io,2)>ymax)ymax=cell(io,2)
  enddo
  bco(1)=xmin ! left-bottom cell corner
  bco(2)=ymin
  xx=xmax-xmin
  psw=600/xx
  do io=1,4
     cell(io,:)=(cell(io,:)-bco(:))*psw
  enddo
  over=0.03*psw
  bbox(1)=-over
  bbox(2)=-over
  bbox(3)=psw+over
  bbox(4)=psw*(ymax-ymin)/xx+over
  print "('%!PS-Adobe-3.0')"
  print "('%%BoundingBox: ',4f6.1)",bbox*xx
  print "('/rgb {setrgbcolor fill} bind def')"
  print "('/L {moveto lineto stroke} bind def')"
  print "('/w {setlinewidth} bind def')"
  print "('/fat {1 0 0 rgb ',f3.1,' 0 360 arc closepath stroke 0 0 0 rgb} bind def')",psw/8.
  print "('/ski {0 0 1 rgb ',f3.1,' 0 360 arc closepath stroke 0 0 0 rgb} bind def')",psw/8.
  print "('/box{moveto lineto lineto lineto closepath} def')"
  print "('1.0 w 0 1 0 rgb ',8f8.2,'  box stroke 0.1 w')",cell(1,:),cell(2,:),cell(3,:),cell(4,:)
  print *,"0 0 0 rgb"
  read(*,*)nn
  do i=1,nn
     read(*,*)iy
     iv=bound5(iy,iba,bar,e)
     iv=ibound(iv)
     pv(i,:)=iv(:)
     p1=ix(iy)
  enddo
  read(*,*)
  do i=1,nn
     read(*,*) it(i,:)
  enddo
  open(101,file='coord.fat')
  open(102,file='coord.ski')
  open(103,file='coord')
  do i=1,nn
     r2=position(pv(i,:),e)-bco ! vertex position
     do k=1,10
        q1=it(i,k)
        if(q1>0) then
           if(i<q1) then ! draw edge, cycle neighb. cells
              if(k<6)s2(:)=e(:,k)
              if(k>5)s2(:)=-e(:,k-5)
              do i1=-1,1
                 do i2=-1,1
                    u2=i1*bar(1,:)+i2*bar(2,:)
                    print "(4f10.2,'  L')",(r2+u2)*psw,(r2+s2+u2)*psw
                 enddo
              enddo
           endif
           k1=k-1 ! start fat : from 72 deg corner
           if(k==1.or.k==6)k1=k1+5
           q2=it(i,k1)
           if(q2>0) then !
              q3=it(q2,k) 
              if(i<q3.and.q3>0) then  ! we have fat
                 if(k<6)s1(:)=e(:,k)/2. ! displacement half-edge
                 if(k>5)s1(:)=-e(:,k-5)/2.
                 if(k1<6)s1(:)=s1(:)+e(:,k1)/2. ! add displace. half other edge
                 if(k1>5)s1(:)=s1(:)-e(:,k1-5)/2.
                 s1=r2+s1   ! center of skinny
!                 s3=fraco(s1,bar)
!                 if(s3(1)<-1.0001)s1=s1+bar(1,:)
!                 if(s3(1)>0.0001)s1=s1-bar(1,:)
!                 if(s3(2)<-1.0001)s1=s1+bar(2,:)
!                 if(s3(2)>0.0001)s1=s1-bar(2,:)
                 print "(2f10.2,'  fat')",s1*psw
                 write(101,"(2f17.8)")s1
                 write(103,"(2f17.8,'  1')")s1
              endif
           endif  ! end fat
           k1=k-2 ! start skinny: from 144 deg corner
           if(k<6.and.k1<1)k1=k1+5
           if(k>5.and.k1<6)k1=k1+5
           q2=it(i,k1)
           if(q2>0) then !
              q3=it(q2,k) 
              if(i<q3.and.q3>0) then  ! we have fat
                 if(k<6)s1(:)=e(:,k)/2. ! displacement half-edge
                 if(k>5)s1(:)=-e(:,k-5)/2.
                 if(k1<6)s1(:)=s1(:)+e(:,k1)/2. ! add displace. half other edge
                 if(k1>5)s1(:)=s1(:)-e(:,k1-5)/2.
                 s1=r2+s1   ! center of skinny
!                 s3=fraco(s1,bar)
!                 if(s3(1)<-1.0001)s1=s1+bar(1,:)
!                 if(s3(1)>0.0001)s1=s1-bar(1,:)
!                 if(s3(2)<-1.0001)s1=s1+bar(2,:)
!                 if(s3(2)>0.0001)s1=s1-bar(2,:)
                 print "(2f10.2,'  ski')",s1*psw
                 write(102,"(2f17.8)")s1
                 write(103,"(2f17.8,'  2')")s1
              endif
           endif  ! end skinny
        endif        
     enddo
  enddo
4444 continue
end program Main
      
