function meshBlackHole(N,M,rMinOverRSch,rMaxOverRSch,numRocketPoints,showTimes)
% usage: meshBlackHole(N,M,rMinOverRSch,rMaxOverRSch)
% N = number of mesh points
% M = Mass of the black hole in G = c = 1 units
% rMinOverRSch = The radius of the cylinder inside the back hole in uses of rSch: rMinOverRSch <= 1
% rMaxOverRSch = Maximum radius of the plot in unit of rSch: rMaxOverRSch > 1
% Where rSch = 2.0*G*M/c^2 is the Schwarzschild Radius
% numRocketPoints = Number of points to use to draw the rocket. If less than three, the rocket is not drawn. If more than 3 the rocket is drawn using 60 points.
% showTimes = 1 to display times on the title line.

G = 1.0;
c = 1.0;

U = 1.0; % Kruskal-Zekeres coordinate to match parabola to cylinder.of radius rMin.

rSch = 2.0*G*M/c^2; % The Schwarzschild Radius

rMin = rMinOverRSch*rSch;
rMax = rMaxOverRSch*rSch;

% The parabola matches to the cylinder at Schwarzschild coordiantes rMin, tMax.
% Find tMax from the value of U above:
temp = sqrt(1.0 - rMin/rSch)*exp(rMin/(2.0*rSch));
tMax = (2.0*rSch/c)*asinh(U/temp);

% Find the height of the cylinder, zMax:
zMax = sqrt(rSch/rMin - 1.0)*c*tMax;

K = zMax - sqrt(4*rSch*rMin);

% set up the x,y grid:
delta = rMax/N;
xVec = -1.0*rMax + (0:(2*N-1))*delta;
yVec = -1.0*rMax + (0:(2*N-1))*delta;

for i = 1:2*N
    for j = 1:2*N
        x = xVec(j); % x value for column j
	y = yVec(i); % y value for row i
	r = sqrt(x*x + y*y);
	phi = atan2(y,x);
	if r >= rMin
	  z(i,j) = sqrt(4.0*rSch*r) + K;
	else
          z(i,j) = 0.0;
	end
    end
end
surf(xVec,yVec,z,1.0*z)
hold
surf(xVec,yVec,-1.0*z,1.0*z)
shading interp
axis off
axis([-1.0*rMax,1.0*rMax,-1.0*rMax,1.0*rMax,-40,40])
view([1.0*rMax 0.6*rMax 1.1*rSch]);
%hidden off
hidden on
%colormap(pink)
colormap(copper)
brighten(-0.3)
alpha(0.9)
%colormap(gray)
%colormap(bone)
%colormap(0.5*bone)
%colormap(hot)
%colormap(jet)
%colormap(winter)
%colormap(cool)
%colormap(hsv)
%colormap(0.3*hot)

% Plot the Schwarzschild radius:

phi0 = (0:N-1)*2.0*pi/N;
xCirc = rSch*cos(phi0);
yCirc = rSch*sin(phi0);
zCirc = ones(N)*(sqrt(4.0*rSch*rSch) + K);
plot3(xCirc,yCirc,zCirc,'r');

plot3(xCirc,yCirc,-1.0*zCirc,'r');

%return;

% Plot the rocket if numRocketPoints >=3

