# HG changeset patch
# User "uncorrelated zombie" <uncorrelated@yahoo.co.jp>
# Date 1721669457 -32400
# Node ID 5bd0f2a15d2a2192efdcd4ddd7e93ef8dd52ed4c

import from https://qiita.com/sakamoti/items/de851e0d07aeef5be310

diff -r 000000000000 -r 5bd0f2a15d2a .hgignore
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Tue Jul 23 02:30:57 2024 +0900
@@ -0,0 +1,4 @@
+syntax: glob
+*.mod
+*.o
+
diff -r 000000000000 -r 5bd0f2a15d2a rkmain.f90
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rkmain.f90	Tue Jul 23 02:30:57 2024 +0900
@@ -0,0 +1,31 @@
+program rkmain
+  use rungekutta !ルンゲクッタ解法をまとめたモジュール
+  implicit none
+  real(kind=dpkind),dimension(3) :: x !状態変数
+  type(rk) :: rktype  !ルンゲクッタ法ソルバに関連する型
+
+  !微分方程式の右辺を計算する関数、解法、時間刻みの指定
+  call rktype%setup(rhside,"rk4",1d-5) 
+  !初期値の指定
+  x=[0.5e0,0.5e0,0.5e0]
+  !時刻0?50秒まで、0.01秒毎に標準出力へ書き出し。
+  !ただし、内部的には刻み幅1e-5秒としている。
+  do 
+    if(rktype%time>5e1)exit
+    call rktype%solver(rktype%time+1e-2,x) !積分実行
+    print*,rktype%time,x
+  enddo
+  !========================================
+  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
+end program
diff -r 000000000000 -r 5bd0f2a15d2a rungekutta.f90
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rungekutta.f90	Tue Jul 23 02:30:57 2024 +0900
@@ -0,0 +1,100 @@
+module rungekutta
+  use,intrinsic :: iso_fortran_env
+  private
+  !---
+  integer,parameter,public :: dpkind=real64 
+  real(kind=dpkind),public,parameter :: PI=3.1415926535897932384626433_dpkind
+  !---
+  type,public :: rk
+    real(kind=dpkind) :: time=0e0 !<險育ｮ鈴幕蟋九°繧峨�ｮ邱冗ｵ碁℃譎る俣
+    real(kind=dpkind) :: dt=1e-4   !<譎る俣蛻ｻ縺ｿ
+    procedure(rhside),pass,pointer :: fun => null() 
+                !<蠕ｮ蛻�譁ｹ遞句ｼ上�ｮ蜿ｳ霎ｺ繧定ｨ育ｮ励☆繧矩未謨ｰ 
+    procedure(rksol),pass,pointer :: solver => null()      
+                !<蠕ｮ蛻�譁ｹ遞句ｼ上�ｮ隗｣豕�
+    contains
+      final :: final_rk !<讒矩��菴薙�ｮ繝輔ぃ繧､繝翫Λ繧､繧ｺ蜃ｦ逅�
+      procedure :: setup !<謨ｰ蛟､遨榊��縺ｮ蛻晄悄蛟､縲∬ｧ｣豕輔ｒ繧ｻ繝�繝医☆繧区焔邯壹″
+  end type
+  !-------------
+  interface
+    !---
+    ! Right hand side function interface
+    function rhside(self,time,wk)
+      import :: dpkind,rk
+      class(rk),intent(inout) :: self !parameter
+      real(kind=dpkind) :: time
+      real(kind=dpkind),dimension(:) :: wk
+      real(kind=dpkind),dimension(size(wk)) :: rhside
+    end function
+    !---
+    ! Solver interface
+    subroutine rksol(self,tnext,wk)
+      import :: dpkind,rk
+      class(rk),intent(inout) :: self
+      real(kind=dpkind),intent(in) :: tnext !output time
+      real(kind=dpkind),dimension(:) :: wk
+    end subroutine
+  end interface
+  !============================--
+  contains
+    !---
+    pure elemental subroutine final_rk(self)
+      type(rk),intent(inout) :: self
+      if(associated(self%fun))nullify(self%fun)
+      if(associated(self%solver))nullify(self%solver)
+    end subroutine
+    !---
+    subroutine setup(self,fun,solver,dt)
+      class(rk),intent(inout) :: self
+      procedure(rhside) :: fun  !蠕ｮ蛻�譁ｹ遞句ｼ上�ｮ蜿ｳ霎ｺ繧定ｨ育ｮ励☆繧矩未謨ｰ
+      character(len=*),intent(in) :: solver !繧ｽ繝ｫ繝舌�ｼ縺ｮ驕ｸ謚�
+      real(kind=dpkind),intent(in),optional :: dt !譎る俣蛻ｻ縺ｿ蟷� 
+      !譎る俣蛻ｻ縺ｿ縺ｮ繧ｻ繝�繝�
+      if(present(dt))then
+        self%dt=dt
+      endif
+      self%fun=>fun
+      !繧ｽ繝ｫ繝舌�ｼ縺ｮ驕ｸ謚�
+      select case(solver)
+      case("rk4")
+        !繝ｫ繝ｳ繧ｲ繧ｯ繝�繧ｿ豕�
+        self%solver=>rk4
+      case("rkf45")
+        !繝ｫ繝ｳ繧ｲ繧ｯ繝�繧ｿ繝ｻ繝輔ぉ繝ｫ繝吶Ν繧ｰ豕�(蛻ｻ縺ｿ蟷�閾ｪ蜍戊ｪｿ謨ｴ蝙�)
+        !self%solver=>rkf45
+      case default
+        write(error_unit,*) "solver must one of the 'rk4, rkf45'"
+      endselect
+    end subroutine
+    !---
+    subroutine rk4(self,tn,wk)
+      class(rk),intent(inout) :: self
+      real(kind=dpkind),intent(in) :: tn !谺｡縺ｫ蜃ｺ蜉帙＠縺溘＞譎る俣
+      real(kind=dpkind),dimension(:) :: wk
+      real(kind=dpkind),dimension(size(wk)) :: k1,k2,k3,k4
+      real(kind=dpkind) :: dt_tmp
+      real(kind=dpkind) :: dthalf
+      associate(time=>self%time,dt=>self%dt)
+        dt_tmp=dt
+        do
+          if(time+dt>tn)dt_tmp=tn-time
+          dthalf=dt_tmp/2e0
+          !!$OMP workshare
+          k1=self%fun(time,wk)
+          k2=self%fun(time+dthalf,wk+k1*dthalf)
+          k3=self%fun(time+dthalf,wk+k2*dthalf)
+          k4=self%fun(time+dt_tmp,wk+k3*dt_tmp)
+          !迥ｶ諷句､画焚繧呈峩譁ｰ
+          wk=wk+dt_tmp*(k1+2e0*k2+2e0*k3+k4)/6e0
+          !!$OMP end workshare
+          !1繧ｹ繝�繝�繝玲ｬ｡縺ｮ隗｣縺瑚ｨ育ｮ励〒縺阪◆縺ｮ縺ｧ縲∵凾髢薙ｒ騾ｲ繧√ｋ
+          time=time+dt_tmp 
+          !譎る俣縺ｮ險育ｮ礼ｯ�蝗ｲ縺荊n繧定ｶ�縺医◆繧芽ｨ育ｮ礼ｵゆｺ�
+          if(time>=tn)then
+            exit
+          endif
+        enddo
+      end associate
+    end subroutine
+end module