program blackhole
real b, ri, rf, m
real dx !!++++++++++++++++++++++++ 적분 미소변위
real pi
real integral(1000,100)
integer i,j
real sm !+++++++++++++++ 더하는것
real theta
b=0.01 !++++++++++++ impact parameter
ri=50 !===================초기반지름
rf=25 !=======================최종반지름
m=10 !=================== 블랙홀질량
sm=0!++++++++++++++++ 합산 초기화
pi=30
do j=1,10,1
ri=ri+(j-1)*(-1)
do i=1,1000,1
dx=(ri-rf)/1000.0
c=(ri-dx*(i-1))
sm=w(c,b,m)*dx+sm
integral(i,j)=sm
if(ri-dx*(i-1)<2*m+0.1) then
c=0
end if
end do
end do
open(unit=3,file='orbit.txt',status='replace')
do j=1,10,1
do i=1,1000,1
write(3,*) (ri-dx*(i-1))*cos(integral(i,j)+pi),(ri-dx*(i-1))*sin(integral(i,j)+pi)
end do
end do
close(3)
end program blackhole
real function w(r,b,m)
real r
real b
real m
w=(1/(b**2)-1/(r*r)*(1-2*m/r))**(0.5)/(r*r)
return
end
