% Damped oscillator model. Increase the photon energy linearly from 0.25 of
% E_B to twice E_B. For energies well below E_B, the electron moves exactly
% in the opposite direction to the E-field (F = -eE). The reradiated field
% is exactly out of phase with the electron motion, hence is in phase with
% the incident radiation. At photon energies well above E_B, the electron
% cannot respond instantaneously and assumes a motion in phase with the
% E-field amplitude. The reradiated field is therefore exactly out of phase
% with the incident radiation and the bound electron therefore dissipates
% energy (i.e., absorption of the incident field occurs). 

clear
close all
pos1 = [0 0 1 1];

% Parameters
E_B = 2.0;  % Binding energy in arbitrary units
dt = 0.02;  % Time step
nx = 200;   % Number of points for wave curves
gamma = 0.1;  % Damping parameter

% Create spatial grid for waves
x = linspace(-10, 10, nx);

% Time vector
t = 0:dt:15;

% Initialize figure
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gcf,'Color','w'); % Black background

% Create video writer object
v = VideoWriter('dampedOscillator.mp4','MPEG-4');
v.FrameRate = 30;
open(v);

% Scattered wave angles 
theta = 0;  % Angle in xy plane
phi = pi/2;    % Angle from xy plane

subplot('position',pos1);

% Calculate scattered wave path
[Y_scatter, X_scatter, Z_scatter] = deal(x*cos(theta), ...
    x*sin(theta)*cos(phi), ...
    x*sin(theta)*sin(phi));

% Calculate displacement direction vector for scattered wave
displacement_dir = [-sin(theta), cos(theta)*cos(phi), cos(theta)*sin(phi)];

for i = 1:length(t)
    % Calculate photon energy (gradually increasing)
    E_photon = E_B * (0.25 + 1.75 * (i/length(t)));

    % Wave parameters
    k = E_photon;  % Wave number proportional to energy
    omega = E_photon;  % Angular frequency

    % Calculate electron response using damped oscillator model
    % Phase relative to driving field follows arctan(γω/(ω₀² - ω²))
    phase = atan2(gamma * omega, E_B^2 - omega^2);

    % Amplitude from damped oscillator model
    amplitude = 1.0 / sqrt((E_B^2 - omega^2)^2 + (gamma * omega)^2);

    % Electron position (oscillating in y-direction)
    electron_pos = [0, -amplitude * cos(omega*t(i) - phase), 0];

    % Calculate waves
    % Incident wave (along x-axis, oscillating in z)
    incident_wave = cos(k*x - omega*t(i));
    Z_incident = -incident_wave;

    % Scattered wave always π out of phase with electron motion
    r = sqrt(X_scatter.^2 + Y_scatter.^2 + Z_scatter.^2);
    scattered_wave = amplitude*cos(k*r - omega*t(i) - phase - pi)./sqrt(r + 0.1);

    % Add oscillation to scattered wave perpendicular to propagation direction
    scattered_wave_displacement_x = scattered_wave .* displacement_dir(1);
    scattered_wave_displacement_y = scattered_wave .* displacement_dir(2);
    scattered_wave_displacement_z = scattered_wave .* displacement_dir(3);

    clf;
    hold on;

    % Plot incident wave
    plot3(x, zeros(size(x)), Z_incident, 'b-', 'LineWidth', 2);

    % Plot scattered wave
    plot3(X_scatter + scattered_wave_displacement_x, ...
        Y_scatter + scattered_wave_displacement_y, ...
        Z_scatter + scattered_wave_displacement_z, 'r-', 'LineWidth', 2);

    % Plot electron
    plot3(electron_pos(1), electron_pos(3), -electron_pos(2), 'ko', ...
        'MarkerFaceColor', 'yellow', 'MarkerSize', 15);

    % Add coordinate system for reference
    plot3([0 2], [0 0], [0 0], 'k-', 'LineWidth', 1);  % x-axis
    plot3([0 0], [0 2], [0 0], 'k-', 'LineWidth', 1);  % y-axis
    plot3([0 0], [0 0], [0 2], 'k-', 'LineWidth', 1);  % z-axis
    text(2.2, 0, 0, 'z');
    text(0, 2.2, 0, 'x');
    text(0, 0, 2.2, 'y');
    set(gca, 'Projection','perspective');

    % Add legend and text
    legend('Incident Wave', 'Scattered Wave', 'Electron', 'Location',... 
        'northeast','FontSize', 20);
    text(-9, 4, 4, sprintf('E_{photon} = %.2f\nE_B = %.2f',...
        E_photon, E_B), 'FontSize', 20);

    % Set view properties
    view(45, 30);  % 3D view
    grid on;
    xlabel('x');
    ylabel('y');
    zlabel('z');
    axis([-10 10 -10 10 -5 5]);
    % Lock them in place
    axis manual
    % axis equal
    axis off
    % Capture frame
    frame = getframe(gcf);
    writeVideo(v, frame);
end

% Close video writer
close(v);
