Damped oscillator model matlab code

% 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);

Download matlab code 2

% Program to plot the absorption response and phase shift of a bound electron as a function of Gamma, the

% damping-term parameter

clear; close all

prompt = 'Lower limit of range of photon energy in units of E_B? [0.2]';

Emin = input(prompt);

if isempty(Emin)

Emin = 0.2;

end

prompt = 'Upper limit of range of photon energy in units of E_B? [2]';

Emax = input(prompt);

if isempty(Emax)

Emax = 2;

end

prompt = 'Lower limit of range of Gamma in units of E_B? [0.05]';

Gmin = input(prompt);

if isempty(Gmin)

Gmin = 0.05;

end

prompt = 'Upper limit of range of Gamma in units of E_B? [0.35]';

Gmax = input(prompt);

if isempty(Gmax)

Gmax = 0.35;

end

set(0,'defaultfigurecolor',[1 1 1]); % Graph background white instead of default grey

figure('ToolBar','none');

vid = VideoWriter('dampedOscillatorGraph.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 30;

open(vid);

for gamma = Gmin:0.00125:Gmax

newplot

hold off

omega = Emin:0.0001:Emax;

omega0 = 1;

amp = ((omega0^2 - omega.^2).^2 + (gamma*omega).^2).^(-1);

delta = atan2(gamma*omega,(omega0^2 - omega.^2));

yyaxis left

plot(omega,amp,'b','LineWidth', 2.0)

set(gca,'YTickLabel',[]);

set(gca,'YTick',[])

ax = gca;

ax.YAxis(1).Color = 'k';

ylabel('Amplitude [arb. units]','color','b');

ylim([0,1.1*max(amp)])

hold on

yyaxis right

plot(omega,delta,'r','LineWidth', 2.0)

ax = gca;

ax.YAxis(2).Color = 'k';

ax.YLabel(1).Color = 'r';

ylabel('Phase [radians]','color','r');

set(gca,'ytick',[0:pi/2:pi]) % where to set the tick marks

set(gca,'yticklabels',{'0','\pi/2','\pi'}) % give them user-defined labels

ylim([0,3.3])

set(gca,'fontsize',18);

set(gca,'TickLength',[0.02, 1]);

set(gca,'linewidth',2);

xlabel('h\nu/E_B');

str1 = 'h\Gamma/E_B = ';

str2 = num2str(gamma,'% 4.3f');

strTot = [str1,str2];

hText = text(0.7*Emax, 2, strTot, 'FontSize',18);

frame = getframe(gcf);

writeVideo(vid,frame);

hold off

end

close(vid)