c  ルンゲ・クッタ法による1次元運動方程式の解法
c      
      real*8 t,x,v,dt
      real*8 t0,x0,v0
      real*8 tlim
      real*8 a,b
      real*8 omega,f0
      real*8 pi
      integer n,nwrite, i
c  t：時間、x：位置、v：速度、dt：時間の刻み幅
c  t0、x0、v0：初期値
c  tlim：時間の最大値
c  omega：強制力の振動数、f0：強制力の振幅
c  alpha：速度に比例する抵抗の係数
c
      common /para/a,b,omega,f0
c
c  出力ファイルの指定
      open(unit=2,file='osc.dat')
      open(unit=3,file='poincare.dat')
      pi=3.1415926535897932385
c  初期値
      t0=0.0                                                            
      x0=0.0
      v0=1.0
c
      a=0.05
      b=1.0
      omega=1.0
      dt=2*pi/1000
c  強制力の振動数、強制力の振幅
      write(6,*) ' f0 ?'
      read(5,*) f0
c  時間の刻み幅
c
      write(6,*) ' tlim/(2*pi) ?'
      read(5,*) tlim      
      tlim=2*pi*tlim
      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,1000).eq.0) then
        write(3,*) x,v
      endif
      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*8 function fx(t,x,v) 
      real*8 t,x,v
      real*8 a,b,omega,f0
      common /para/a,b,omega,f0
      fx=v                                           
      return                                                            
      end                                                               
*********************************************************************** 
c  vの導関数
      real*8 function fv(t,x,v)
      real*8 t,x,v
      real*8 a,b,omega,f0
      common /para/a,b,omega,f0
      fv=-a*v-b*x**3+f0*cos(omega*t)
      return                                                            
      end                                                               
*********************************************************************** 
c   1次元運動方程式を解くためのルンゲ・クッタ法のサブルーチン
      subroutine mtn1d(t,x,v,dt)                                       
      real*8 t,x,v,dt
      real*8 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