*********************************************
*     etude.for 自励系微分方程式の解曲線群  *
*     (c) Kenrou Adachi, 2004/11/15         *
*     fl /FPi etude.for /link graphics.lib  *
*********************************************

$notruncate

****** 主ルーチン *****

      program etude

      real a, b, eps, e, h, sh
      integer n, l, i

*     *** データの入力 ***
*     *** call inputdata(a, b, eps, n, l) ***

      a = 12.8
      b = 6.4
      eps = 0.2
      n = 3
      l = 100

*     *** グラフィックスシステムの初期化 ***

      call ginit

*     *** 座標軸の描画 ***

      call drawaxes

*     *** 時間刻み幅の設定 ***

      sh = (a / 320.0) * (a /320.0) + (b / 200.0) * (b / 200.0)
      h = 2.0 * sqrt(sh)

*     *** 解曲線群描画 ***

      do 10 i = 1, n
         call drawing(a, b, eps, h, l)
 10   continue

*     *** 処理終了 ***

      write(*,*) '終了しました。'
      call fkeyon
      call cursoron

      stop
      end


***** 解曲線群描画ルーチン *****

      subroutine drawing(a, b, eps, h, l)
      real    a, b, eps, h
      integer l

      integer j, count
      real    e1, e2, e, rnd, x(40), y(40), xx(40), yy(40), rand

*     *** 初期点をランダムに生成 ***

      j = 1
      do 20 e1 = -1.0, 0.9, 0.25
         do 30 e2 = -1.0, 0.9, 0.4
            rnd = rand()
            xx(j) = a * (e1 + 0.25 * rnd)
            rnd = rand()
            yy(j) = b * (e2 + 0.4 * rnd)
            j = j + 1
 30      continue
 20   continue

*     *** 40 本の積分曲線を描く ***

*     * 時間の向きを正に *

      e = 1.0

 40   continue

        do 50 j = 1, 40
           x(j) = xx(j)
           y(j) = yy(j)
 50     continue

        if (e .gt. 0) then
           call setcolor(4)
        else
           call setcolor(2)
        endif

        do 60 count = 1, l
           call repeat(x, y, a, b, eps, h, e)
 60     continue

*     * 時間の向きを逆転してもう一度 *

        e = -e

*     * 逆転は一度だけ *

      if (e .lt. 0) goto 40

      return
      end


***** 微小積分曲線素描画ルーチン *****

      subroutine repeat(x, y, a, b, eps, h, e)
      real x(40), y(40), a, b, eps, h, e

      integer j, iwx, iwy, iwx1, iwy1
      real    x0, y0, xr0, x1, y1, xr1, wx, wy, wx1, wy1
      real    u, v, u1, v1, q, qh

      do 70 j = 1, 40
         x0 = x(j)
         y0 = y(j)

         if ((abs(x0) .gt. (2.0 * a) .and. (abs(y0) .gt. (2.0 * b)) goto 70

         xr0 = sqrt(x0 * x0 + y0 * y0)

         if (xr0 .lt. h) goto 70

         u = fx(x0, y0, xr0)
         v = gy(x0, y0, xr0)
         q = sqrt(u * u + v * v)

         if (q .lt. eps) goto 70

         qh = h / q
         x1 = x0 + 0.5 * e * u * qh
         y1 = y0 + 0.5 * e * v * qh
         xr1 = sqrt(x1 * x1 + y1 * y1)
         u1 = fx(x1, y1, xr1)
         v1 = gy(x1, y1, xr1)

         x(j) = x0 + e * u1 * qh
         y(j) = y0 + e * v1 * qh

         wx = 320.0 * x0 / a
         wy = 200.0 * y0 / b
         wx1 = 320.0 * x(j) / a
         wy1 = 200.0 * y(j) / b

         iwx = int(wx)
         iwy = int(wy)
         iwx1 = int(wx1)
         iwy1 = int(wy1)

         call line(iwx + 320, 200 - iwy, iwx1 + 320, 200 - iwy1)

 70   continue

      return
      end


***** 座標軸描画ルーチン *****

      subroutine drawaxes

      call setcolor(1)
      call line(0, 200, 639, 200)
      call line(320, 0, 320, 399)

      return
      end


***** 自励系微分方程式の定義 *****
***** dx/dt = fx(x, y, r), dy/dt = gy(x, y, r) *****

      real function fx(x, y, r)
      real x, y, r

      fx = cos(x) - x * sin(y)

      return
      end


      real function gy(x, y, r)
      real x, y, r

      gy = x * sin(y) + y * sin(x)

      return
      end


***** 疑似乱数の発生 *****
***** 大駒誠一 著、改訂 FORTRAN 77、サイエンス社 より *****

      real function rand()

      integer c, k, r, m
      real    MAXR
      parameter (c = 656329, k = 163, m = 12518383, MAXR = 12518383.0)
      save r
      data r/12345678/

      r = mod(k * r + c, m)
      rand = real(r) / MAXR

      return
      end