      PROGRAM invasion
*     侵略型パーコレーション・クラスターの生成
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER IR_, Lx, Ly, perx(5000), pery(5000), nper
      CALL GWopen(IR_, 0)
      CALL initial(Lx,Ly)
      CALL assign(site,perx,pery,nper,Lx,Ly)
      CALL invade(site,perx,pery,nper,Lx,Ly)
      CALL average(site,Lx,Ly)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(Lx,Ly)
      IMPLICIT NONE
      REAL xwin, ywin, GWaspect
      REAL dmy, rnd
      INTEGER IR_, x, y, Lx, Ly
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'length in y direction = '
      READ(*,*) Ly
      Lx = 2*Ly
      CALL GWvport(IR_, 0.01*GWaspect(1), 0.2, 0.65*GWaspect(1), 0.99)
      CALL compute_aspect_ratio(REAL(Lx),xwin,ywin)
      CALL GWindow(IR_, -0.5, -0.5, xwin-0.5, ywin-0.5)
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      CALL GWrect(IR_, 1.0, 1.0, Lx+1.0, Ly+1.0)
      DO y = 1, Ly
         DO x = 1, Lx
            CALL GWsetpxl(IR_, REAL(x)+0.5, REAL(y)+0.5, 0)
         END DO
      END DO
      END
C
      SUBROUTINE assign(site,perx,pery,nper,Lx,Ly)
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      REAL rnd
      INTEGER IR_, x, y, nper, perx(5000), pery(5000), Lx, Ly
      DO y = 1, Ly
         site(1,y) = 1                      ! 最初の列を占有する
         CALL GWsrect(IR_, 1.0, REAL(y), 
     +                     2.0, REAL(y)+1.0, 16)
      END DO
*     残った格子点に乱数を割り当てる
      DO y = 1, Ly
         DO x = 2, Lx
            site(x,y) = rnd(0)
         END DO
      END DO
*     第2列にある格子点が最初の周辺の点となる
*     周辺の点であるなら site(x,y) は 2 より大きい
      x = 2
      nper = 0
      DO y = 1, Ly
         site(x,y) = 2 + site(x,y)
         nper = nper + 1                    ! 周辺の点の数
*        周辺の点をソートして周辺の点の表の順序付けを行う
*       
         CALL insert(site,perx,pery,nper,x,y)
      END DO
      END
C
      SUBROUTINE insert(site,perx,pery,nper,x,y)
*     線形探索または2分法探索のルーチンを呼ぶ
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER ilist, ninsert, nper, x, y, perx(5000), pery(5000)
      CALL binary_search(site,perx,pery,nper,x,y,ninsert)
*     小さい乱数の格子点のすべてをつぎの配列要素へ移す
      DO ilist = nper, ninsert + 1, - 1
         perx(ilist) = perx(ilist-1)
         pery(ilist) = pery(ilist-1)
      END DO
      perx(ninsert) = x                     ! 新しい格子点を表に挿入する
      pery(ninsert) = y
      END
C
      SUBROUTINE binary_search(site,perx,pery,nper,x,y,ninsert)
*     表を半分に分けそのどちらにその乱数が属するかを調べ，
*     その乱数の位置が定まるまでこの分割を続ける
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER nmid, xmid, ymid, nper, x, y, ninsert, nfirst, nlast, 
     +        perx(5000), pery(5000)
      LOGICAL Ltmp1_
      nfirst = 1                            ! 表の始め
      nlast = nper - 1                      ! 表の終わり
      IF(nlast .LT. 1) nlast = 1
      nmid = (nfirst + nlast)/2             ! 表の中央
*     新しい数が表のどちらに属するかを決める
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         IF(nlast - nfirst .LE. 1) THEN
            ninsert = nlast
            RETURN
         END IF
         xmid = perx(nmid)
         ymid = pery(nmid)
         IF(site(x,y) .GT. site(xmid,ymid)) THEN
            nlast = nmid                    ! 上半分を検索
         ELSE
            nfirst = nmid                   ! 下半分を検索
         END IF
         nmid = (nfirst + nlast)/2
      END DO
      END
C
      SUBROUTINE linear_search(site,perx,pery,nper,x,y,ninsert)
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER iper, xperim, yperim, nper, x, y, ninsert,
     +        perx(5000), pery(5000)
      IF(nper .EQ. 1) THEN
         ninsert = 1
      ELSE
         DO iper = 1, nper - 1
            xperim = perx(iper)
            yperim = pery(iper)
            IF(site(x,y) .GT. site(xperim,yperim)) THEN
               ninsert = iper               ! 新しい格子点を挿入
               RETURN
            END IF
         END DO
      END IF
      ninsert = nper
      END
C
      SUBROUTINE invade(site,perx,pery,nper,Lx,Ly)
*     nx と ny は最隣接点に向かうベクトル成分
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER IR_, nx(4), ny(4), perx(5000), pery(5000),
     +       Lx, Ly, i, x, xnew, y, ynew, nper, DP_, DATA_(8)
      LOGICAL Ltmp1_
      DATA DP_/1/, DATA_/1,0,-1,0,0,1,0,-1/
      DO i = 1, 4
         nx(i) = DATA_(DP_)
         DP_ = DP_ + 1
         ny(i) = DATA_(DP_)
         DP_ = DP_ + 1
      END DO
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         x = perx(nper)
         y = pery(nper)
         nper = nper - 1
*        占有された点に印をつける．その点は周辺の点ではなくなる
         site(x,y) = site(x,y) - 1
         CALL GWsrect(IR_, REAL(x), REAL(y), 
     +                     REAL(x)+1.0, REAL(y)+1.0, 16)
         DO i = 1, 4
            xnew = x + nx(i)
            ynew = y + ny(i)
            IF(ynew .GT. Ly) THEN
               ynew = 1
            ELSE IF(ynew .LT. 1) THEN
               ynew = Ly
            END IF
            IF(site(xnew,ynew) .LT. 1) THEN
               site(xnew,ynew) = site(xnew,ynew) + 2
               nper = nper + 1
               CALL insert(site,perx,pery,nper,xnew,ynew)
            END IF
         END DO
         Ltmp1_ = .NOT.(x .GE. Lx)
      END DO
      END
C
      SUBROUTINE average(site,Lx,Ly)
*     確率密度 P(r) を計算する
      IMPLICIT NONE
      REAL*8 site(0:200,0:100), dr, r, P(0:20)
      INTEGER Lmax, Lmin, n, nbin, nr(0:20), occupied,
     +       ibin, x, y, Lx, Ly
      DATA occupied/0/, P/21*0.0D0/, nr/21*0/
      Lmin = Lx/3
      Lmax = 2*Lmin
      n = (Lmax - Lmin + 1)*Ly              ! 格子の中央の半分にある格子点数
      dr = 0.05D0
      nbin = 1/dr
      DO x = Lmin, Lmax
         DO y = 1, Ly
            ibin = nbin*(mod(site(x,y),1))
            nr(ibin) = nr(ibin) + 1
            IF(site(x,y) .GE. 1 .AND. site(x,y) .LT. 2) THEN
               occupied = occupied + 1      ! 占有された格子点の総数
               P(ibin) = P(ibin) + 1
            END IF
         END DO
      END DO
      WRITE(*,*) '# occupied sites =', occupied
      WRITE(*,'(2A13)') ' r', '  P(r)'
      WRITE(*,*)
      DO ibin = 0, nbin
         r = dr*ibin
         IF(nr(ibin) .GT. 0) WRITE(*,'(2F13.6)') r, P(ibin)/nr(ibin)
      END DO
      WRITE(*,*)
      END
