clear all
clc
% Define the pde function
function ut = pde(t, u)
% Define problem parameters
global x1 xu n ncall gamma_k;
global coeff1 coeff2 coeff3 coeff4 TT1 TT2 omega1 omega2;
global dx1 dx2 g1 g2 v gB tau_s trise;
x1 = 0; % Lower boundary of spatial domain
xu = 15; % Upper boundary of spatial domain
n = 190; % Number of discretization points
gB = 5*10^(-11); %Brilluoin gain factor
v = 0.2e9;
tau_s = 10e-9;
z = linspace(x1,xu,n);
gamma_k = 1/tau_s ; % Some parameter
coeff1 = 2; % Coefficient 1 for outer region
coeff2 = 1; % Coefficient 2 for inner region
coeff3 = 2; % Coefficient 3 for outer region
coeff4 = 1; % Coefficient 4 for inner region
omega1 = 0; % Omega1 value for outer region
omega2 = 2*pi*12.8e9 ; % Omega2 value for inner region
dx1 = 2; % Inner region starting point
dx2 = 12; % Inner region ending point
g1 = v*gamma_k*omega1*gB; % g1 value for outer region
g2 = v*gamma_k*omega2*gB; % g2 value for inner region
TT1 = (omega2)^2 - (omega1);
TT2 = (omega2)^2 -(omega2);
trise = 0.1e-9;
Ps = 10e-3;
Pb = 5e-3;
Aeff = 50e-6;
Es = (Ps/Aeff)^(1/2);
Eb = (Pb/Aeff)^(1/2);
tstart = 0;
tend = 90e-9;
t = [tstart ,tend];
% Calculate the time-dependent shape of the pulse waveform A(t)
t1 = (t - tau_s/2) / (trise/0.45);
t2 = (t + tau_s/2) / (trise/0.45);
A_t = sqrt((tanh(t1) - tanh(t2))/2);
% Calculate Es at z = 0 using the equation
Es_z0 = (Es - Eb) * A_t + Eb;
% Define the pulse function
pulse = @(t) A_t;
% Spatial index
for ii = 1:n
% Set initial conditions
if ii == 1
y1(1) = pulse(t);
y2(1) = 5e-3;
y3(1) = 1e-11;
y4(1) = 0;
y5(1) = 0;
y6(1) = 0;
end
if ii == n
y1(n) = pulse(tend);
y2(n) = 5e-3;
end
u0 =[y1(1); y2(1); y3(1); y4(1); y5(1); y6(1)];
% Set solver options
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8); % Adjust tolerances as needed
% Initialize function call counter
ncall = 0;
[t, u] = ode23t(@pde, t, u0, options);
i1 = i(dx1<=i & i<=dx2);
i2 = i(i<dx1 | dx2<i);
coeff = zeros(n,1);
coeff(i1) = coeff2;
coeff(i2) = coeff1;
coeffsin = zeros(n,1);
coeffsin(i1) = coeff4;
coeffsin(i2) = coeff3;
omega = zeros(n,1);
omega(i1) = omega2;
omega(i2) = omega1;
g = zeros(n,1);
g(i1) = g2;
g(i2) = g1;
% Reshape u into column vectors
u = reshape(u,n,6);
y1 = u(:,1); % As
y2 = u(:,2); % Ap
y3 = u(:,3); % Aa
y4 = u(:,4); % dAa/dt
y5 = u(:,5); % Phia
y6 = u(:,6); % dPhia/dt
% BCs
y1(1) = pulse(t);
y2(n) = 5e-3;
y1x = -super(xl, xu, n, y1, 1)';
y2x = super(xl, xu, n, y2, -1)';
% Derivatives in t
y1t = y1x - y2.*y3.*coeff;
y1t(1) = 0; % boundary condition
y2t = y2x + y1.*y3.*coeff;
y2t(n) = 0; % boundary condition
y3t(i) = y4; % y3t = dAa/dt
y5t(i) = y6; % y5t = dPhia/dt y6t = d2Phia/dt^2
y4t = (-g.*y1.*y2.*coeffsin - (2*gamma_k*y4 + 2*y3.*omega.*y6) + y3.*y6.*y6 - TT.*y3);
y6t = (-g.*y1.*y2.*coeffsin + 2*omega.*(y4 + gamma_k.*y3) - 2*y6.*(y4 + y3.*gamma_k))/y3;
% Combine derivatives into ut column vector
ut = [y1t; y2t; y3t; y4t; y5t; y6t];
% Increment calls to pde
ncall = ncall + 1;
end
y1_solution = u(:,1:n)';
figure;
plot(z, y1_solution);
xlabel('Spatial Domain (z)');
ylabel('Stokes Field (y1)');
title('Stokes Field Profile');
grid minor;
end
