% Nutsinee Kijbunchoo March 14, 2018
% Fabry-Perot cavity simulator
% Goal: To calculate transmitted power as a fuction of cavity length
% Use LHO VOPO parameters (March 28)

clear all
close all

R1 = .98 ; % green input coupler reflectivity (M1)
%R1 = .875 ; % green input coupler reflectivity (M1)
R2 = .9964 ; % final mirror in the bowtie path (M4)
T1 = 1-R1 ; 
T2 = 1-R2 ; 
t1 = sqrt(T1);
t2 = sqrt(T2);
r1 = sqrt(R1);
r2 = sqrt(R2);

c = 299792458;

lcr = 532e-9 ; %pump light
%lsb = c/(80e6);

fcr = 0;
fsb = 80e6;

%n_ktp = 1.8296; % crystal index of refraction for 1064. For green it's 1.8868
%d_crys = 10e-3; %SHG crystal length
%l = 0.05-0.01; %SHG Cavity length - crystal length
%L = l+(n_ktp*d_crys); %cavity total path length for 1064
%L = 0.279 ; %OPO optical path length from Stefszky's table 7.1
L = 2*0.175; %this makes FSR closer to the vopo final design document
%this is not the correct L. Calculate backward using 

FSR = c/(L); %Hz would have been 2L for linear cavity

%gamma = 0.05/2; %arbituary
gamma =0.094;
%gamma = 0.1;
%delt_l = logspace(-532e-9,5,532e-9);
delt_l = -600e-9:1e-10:600e-9;
%delt_l = 0;

%Ein = sqrt(7.4e-3); %OPO threshold power
Pin = 3e-3; %input power Watt
Ein = sqrt(Pin);

%Fin = 

%% Circ
%E circulated carrier
Ecirccr = Ein.*t1./(1-r1.*r2.*exp(-2.*i.*(pi.*fcr./FSR + 2*pi.*delt_l./lcr) ));

%E circulated SB
Ecircsb = i*gamma/2.*Ein.*t1./(1-r1.*r2.*exp(-2.*i.*(pi.*fsb./FSR + 2.*pi.*delt_l./lcr) ));
neg_Ecircsb = i*gamma/2.*Ein.*t1./(1-r1.*r2.*exp(-2.*i.*(-pi.*fsb./FSR + 2.*pi.*delt_l./lcr) ));

%% Trans
%E trans carrier
Etcr = Ein*t1.*t2.*exp( -i*(pi*fcr/FSR + 2*pi*delt_l/lcr) )...
        ./ (1-r1*r2*exp(-2*i*(pi*fcr/FSR + 2*pi*delt_l/lcr) ) );

%E trans sideband (+)
Etsb = Ein*i*gamma/2.*t1.*t2.*exp( -i.*(2.*pi.*delt_l./lcr))...
        ./ (1-r1.*r2.*exp(-2.*i.*(pi.*fsb./FSR + 2.*pi.*delt_l./lcr) ) );
%E trans sideband (-)
neg_Etsb = Ein*i*gamma/2.*t1.*t2.*exp( -i.*(2.*pi.*delt_l./lcr))...
        ./ (1-r1.*r2.*exp(-2.*i.*(-pi.*fsb./FSR + 2.*pi.*delt_l./lcr) ) );

%% Refl
%E reflected carrier
Ercr = -r1.*Ein + t1.*r2.*Ecirccr.*exp(-2.*i.*(-pi.*fcr./FSR + 2.*pi.*delt_l./lcr));
Ersb = -r1.*i*(gamma/2).*Ein + t1.*r2.*Ecircsb.*exp(-2*i*(pi.*fsb./FSR + 2.*pi.*delt_l./lcr));
neg_Ersb = -r1.*i*(gamma/2).*Ein + t1.*r2.*neg_Ecircsb.*exp(-2*i*(-pi.*fsb./FSR + 2.*pi.*delt_l./lcr));

%Ersb = -r1.*gamma.*Ein + t1.*Ecircsb;
%neg_Ersb = -r1.*gamma.*Ein + t1.*neg_Ecircsb;

%% Conjugate
Etcr_cj = conj(Etcr);
Etsb_cj = conj(Etsb);
neg_Etsb_cj = conj(neg_Etsb);

neg_Etsb_cj = conj(neg_Etsb);
Ecirccr_cj = conj(Ecirccr);
Ecircsb_cj = conj(Ecircsb);

neg_Ecircsb_cj = conj(neg_Ecircsb);
Ercr_cj = conj(Ercr);
Ersb_cj = conj(Ersb);
neg_Ersb_cj = conj(neg_Ersb);


%PRF = exp(i*-0.6)*((Etcr_cj.*neg_Etsb) + (Etcr_cj.*Etsb) + (Etcr.*neg_Etsb_cj) + (Etcr.*Etsb_cj));

%remove the mid two terms (why? don't remember...)
%PRF = exp(i*-0.6)*((Etcr_cj.*neg_Etsb) + (Etcr.*Etsb_cj)); %exp(i*theta) gives phase offset

PRF = ((Ercr_cj.*neg_Ersb) + (Ercr.*Ersb_cj));
PDC = (Ercr.*Ercr_cj)+(Ersb.*Ersb_cj)+(neg_Ersb.*neg_Ersb_cj);

%calculate phase where Q is flat at 0
a=log(0.00001./((Ercr_cj.*neg_Ersb) + (Ercr.*Ersb_cj)))/i;
%zero at a(6000)

Psb = (Ersb.*Ersb_cj)+(neg_Ersb.*neg_Ersb_cj); %sideband power
Pcirc = (Ecirccr.*Ecirccr_cj); %Power inside cavity. SB power negligible

Ptrans = (Etcr.*Etcr_cj)+(Etsb.*Etsb_cj)+(neg_Etsb.*neg_Etsb_cj);

figure(30)
subplot(3,1,1)
plot(delt_l, PDC, 'r')
xlabel('Frequency (Hz)'), ylabel('Power (W)')
grid on
legend('DC power')

subplot(3,1,2)
plot(delt_l, real(PRF), 'b') %plot RF power real part
%plot(delt_l, Psb, 'b') %plot RF power real part

xlabel('Cavity Length (m)'), ylabel('Power (W)')
grid on
legend('RF power (Re or I)')

subplot(3,1,3)
%plot(delt_l, imag(PRF), 'k') %plot RF power imag part
plot(delt_l, Psb, 'k') %plot RF power imag part

grid on
legend('RF power (Im or Q)')

figure(31)
subplot(2,1,1);
plot(delt_l, Pcirc);
grid on;
xlim([-1e-8 1e-8])
legend('circulated power (W)');
subplot(2,1,2);
plot(delt_l, imag(1./(1-r1.*r2.*exp(-2.*i.*(pi.*fcr./FSR + 2*pi.*delt_l./lcr))))); %Img part of Pcirc
grid on;
xlim([-1e-8 1e-8])
legend('phase (rad)');


% To-do 
% calculate linewidth


