clc
clear all
tspan = [0 5e-9];
y0 = [0.2,0.2,0.2];
[t,y] = ode45(@rate_fano1,tspan,y0);
size(t);
t=t*1e9;
y = abs(y);
params
pt=2.*epsilon0.*ref_index.*c.*(abs(y(:,2)+p.*y(:,3))).*2;
pc=2.*epsilon0.*ref_index.*c./(gamma_p./gamma_c).*(abs(y(:,3))).*2;
figure(1)
plot(t, pt, t, pc,'.-');
xlabel('time [ns]','FontSize',14);
ylabel('Arbitrary units','FontSize',14);
set(gca,'FontSize',14);
legend('pt','pc', 'Location','SE')
x=[0 deltac/gamma_T];
figure(2)
plot(x, pt);
xlabel('time [ns]','FontSize',14);
ylabel('Arbitrary units','FontSize',14);
set(gca,'FontSize',14);
legend('pt','pc', 'Location','SE')
figure(3)
plot(t, y(:,1)./N_0);
xlabel('time [ns]','FontSize',14);
ylabel('Arbitrary units','FontSize',14);
set(gca,'FontSize',14);
legend('N', 'Location','SE')
figure(4)
plot(t, y(:,2));
xlabel('time [ns]','FontSize',14);
ylabel('Arbitrary units','FontSize',14);
set(gca,'FontSize',14);
legend('A+', 'Location','SE')
figure(5)
plot(t, y(:,3));
xlabel('time [ns]','FontSize',14);
ylabel('Arbitrary units','FontSize',14);
set(gca,'FontSize',14);
legend('A-', 'Location','SE')
function yx = rate_fano1(t,y)
params
r_R=(-p*gamma_c)/(1i*(deltac)+gamma_T);
sigmaa=((2.*epsilon0.*ref_index.*ref_indexg)./(hbar.*w_s)).*(1+(abs(r_R))).*(1+(abs(r_R)))./(conf.*g_n.*(y(1)-N_0))
yx=zeros(3,1);
yx(1)=I./(e.*V_a)-y(1)./tau1-(V_g.*g_n.*(y(1)-N_0).*sigmaa.*(abs(y(2)))^2)./V_m;
yx(2)=1/2*(1-1i*henry).*conf.*V_g.*g_n.*(y(1)-N_ss).*y(2)+gamma_L.*((y(3)./(r_R)-y(2)));
yx(3)=(-1i.*deltac-gamma_T).*y(3)-p.*gamma_c.*y(2);
c = 2.99792458e8;
e = 1.6021766208e-19;
h = 6.626068e-34;
h_eV = h/e;
hbar = h/(2*pi);
hbar_ev = hbar/e;
p=-1;
L = 5e-6;
A = 0.21e-12;
conf = 0.5;
conf_NC =0.3;
V_a=A*L;
V_m = V_a/conf;
V_NC=0.24e-18;
ref_index = 3.5;
ref_indexg=3.5;
V_g = c/ref_index;
tau1 = 0.5e-9;
g_n = 5e-21;
N_0=1e24;
N_ss=1e24;
I =5e-3;
henry_i=1000;
henry = 1;
Qt = 500;
Qi = 14300;
Qp = 10000;
lambda0 = 1554e-9;
deltalambda = 0.01e-9;
lambda_start = 1540e-9;
lambda_end = 1560e-9;
lambda = lambda_start:deltalambda:lambda_end;
omega = 2*pi*c./lambda;
omega0 = omega(find(lambda==lambda0));
gamma_i = omega0/(2*Qi);
gamma_p = omega0/(2*Qp);
gamma_T = omega0/(2*Qt);
gamma_c = gamma_T-gamma_i-gamma_p;
gamma_L=V_g/(2*L);
rB = 0.2 ;
tB = sqrt(1-rB^2);
w_r=omega0;
w_s=1.2121*10^15;
w_c=1.2125*10^15;
epsilon0 = 8.854187817e-12;
rho=(2*epsilon0*ref_index*c)./gamma_c.*hbar_ev.*w_r;
r_L=1;
deltac=1.2125*10^15-1.2121*10^15;