format long;
fflush(stdout);
hold off;

M0 =  input(" Please enter the mass of the black-hole in Kilograms  -->  " );

if (M0<= 0)
sprintf(" Invalid data. The program is exiting.  ")
exit
end

sprintf(" Schwarzschild radius of your black-hole is  =  ")

Rs = (2*6.6732*(10^(-11))*M0)/(3*10^8)^2

r00 = input(" Please enter initial radial coordinate in meters. Your particle's starting position should be outside the Schwarzchild radius. -->  ");

if (r00<Rs)
sprintf (" Particle is inside the Schwarzschid radius. This program is not valid for such initial conditions. Hence exiting. ")
exit
end

if (r00==Rs)
sprintf (" Particle is on the Schwarzschid radius. This program is not valid for such initial conditions. Hence exiting. ")
exit
end


p0 = input(" Please enter initial azimuthal coordinate (a pure number) -->  ");

x = 3*10^8;
r0 = r00/(Rs);
y = (1-(1/r0));


tr00 = input(" Please enter initial radial velocity in meters/sec -->  ");
tr0 = tr00/(x);

validangvel = (x/Rs)*(1/r0)*sqrt(y - (((tr0)^2)/(y)))

tp00 = input(" Please enter initial azimuthal coordinate veocity in sec^-1. { A valid azimuthal velocity should lie in the interval (-validangvel,validangvel).} -->  ");

if (tp00 >= validangvel | tp00 <= -validangvel)

sprintf(" Positivity of energy (hence also the requirement that the test particle's velocity should be less than the speed of light) is not satisfied by the space-time position and coordinate velocity of the test-particle. Hence exiting." )
exit
end
  

tp0 = (tp00*Rs)/(x);

f = sqrt(((tr0)^2/(y)) + (r0*tp0)^2)

sprintf(" Your test particle is moving in the Schwarzchild space-time of the given black-hole at a velocity of 'f' times the speed of light" )


if (tp00 != 0)
n = input(" Please enter what is the maximum number of revolutions of the black-hole you want to see  -->  ");
end
 
l = input ("Please enter how many times the Schwarzschild radius due you want to track the particle  -->  ");

if (l<0)
sprintf(" Invalid Data. Hence exiting")
exit
end

sprintf("The particle's motion will be tracked till just before it reaches the Event Horizon (if it does)")       

function x = B(y)
x = 1-(1/(y));
end

B0 = B(r0);
global j = (((r0)^2)/B0)*(tp0);
global er = (1/B0)-((1/((B0)^3))*((tr0)^2))-((j/(r0))^2);


function x = rhsr(y)
global j;
global er;
z = B(y);
x = (1/z)-(er)-(j^2/y^2);
end

function x = rhsp(y)
global j;
z = B(y);
x = (j*z)/(y^2);
end

function x = rhs(y)
z = B(y);
a = rhsr(y);
b = rhsp(y);

if ((a>=0) & (b!=0))
x = ((z^(3/2))*sqrt(a))/(b);
else
x = "Instability";
end

end



rr = 1;
pp = [0:0.01:2*pi];
x = rr*cos(pp);
y = rr*sin(pp);
plot (x,y,"*r");
hold on;
r = r0;
p = p0;


h = 0.01;


if (tp0!=0)

while ((abs(p-p0) <= 2*n*pi) & (1<r) & (r<l*r0)  )

z1 = rhs(r);

if (z1 != "Instability")
k1 = h * z1;
else
pause;
end

z2 = rhs(r+(k1/2));

if (z2 != "Instability")
k2 = h * z2;
else
pause;
end

z3 = rhs(r+(k2/2));


if (z3 != "Instability")
k3 = h * z3;
else
pause;
end

z4 = rhs (r+k3);

if (z4 != "Instability")
k4 = h * z4;
else
pause;
end

r = r + (k1/6) + (k2/3) + (k3/3) + (k4/6);

plot((r)*cos(p),(r)*sin(p),"*r");
hold on;
p+=h;
end


else



while (1<r & r<l*r0 )
z = B(r);

z1 = rhsr(r);

if (z1 > 0)
k1 = h *sqrt(z1)*(z^(3/2));
else
pause;
end

z2 = rhsr(r+(k1/2));

if (z2 >0 )
k2 = h *sqrt(z2)*(z^(3/2));
else
pause;
end

z3 = rhsr(r+(k2/2));

if (z3 > 0 )
k3 = h *sqrt(z3)*(z^(3/2));
else
pause;
end

z4 = rhsr(r+k3);

if (z4 > 0)
k4 = h *sqrt(z4)*(z^(3/2));
else
pause;
end

r = r + (k1/6) + (k2/3) + (k3/3) + (k4/6);

plot((r)*cos(p),(r)*sin(p),"*r");
hold on;
end

end