!=========================================
!　　　ルンゲクッタ　メインプログラム
!　　　　・問題の設定をこのxrunge.f90で行う。
!　　　　・ルンゲクッタステップは
!　　　　　サブルーチン"runge"として
!　　　　　runge.f90に記述されている。
!　　　　　メインプログラム内で呼び出して使う。
!=========================================
program main
  implicit none !暗黙の型宣言禁止

  integer(kind=4),parameter:: neq=1 !ここでは方程式の本数を１に
  integer(kind=4),parameter:: nmax=5000 !タイムステップ数　ここでは5000に

  real(kind=8),parameter:: eps=1.0d-10

  external ex2 !関数を宣言。下の方にある。

  real(kind=8),dimension(neq):: y0,yn !Nに相当。ynにはrungeで更新された値が入る。
  real(kind=8) x0,xn,dx !タイムステップ関係を定義

  real(kind=8),dimension(neq,12):: work

  integer(kind=4) ierr,i,n,iseed
  character(len=20) filename !データファイル名を格納する場所
!------------------　ここまで宣言　-----------------------

  n=0
  x0=0.0d0 !ステップの初期値
  dx=1.0d-2 !時間刻み幅。タイムステップ数に合わせて変更。
  y0(1)=1.0d5 !yの初期値。いまは１０万人としている。
!----------------　ここまで初期値設定　-------------------------

  write(*,*) "filename ?" !標準出力（PC画面）にコレを表示
  read(*,*) filename !ユーザーが打ち込んだ文字列をfilenameに代入

!-------------　ここまで標準入出力操作　----------------------------

  open(10,file=filename) !filenameを装置番号"10"として開く
  write(*,'(i5,3f20.10)') n,x0,y0(1) !初期値を標準出力に表示
  write(10,'(i5,3f20.10)') n,x0,y0(1) !初期値を書き出し
  
  do n=1,nmax !整数変数nを1からnmaxまで変えつつループ開始
     xn=x0+dx !ステップ幅dxで時間を進める
     call runge(neq,ex2,x0,xn,y0,yn,eps,ierr,work) !ルンゲクッタ呼び出し
     x0=xn !
     do i=1,neq !ルンゲによって更新されたyをy0に代入
        y0(i)=yn(i) !rungeで得られたynをy0に代入。iは式の区別。
     end do
     write(*,'(i5,3f20.10)') n,x0,y0(1) !更新値を標準出力に表示。
     write(10,'(i5,3f20.10)') n,x0,y0(1) !更新値を書き出し。
  end do
  close(10) !開いた装置番号"10"を閉じる。
!---------------- ここまで　計算とファイル書き出し -------------------------
  
end program main









!=========================================
!　　　関数　"ex2"
!　　　　・常微分方程式
!          dy/dx = F
!　　　　　の右辺を記述。
!　　　　・サブルーチンrunge内で使われる。
!　　　　・連立の例は他の人のを参照のこと。
!=========================================
real(kind=8) function ex2(i,x,y) !関数を記述
  implicit none !暗黙の型宣言禁止
  
  integer(kind=4),parameter:: neq=1 !方程式の個数
  
  integer(kind=4) i !
  real(kind=8) x
  real(kind=8),dimension(neq):: y

  real(kind=8) fmon,fmot,B,q
  real(kind=8),parameter:: alpha = 0.04d0
  real(kind=8),parameter:: Nmax = 1.0d6
  real(kind=8),parameter:: q0 = 1.0d4
  real(kind=8),parameter:: c1 = 1.1*q0/(1.0d5)**2
!-------------　ここまで宣言　----------------------------
  
    fmon = c1 * (1.0d0-alpha) * y(1)
    fmot = (1.0d0-alpha) * y(1)
    B = (Nmax - y(1)) / Nmax
      q = fmon * fmot
        ex2 = (q-q0)*b
!---------------　ここまで関数--------------------------
  return
end function ex2