if numRocketPoints > 0
 % Set up to compute times for fall into supermassive BH at center of Milkyway.
 rSchSun = 2.95; % Schwarzschild radius of the Sun in kilometers
 rSchMWSBH = 4.1e6*rSchSun; % Schwarzschild radius of supermassive BH at center of Milkyway.
 cKmPerSec = 2.99e5; % Speed of light in km/s

 rocketInitialRFactor = 0.9; % Used to give the inital position of the rocket.
 rocketScaleSize = 0.5; % The rocket initially extends from r0 to (1+rocketScaleSize)*r0.
 r0 = rocketInitialRFactor*rSch; % Starting r coordinate of the nose of the rocket.

 % Find the initial E.F. time coordinate tPrime0, when r = r0 the initial position of the nose of the rocket.
 % The initial U coordinate of the rocket is U0Rocket.
 U0Rocket = 0.8;
 rMin0 = r0;
 temp = sqrt(1.0 - rMin0/rSch)*exp(rMin0/(2.0*rSch));
 tSch0 = (2.0*rSch/c)*asinh(U0Rocket/temp);
 tPrime0 = tSch0 + (rSch/c)*log(1.0 - rMin0/rSch); % This is the E.F. time the rocket was at r0.

 % Find the initial E.F. time coordinate for the current slice, where the parabola meets the cylinder.
 tPrimeMax = tMax + (rSch/c)*log(1.0 - rMin/rSch);

 % Find the coordinate time inside the supermassive BH.
 coodTimeInsideBH0 = rocketInitialRFactor*rSchMWSBH/cKmPerSec;
 coodTimeInsideBH = coodTimeInsideBH0 - rMinOverRSch*rSchMWSBH/cKmPerSec;
 rMinMWSBH0 = rocketInitialRFactor*rSchMWSBH;
 % Find the proper time at the center of the supermassive BH (along Schwarzschild = constant line, e.g., z = 0 in the embedding diagram).
 rMinMWSBH = rMinOverRSch*rSchMWSBH;
 rRMWSBH = rMinMWSBH + (0:99)*(rMinMWSBH0 - rMinMWSBH)/99.0;
 integrandTimeInsideBH = ones(1,100)./sqrt(rSchMWSBH*ones(1,100)./rRMWSBH - 1.0);
 tInsideBH = trapz(rRMWSBH,integrandTimeInsideBH)/cKmPerSec;

end

if numRocketPoints > 3
 % The rocket is drawn using 60 points.

 % Starting the rocket from these initial conditions,
 for k=1:20
  deltaPhi = ((k-1)/19.0)*.2;
  r = (1.0 + (rocketScaleSize/0.2)*deltaPhi)*r0;
  tPrime = tPrime0;
  phi0 = pi + deltaPhi;
  % Next, integrate, using steps of deltaTau;
  deltaTau = 0.01;
  Tau = 0;
  % the radial geodesic equations until we
  % either get to the current E.F. time or the current rMin.
  while (tPrime < tPrimeMax) & (r > rMin)
   r = r - sqrt(rSch/r)*c*deltaTau;
   tPrime = tPrime + deltaTau*(1.0 - (rSch/r)^(3.0/2.0))/(1.0 - rSch/r);
   Tau = Tau + deltaTau;
  end
  if k == 1
     rocketProperTime = (2.0/3.0)*(rSchMWSBH/cKmPerSec)*( (r0/rSch)^1.5 - (r/rSch)^1.5 ); % save the proper time that has elapsed on the rocket
  end
  if r >= rMin
    xLine(k) = r*cos(phi0);
    yLine(k) = r*sin(phi0);
    zLine(k) = sqrt(4.0*rSch*r) + K;
  else
    tSch = tPrime - (rSch/c)*log(1.0 - r/rSch);
    xLine(k) = rMin*cos(phi0);
    yLine(k) = rMin*sin(phi0);
    zLine(k) = sqrt(rSch/r - 1.0)*c*tSch;
  end
 end

 deltaPhiMax = deltaPhi;
 for k=21:40
  deltaPhi = ((k-20-1)/19.0)*.4;
  r = (1.0 + (rocketScaleSize/0.2)*deltaPhiMax)*r0;
  tPrime = tPrime0;
  phi0 = pi + deltaPhiMax - deltaPhi;
  % Next, integrate, using steps of deltaTau;
  deltaTau = 0.01;
  Tau = 0;
  % the radial geodesic equations until we
  % either get to the current E.F. time or the current rMin.
  while (tPrime < tPrimeMax) & (r > rMin)
   r = r - sqrt(rSch/r)*c*deltaTau;
   tPrime = tPrime + deltaTau*(1.0 - (rSch/r)^(3.0/2.0))/(1.0 - rSch/r);
   Tau = Tau + deltaTau;
  end
  if r >= rMin
    xLine(k) = r*cos(phi0);
    yLine(k) = r*sin(phi0);
    zLine(k) = sqrt(4.0*rSch*r) + K;
  else
    tSch = tPrime - (rSch/c)*log(1.0 - r/rSch);
    xLine(k) = rMin*cos(phi0);
    yLine(k) = rMin*sin(phi0);
    zLine(k) = sqrt(rSch/r - 1.0)*c*tSch;
  end
 end

 % Here r is the r coodrinate of the back of the rocket, and tPrimeMax is the E.F. time of this slice.
 % If this is outside the horizon, find the Schwarschild time for fall into super massive BH as center of the galaxy.
 if r > rSch

    % Find the Schwarzschild time inside the BH where the parabola and cylinder join for the super massive BH.
    % Convert this to the E.F. time of this slice for the super massive BH.
    temp = sqrt(1.0 - rMin/rSch)*exp(rMin/(2.0*rSch));
    tSchMWSBM = (2.0*rSchMWSBH/cKmPerSec)*asinh(U0Rocket/temp);
    tPrimeMWSBM = tSchMWSBM + (rSchMWSBH/cKmPerSec)*log(1.0 - rMin/rSch);

    % Now find the Schwarzschild time for the back of the rocket for this same E.F. time (since the parabola is a constant E.F. slice).
    tOutSideBH = tPrimeMWSBM - (rSchMWSBH/cKmPerSec)*log(r/rSch - 1.0);
 else
     tOutSideBH = -1; % time is infinite.
 end

 for k=41:60
  deltaPhi = ((k-40-1)/19.0)*.2;
  r = (1.0 + (rocketScaleSize/0.2)*(deltaPhiMax - deltaPhi))*r0;
  tPrime = tPrime0;
  phi0 = pi - deltaPhiMax + deltaPhi;
  % Next, integrate, using steps of deltaTau;
  deltaTau = 0.01;
  Tau = 0;
  % the radial geodesic equations until we
  % either get to the current E.F. time or the current rMin.
  while (tPrime < tPrimeMax) & (r > rMin)
   r = r - sqrt(rSch/r)*c*deltaTau;
   tPrime = tPrime + deltaTau*(1.0 - (rSch/r)^(3.0/2.0))/(1.0 - rSch/r);
   Tau = Tau + deltaTau;
  end
  if r >= rMin
    xLine(k) = r*cos(phi0);
    yLine(k) = r*sin(phi0);
    zLine(k) = sqrt(4.0*rSch*r) + K;
  else
    tSch = tPrime - (rSch/c)*log(1.0 - r/rSch);
    xLine(k) = rMin*cos(phi0);
    yLine(k) = rMin*sin(phi0);
    zLine(k) = sqrt(rSch/r - 1.0)*c*tSch;
  end
 end

 fill3(xLine,yLine,zLine,'k');

 if showTimes == 1
  if tOutSideBH > -1
    timeString = sprintf('TIME INSIDE = %4.1f sec.     ROCKET TIME = %4.1f sec.    TIME OUTSIDE = %5.1f sec.',tInsideBH,rocketProperTime,tOutSideBH);
  else
    timeString = sprintf('TIME INSIDE = %4.1f  sec.    ROCKET TIME = %4.1f  sec.   TIME OUTSIDE = INFINITY',tInsideBH,rocketProperTime);
  end
  title(timeString);
 end

 return;

