Program blackhole
      Implicit None
      Real*8 EPS_BH,EPS_TURN                                                                      //各変数、定数の設定
      Real*8 fugo,dr,xold,x,y
      Real*8 t,x0,x1,x2,h,b,t0,tf
      Parameter(t0=0.0D0,tf=10.0D0)
      EPS_BH =0.00001D0
      EPS_TURN=0.0000001D0
     
      fugo=-1.0D0
      b=8.5D0
      x0=10.0D0
      xold=x0
      x1=x0
      x2=dasin(b/x0)

      
      do t=t0,tf,h
         Do while((x1 .le. x0) .and. (x1 .ge. (1.0D0+EPS_BH)))    　　　　　　//x1がx0以下、ある数値以上の場合のみ、計算する
            h=0.1D0*sqrt((x1**3-(x1-1.0D0)*(b**2))/(x1**3))                        //ある数値以下になった場合、ブラックホールに落ちたとみなす
            Call rk4(t,x1,x2,h,fugo,b) 　　　　　　　　　　　　　　　　　　　　　　//下のサブルーチンを呼び出す
            dr=x1-xold
            if(abs(dr).lt.EPS_TURN) fugo=1.0D0
            xold=x1
            x=x1*cos(x2)
            y=x1*sin(x2)
            Write(*,*) x1,x2,fugo      
         end do
      end do
      End

      Subroutine rk4(t,x1,x2,h,fugo,b)
      Implicit None
      real*8 t,x1,x2,h,fugo,b
      real*8 k1,k2,k3,k4
      real*8 l1,l2,l3,l4
      real*8 diff_eq1,diff_eq2
      
      k1=diff_eq1(t,x1,x2,fugo,b)　　　　　　　　　　　　　　　　　　　//ルンゲ・クッタ法
      l1=diff_eq2(t,x1,x2,b)
      k2=diff_eq1(t+0.5D0*h,x1+0.5D0*h*k1,x2+0.5D0*h*l1,fugo,b)
      l2=diff_eq2(t+0.5D0*h,x1+0.5D0*h*k1,x2+0.5D0*h*l1,b)
      k3=diff_eq1(t+0.5D0*h,x1+0.5D0*h*k2,x2+0.5D0*h*l2,fugo,b)
      l3=diff_eq2(t+0.5D0*h,x1+0.5D0*h*k2,x2+0.5D0*h*l2,b)
      k4=diff_eq1(t+h,x1+h*k3,x2+h*l3,fugo,b)
      l4=diff_eq2(t+h,x1+h*k3,x2+h*l3,b)
      
      x1 = x1+(1/6.D0)*h*(k1+2.0D0*k2+2.0D0*k3+k4)
      x2 = x2+(1/6.D0)*h*(l1+2.0D0*l2+2.0D0*l3+l4)
      
      Return
      End
      
      Real*8 function diff_eq1(t,x1,x2,fugo,b)
      Implicit None
      real*8 t,x1,x2,fugo
      real*8 b
     
      diff_eq1=fugo*(1.0D0-(1.0D0/x1))*　　　　　　　　　　//微分方程式１
     &sqrt(1.0D0-(b**2/x1**2)+(b**2/x1**3))
     
      return
      end
      
      Real*8 Function diff_eq2(t,x1,x2,b)
      Implicit none
      real*8 t,x1,x2
      real*8 b
     
      diff_eq2=(1.0D0-(1.0D0/x1))*(b/(x1**2))　　　　　　　　　　　　　　//微分方程式２
     
      Return
      end