clear all
close all

a.W21 = 121; %s^-1
a.W31 = 118; %s^-1
a.W32 = 2000; %s^-1
a.W41 = 51; %s^-1
a.W42 = 57; %s^-1
a.W43 = 1.7e5; %s^-1
a.h = 6.62e-34; %Js
a.sig_p = 4.2e-20; %cm^2
a.f1_p = 0.24;
a.f2_p = 0.22;
a.sig_l = 3.1e-20; %cm^2
a.f1_l = 0.0215;
a.f2_l = 0.2090;
a.lambda_p = 1532e-7; %cm         CHANGE THESE TO CM
lambda_l = 1645e-7; %cm
Ner = 0.7e20; %cm^3
a.Cup = 2e-18; %cm^3s^-1
a.A = 0.012; %cm^2
a.L_cry = 3; %cm
a.L_res = 9; %cm
a.c = 3e10; %cms^-1
a.Pin = 60; %w
laser_mode_area = 0.05; %cm^3

a.Toc = 0.99; %for the first time, different second time
a.phi_sp = 0;

%%%%Pout = phi*Toc

N0 = [0.7e20,0,0,0,0];

Er_Pop = @(t,N) laser_model(N,a);

[t,N] = ode45(Er_Pop,[0,10e-3],N0);
%run ode solver first, with N1 = 0.7e20, other 0 for 10ms
N1 = N(length(t),:);
a.Toc = 0.6;
N1(5) = 10e10;
[t,N] = ode45(Er_Pop,[0,100e-7],N1);
%use these results as input, run again for 100ns

plot(t,N(1),t,N(2),t,N(3),t,N(4));





function Ndot = laser_model(N,a)

%unpack structure a
phi = N(5);
sig_l = a.sig_l;
f1_l = a.f1_l;
f2_l = a.f2_l;
Pin = a.Pin;
lambda = a.lambda_p;
sig_p = a.sig_p;
f1_p = a.f1_p;
f2_p = a.f2_p;
h = a.h;
c = a.c;
A = a.A;
Toc = a.Toc;
L_res = a.L_res;
L = a.L_cry;
W21 = a.W21;
W31 = a.W31;
W41 = a.W41;
Cup = a.Cup;
W42 = a.W42;
W32 = a.W32;
W43 = a.W43;
phi_sp = a.phi_sp;

D = c*(exp(2*sig_l*L*(f2_l*N(2)-f1_l*N(1)))-1);
Rp = Pin*lambda*((1-exp(sig_p*L*(f2_p*N(2)-f1_p*N(1))))/(h*c*L*A));
rc = -c*log(1-Toc)/(2*L_res);

Ndot(1) = -Rp + phi*D/(2*L) + N(2)*W21 + N(3)*W31 + N(4)*W41 + Cup.*N(2).^2;
Ndot(2) = Rp - phi*D/(2*L) - N(2)*W21 + N(3)*W32 + N(4)*W42 - 2*Cup.*N(2).^2;
Ndot(3) = -N(3)*W31-N(3)*W32+N(4)*W43;
Ndot(4) = -N(4)*W41-N(4)*W42-N(4)*W43+Cup.*N(2).^2;
Ndot(5) = phi*D/(2*L_res) - rc*phi + phi_sp;

Ndot = transpose(Ndot);

end