elseif numRocketPoints == 3

 % Starting the rocket from these initial conditions,
 k = 1;
 r = r0;
 tPrime = tPrime0;
 phi0 = pi;
 % Next, integrate, using steps of deltaTau;
 deltaTau = 0.01;
 Tau = 0;
 % the radial geodesic equations until we
 % either get to the current E.F. time or the current rMin.
 while (tPrime < tPrimeMax) & (r > rMin)
  r = r - sqrt(rSch/r)*c*deltaTau;
  tPrime = tPrime + deltaTau*(1.0 - (rSch/r)^(3.0/2.0))/(1.0 - rSch/r);
  Tau = Tau + deltaTau;
 end
 if k == 1
     rocketProperTime = (2.0/3.0)*(rSchMWSBH/cKmPerSec)*( (r0/rSch)^1.5 - (r/rSch)^1.5 ); % save the proper time that has elapsed on the rocket
 end
 if r > rMin
   xLine(k) = r*cos(phi0);
   yLine(k) = r*sin(phi0);
   zLine(k) = sqrt(4.0*rSch*r) + K;
 else
   tSch = tPrime - (rSch/c)*log(1.0 - r/rSch);
   xLine(k) = rMin*cos(phi0);
   yLine(k) = rMin*sin(phi0);
   zLine(k) = sqrt(rSch/r - 1.0)*c*tSch;
 end

 % Repeat with:
 k = 2;
 r = (1.0+rocketScaleSize)*r0;
 tPrime = tPrime0;
 phi0 = pi + 0.2;
 % Next, integrate, using steps of deltaTau;
 if rMinOverRSch == 0.81332555
   % The rMinOverRSch == 0.81332555 and deltaTau = 0.00000001 are special
   % values to get within 1.e-10 of the horizon to compute tOutSideBH there.
   deltaTau = 0.00000001;
 else
   deltaTau = 0.01;
 end
 Tau = 0;
 % the radial geodesic equations until we
 % either get to the current E.F. time or the current rMin.
 while (tPrime < tPrimeMax) & (r > rMin)
  r = r - sqrt(rSch/r)*c*deltaTau;
  tPrime = tPrime + deltaTau*(1.0 - (rSch/r)^(3.0/2.0))/(1.0 - rSch/r);
  Tau = Tau + deltaTau;
 end
 if r >= rMin
   xLine(k) = r*cos(phi0);
   yLine(k) = r*sin(phi0);
   zLine(k) = sqrt(4.0*rSch*r) + K;
 else
   tSch = tPrime - (rSch/c)*log(1.0 - r/rSch);
   xLine(k) = rMin*cos(phi0);
   yLine(k) = rMin*sin(phi0);
   zLine(k) = sqrt(rSch/r - 1.0)*c*tSch;
 end

 % If the back of the rocket is still outside the BH, find the Schwarzschild time for fall into super massive BH as center of the galaxy.
 % Here r is the r coodrinate of the back of the rocket, and tPrimeMax is the E.F. time of this slice.
 % If this is outside the horizon, find the Schwarschild time for fall into super massive BH as center of the galaxy.
 if r > rSch

    % Find the Schwarzschild time inside the BH where the parabola and cylinder join for the super massive BH.
    % Convert this to the E.F. time of this slice for the super massive BH.
    temp = sqrt(1.0 - rMin/rSch)*exp(rMin/(2.0*rSch));
    tSchMWSBM = (2.0*rSchMWSBH/cKmPerSec)*asinh(U0Rocket/temp);
    tPrimeMWSBM = tSchMWSBM + (rSchMWSBH/cKmPerSec)*log(1.0 - rMin/rSch);

    % Now find the Schwarzschild time for the back of the rocket for this same E.F. time (since the parabola is a constant E.F. slice).
    tOutSideBH = tPrimeMWSBM - (rSchMWSBH/cKmPerSec)*log(r/rSch - 1.0);
 else
     tOutSideBH = -1; % time is infinite.
 end

 % Repeat with:
 k = 3;
 r = (1.0+rocketScaleSize)*r0;
 tPrime = tPrime0;
 phi0 = pi - 0.2;
 % Next, integrate, using steps of deltaTau;
 deltaTau = 0.01;
 Tau = 0;
 % the radial geodesic equations until we
 % either get to the current E.F. time or the current rMin.
 while (tPrime < tPrimeMax) & (r > rMin)
  r = r - sqrt(rSch/r)*c*deltaTau;
  tPrime = tPrime + deltaTau*(1.0 - (rSch/r)^(3.0/2.0))/(1.0 - rSch/r);
  Tau = Tau + deltaTau;
 end
 if r >= rMin
   xLine(k) = r*cos(phi0);
   yLine(k) = r*sin(phi0);
   zLine(k) = sqrt(4.0*rSch*r) + K;
 else
   tSch = tPrime - (rSch/c)*log(1.0 - r/rSch);
   xLine(k) = rMin*cos(phi0);
   yLine(k) = rMin*sin(phi0);
   zLine(k) = sqrt(rSch/r - 1.0)*c*tSch;
 end

 fill3(xLine,yLine,zLine,'k');

 if showTimes == 1
  if tOutSideBH > -1
    timeString = sprintf('TIME INSIDE = %4.1f sec.     ROCKET TIME = %4.1f sec.    TIME OUTSIDE = %5.1f sec.',tInsideBH,rocketProperTime,tOutSideBH);
  else
    timeString = sprintf('TIME INSIDE = %4.1f  sec.    ROCKET TIME = %4.1f  sec.   TIME OUTSIDE = INFINITY',tInsideBH,rocketProperTime);
  end
  title(timeString);
 end

 return;
end
