old computational project

clear variables
close all
load ndyag
solar = csvread('ASTMG173_Solar_data.csv',2,0,'A3..B2004');
%808nm occurs at row 649
solar_wav = solar(:,1);
solar_em = solar(:,2);
%Constants
c = 3e8;
k = 1.3806505e-23;
h = 6.63e-34;
hbar = h/2*pi;

%Let's try and create the Pedrotti cavity first!
R_1 = 2;                 %Radius of mirror one (m)
R_2 = [1 2 3 10];        %Radius of mirror two (m)
d = 0.5;                 %Cavity length (m)
Pp = 700;                %Input power (full spectrum) (W)
r_1 = 0.99;              %Reflectivity of first mirror (i.e. 99%)
r_2 = 0.95;              %Reflectivity of second mirror (i.e. 90%)
ndyag_wav = NdYAG(:,1);  %Wavelengths of Nd:YAG spectrum
ndyag_abs = NdYAG(:,2);  %Corresponding absorption coefficients
irr_808 = solar_em(649)/sum(solar_em(:));   %Total solar irradiance required
pump_wav = 808e-9;       %Wavelength of pump (m)
pump_fre = c/pump_wav;   %Frequency of pump (Hz)
N = 10e20;               %Number of Nd ions/cm^3
sig = 7.6e-19;           %Stim emission cross section for 808nm transition (cm^2)
tau = 230e-6;            %Fluorescent lifetime of Nd:YAG (s)
gamma = (1-r_1)+(1-r_2); %Example fractional round trip loss
r = 1.5e-3;              %Radius of Nd:YAG crystal
L = 20e-3;               %Length of Nd:YAG crystal
vol = pi*r^2*L;          %Volume of Nd:YAG crystal (cm^3)
Pp_thr = (gamma*h*pump_fre*pi)/(2*sig*tau);
Pp_thr_full = Pp_thr/irr_808; %Pump threshold of solar energy if only relying on 808nm part
N_thr = gamma/(2*sig*L);
Rp_thr = N_thr/tau;      %Threshold pump rate
out_wav = 1064e-9;       %Wavelength of laser
out_fre = c/out_wav;     %Frequency of laser
cols = {'r-','b-','y-','g-','p-'};
mirrorcols = {'r--','b--','y--','g--','p--'};
I_0 = 1;                 %Maximum intensity of laser


%stable = zeros(1, length(R_2));
z_0 = zeros(1, length(R_2));                %Empty array for z_0
z_1 = zeros(1, length(R_2));                %Empty array for z_1
z_2 = zeros(1, length(R_2));                %Empty array for z_2

for i = 1:length(R_2)
stable(1,i) = (1-(d/R_1))*(1-(d/R_2(i)));
z_0(1,i) = sqrt((d*(R_1-d)*(R_2(i)-d)*(R_1+R_2(i)-d))/(R_1+R_2(i)-2*d)^2); %Calculates z_0 for a given R_2
z_1(1,i) = -R_1/2 + sqrt((R_1^2-(4*z_0(1,i)^2))/4);                        %Calcualtes z_1 for a given R_2
z_2(1,i) = R_2(i)/2 - sqrt((R_2(i)^2-(4*z_0(1,i)^2))/4);                   %Calculates z_2 for a given R_2
end

for i = 1:length(z_0)
    z(i,:) = linspace(z_1(1,i),2,20000);
end

w_0sq = (z_0*out_wav)/pi;                       %Radius of beam waist (squared)
w_0 = sqrt(w_0sq);                              %Radius of beam waist
w = zeros(1,length(z));                         %For calculating beam radius at different locations
R = zeros(1,length(z));                         %For calculating wavefront radius of curvature at different locations
all_w = cell(1,length (R_2));                   %Stores radius arrays
for i = 1:length(w_0)
    rho(i,:) = linspace(0,2*w_0(1,i),20);       %Different values for distance from beam centre
end

I = zeros(length(rho),length(z));
all_I = cell(1,length(R_2));

for j = 1:length(R_2)
    rho_loop = rho(j,:);
    for i = 1:length(z)
        for k = 1:length(rho)
        w(1,i) = sqrt(w_0sq(1,j)*(1+(z(1,i)^2/z_0(1,j)^2)));    %Calculates radius along direction of propagation

        R(1,i) = z(1,i)*(1+(z(1,i)^2/z_0(1,j)^2));              %Calculates wavefront radius of curvature along direction of propagation
        I(k,i) = I_0*(w_0(1,j)/w(1,i))^2*exp(-2*rho_loop(k)^2/w(1,i)^2);%Calculates relative intenstiy along direction of propagation
        end
    end
    all_w{j} = w;
    all_w_down{j} = -w;
    all_I{j} = I;
    all_R{j} = R;
end

figure; hold on                 %For plotting irradiance profiles
for i = 1:length(R_2)
I = all_I{i};
plot(z(1,:),I(1,:),cols{i})
title('Irradiance profile of beam (\rho = 0)')
legend('R_2 = 1m','R_2 = 2m', 'R_2 = 3m','R_2 = 10m')
xlabel('Distance from beam waist (m)')
ylabel('Intensity relative to I_0')
end
y_z1 = [-1.5e-3 1.5e-3];
y_z2 = [-1.5e-3 1.5e-3];

figure; hold on                 %For plotting beam trajectories
for i = 1:length(w_0)
    x_z1 = [z_1(1,i) z_1(1,i)];
    x_z2 = [z_2(1,i) z_2(1,i)];
    plot(z(i,:),all_w{i},cols{i},z(i,:),-all_w{i},cols{i})
    plot(x_z1, y_z1,mirrorcols{i})
    plot(x_z2, y_z2,mirrorcols{i})
end

title('Beam waist of laser')
xlabel('Distance (m)')
ylabel('Beam radius (m)')
legend('R_2 = 1','R_2 = 2','R_2 = 5','R_2 = 10')

w0_ch = min(w_0);               %Chooses the minimum beam waist
theta_ff = out_wav/(pi*w0_ch);  %Finds far field divergence of minimum beam waist
Pp_thr = Pp_thr*(w0_ch)^2;      %Finds threshold pump power (808nm only)
Pp_thr_full = Pp_thr_full*(w0_ch)^2; %Finds threshold pump power (full spectrum)
Pp_fact = Pp/Pp_thr_full;       %Finds factor by which input power is over threshold
Pp_808 = Pp_thr*Pp_fact;        %Finds input pump power (808nm only)

Rp = Pp_808/(h*pump_fre*pi*w0_ch^2*L); %Finds pump rate
F = (2*L)/gamma*(Rp - Rp_thr);  %Intracavity photon flux
T = 1 - r_2;                    %Percentage of light exiting cavity as beam
wl = 4.5e-4;                    %Estimation of pump
Po = (F/2)*(1-r_2)*(pi*wl^2)*(h*out_fre); %Calculates output power