--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rkmod.f90	Sun Jul 21 03:56:49 2024 +0900
@@ -0,0 +1,39 @@
+module rkmod
+  use rungekutta !ルンゲクッタ解法をまとめたモジュール
+  implicit none
+  type(rk) :: rktype  !ルンゲクッタ法ソルバに関連する型
+  !========================================
+  contains
+    !微分方程式の右辺を定義
+    function rhside(self,time,wk)
+      class(rk),intent(inout) :: self !parameter
+      real(kind=dpkind) :: time
+      real(kind=dpkind),dimension(:) :: wk
+      real(kind=dpkind),dimension(size(wk)) :: rhside
+      real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
+      rhside(1)=a*(wk(2)-wk(1))
+      rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
+      rhside(3)=wk(1)*wk(2)-c*wk(3)
+    end function
+    !Rから呼び出すサブルーチン
+    subroutine calc(x, max_t, r)
+      implicit none
+      real(kind=dpkind),dimension(3),intent(in) :: x !状態変数
+      real(kind=dpkind),intent(in) :: max_t
+      real(kind=dpkind),allocatable,dimension(:,:),intent(out) :: r
+      real(kind=dpkind) :: zero = 0
+      integer i, max_i
+      !微分方程式の右辺を計算する関数、解法、時間刻みの指定
+      call rktype%setup(rhside,"rk4",1d-5) 
+      max_i = int(max_t / 0.01)
+      allocate(r(0:max_i, 4))
+      !時刻max_t秒まで、0.01秒毎に配列に記録。
+      !ただし、内部的には刻み幅1e-5秒としている。
+      r(0, :) = [zero, x] ! 開始点
+      do i = 1, max_i
+        call rktype%solver(rktype%time+1e-2,x) !積分実行
+        r(i, :) = [rktype%time, x]
+      end do
+    end subroutine
+end module
+