My plot for a photon traveling around a black hole isn't updating.

I am currently working on creating a program that plots the trajectory of a photon travelling around a 2-D black hole but am having troubles. My figure isn't updating with the trajectory of the photon. I have the program set up where the user can specify what the mass of the black hole is, where the initial position of the photon is, the direction of travel and the number of photons. I am still somewhat new to Matlab, so I understand that the program is far from optimized and is rudimentary, however, any help would be beneficial.

Here is the program:

%% Preclose everything
clear
clc
close all
%% Constants Needed
c = 299792458; % Speed of light in m/s
G = 6.67e-11; % Gravitational Constant
M_sun = 1.989e30; % Mass of the sun in kilograms
M_solar_mass = input('Enter the mass of the black hole in solar masses: ');
M_blackhole = M_solar_mass * M_sun; % Mass of the black hole in kilograms
rs = (2*G*M_blackhole)/(c*c); % Schwarzschild Radius in meters
rs_km = rs/1000; %Schwarzschild Radius in kilometers

%% Code to Generate Black Hole
theta = linspace(0, 2*pi, 100);
x_bh = rs_km*cos(theta);
y_bh = rs_km*sin(theta);
fill(x_bh, y_bh, 'k');
xlabel('x (km)');
ylabel('y (km)');
title('Photon Trajectory');
hold on;
%% Code to generate Accretion Disk
x_accretiondisk = 3*rs_km*cos(theta);
y_accretiondisk = 3*rs_km*sin(theta);
plot(x_accretiondisk, y_accretiondisk, "r", 'LineWidth', 2)

%% Code to generate Photon's "Unstable Orbit"
x_photonorbit = 1.5*rs_km*cos(theta);
y_photonorbit = 1.5*rs_km*sin(theta);
plot(x_photonorbit, y_photonorbit, "g", 'LineWidth', 2)
hold on;
%% Code to Generate an array of photons approaching
% Get initial position of "centered" vector
InitialPosition = input('Enter the initial position [x0 y0] in km: ');
% Ensure that "centered" photon is not inside Schwarzschild Radius
center = [0 0];
dist1 = norm(InitialPosition - center);
if dist1 < rs_km
    error('Photons will start inside the black hole. Please select a different initial position.');
end

%Rather than entering a velocity vector, the photon will travel at the
%speed of light and the user will enter the angle at which it will travel
angle_degrees = input('Enter the angle the photons will travel in degrees: ');
angle_radians = deg2rad(angle_degrees);
velocity_x = -sin(angle_radians);
velocity_y = cos(angle_radians);
InitialVelocity = c * [velocity_x, velocity_y];

% Distance between photon's in km
d = rs_km / 10;

%Calculate distance traveled in one time step
dist2 = norm(InitialVelocity);

%Normalilze Velocity Vector
InitialVelocity_norm = InitialVelocity/dist2;

%Create an array of photons
num_photons = input('How many photons would you like to send? (Note - the more photons, the longer the program will take to run): '); % Can be changed to increase/decrease time
photon_pos = zeros(num_photons, 2);
for n = 1:num_photons
    photon_pos(n,1:2) = InitialPosition + n*d*InitialVelocity_norm;
end
photon_vel = repmat(InitialVelocity_norm, num_photons, 1);
plot(photon_pos(:,1), photon_pos(:,2), '.', 'LineWidth', 2);

%% Plot the trajectory of the photons
% Calculate the radius of the "Innermost Stable Circular Orbit (ISCO)
r_isco = 6*rs_km;
% Velocity of the Photon at the ISCO in km/s
v_isco = sqrt(G*M_blackhole/r_isco)/1000;
%The orbital period of a photon at the ISCO (needed to determine how long
%to run the program)
T_isco = 2*pi*r_isco/v_isco;

% Define the time step and total time
t_final = 1.5*T_isco;
dt = t_final/100;

% Initialize the plot
h = figure;
hold on
h_traj = plot(nan, nan);
h_time = text(0.1, 0.9, sprintf('Time: %.2f s', 0), 'Units', 'normalized');

for n = 1:num_photons
    pos = photon_pos(n,:);
    vel = photon_vel(n,:);

    tspan = [0 t_final];

    options = odeset('InitialStep', dt);
    [t, y] = ode45(@(t,y) geodesic(t, y, M_blackhole), tspan, [pos, vel], options);

    photon_traj{n} = y(:,1:2);
    set(h_traj, 'XData', photon_traj{n}(:,1), 'YData', photon_traj{n}(:,2))
    
    % Update the time text and refresh the figure
    set(h_time, 'String', sprintf('Time: %.2f s', t(end)))
    drawnow
    
    % Pause briefly to slow down the animation
    pause(0.01)
end


% NOTE: First time running the program, it crashed due to memory allocation
% from the ODE45 Runge Kutta Method. A more streamline and rudimentary
% approach is needed, perhaps the Leapfrog algorithm or the Verlet
% algorithm. Further trials are needed with various ODE's 

Here is the function "geodesic.m" that I am also using:

function dy = geodesic(t, y, M_blackhole)
    x = y(1);
    yy = y(2);
    vx = y(3);
    vy = y(4);

    c = 299792458;
    G = 6.67e-11;
    rs = (2*G*M_blackhole)/(c^2);
    alpha = 1-rs/x;
    beta = -rs/(2*x);
    
    dy  = zeros(4,1);
    dy(1) = vx;
    dy(2) = vy;
    dy(3) = alpha*beta*(vy^2 - vx^2) - (G*M_blackhole/x^2) * vx;
    dy(4) = alpha*beta*(vx^2 - vy^2) - (G*M_blackhole/x^2) * vy;
end

[Torsten]

Instead of the input lines, could you give values to the parameters expected there as user inputs ?

And could you mark the line in your code where you expect a plot to be updated ? Because I only see a "drawnow" command in your loop, but no corresponding "plot" command.

Answer 1

[Constantino Carlos Reyes-Aldasoro]

Is h_time being updated? if yes, then it may be due to the initial display with h_traj = plot(nan, nan);

try not using nans but a valid plot and that may set the correct parameters and then the plot may be updated each cycle of the for.