c  ルンゲ・クッタ法による1次元運動方程式の解法
c      速度に比例する抵抗がある調和振動子
      real t,x,v,dt
      real t0,x0,v0
      real tlim
      real alpha
      real pi
      integer n,nwrite, i
c  t：時間、x：位置、v：速度、dt：時間の刻み幅
c  t0、x0、v0：初期値
c  tlim：時間の最大値
c  alpha：速度に比例する抵抗の係数
c
      common /para/alpha
c
c  出力ファイルの指定
      open(unit=2,file='harm.dat')
c  初期値
      t0=0.0                                                            
      x0=0.0
      v0=1.0
c  速度に比例する抵抗の係数
      write(6,*) ' alpha ?'
      read(5,*) alpha
c  時間の刻み幅
      write(6,*) ' dt ?'
      read(5,*) dt
c
      pi=3.1415926535897932385
      tlim=6.0*pi
      n=nint(tlim/dt)
      nwrite=n/10
      if(nwrite.eq.0) nwrite=1
c
      i=0
      t=t0                                                              
      x=x0                        
      v=v0                        
   10 continue
      i=i+1
c  ルンゲ・クッタ法ルーチンを呼び出す
      call mtn1d(t,x,v,dt)                                               
c  ファイルに結果を書き出す
      write(2,601) t,x,v                                     
      if(mod(i,nwrite).eq.0) then
c  画面に結果を表示
        write(6,601) t,x,v                                   
      endif
      if(t.gt.tlim) then
        write(6,*) ' Number of Steps= ',i
        stop
      endif
c 次のステップへ
      go to 10  
c                                                        
  601 format(1h ,5(1pe12.5,2x))                                         
c
      end                                                               
*********************************************************************** 
c  xの導関数
      real function fx(t,x,v) 
      real t,x,v
      real alpha
      common /para/alpha
      fx=v                                           
      return                                                            
      end                                                               
*********************************************************************** 
c  vの導関数
      real function fv(t,x,v)
      real t,x,v
      real alpha
      common /para/alpha
      fv=-x-alpha*v
      return                                                            
      end                                                               
*********************************************************************** 
c   1次元運動方程式を解くためのルンゲ・クッタ法のサブルーチン
      subroutine mtn1d(t,x,v,dt)                                       
      real t,x,v,dt
      real h1,kx(4),kv(4),kx1,kv1,x1,v1,x2,v2,x3,v3
c
      h1=dt/2.0                                                         
      kx(1)=fx(t,x,v)*dt                                              
      kv(1)=fv(t,x,v)*dt                                              
      x1=x+kx(1)/2.0
      v1=v+kv(1)/2.0
c
      kx(2)=fx(t+h1,x1,v1)*dt                                          
      kv(2)=fv(t+h1,x1,v1)*dt                                          
      x2=x+kx(2)/2.0
      v2=v+kv(2)/2.0
c
      kx(3)=fx(t+h1,x2,v2)*dt                                          
      kv(3)=fv(t+h1,x2,v2)*dt                                          
      x3=x+kx(3)
      v3=v+kv(3)
c
      kx(4)=fx(t+dt,x3,v3)*dt                                          
      kv(4)=fv(t+dt,x3,v3)*dt                                          
c
      t=t+dt                                                            
      kx1=(kx(1)+2.0*kx(2)+2.0*kx(3)+kx(4))/6.0              
      kv1=(kv(1)+2.0*kv(2)+2.0*kv(3)+kv(4))/6.0              
      x=x+kx1
      v=v+kv1
c
      return
      end