% Nutsinee Kijbunchoo Aug23, 2018
% bowtie cavity simulator
% Goal: To calculate transmitted power as a fuction of cavity length
% Use LHO VOPO parameters 

%clear all
%close all

R1 = .98 ; % green input coupler reflectivity (M1)
%R1 = 0.87298077; %for red, testing code against sheila's alog40500

R2 = .99998293 ; %mirror#11 final mirror in the bowtie path. E1800011
%R2 = 0.99681; %tweaked
%R2 = 0.99818636; %red, testing code against sheila's alog40500
R3 = 0.99998067; %mirror#6
R4 = 0.99998159; %mirror#8
%Rsum = R2+R3+R4 

T1 = 1-R1; 
T2 = 1-R2; 
T3 = 1-R3;
T4 = 1-R4;

t1 = sqrt(T1);
t2 = sqrt(T2);
t3 = sqrt(T3);
t4 = sqrt(T4);
r1 = sqrt(R1);
r2 = sqrt(R2);
r3 = sqrt(R3);
r4 = sqrt(R4);
rloss = 0.99467; %exp(2*alpha*round trip length) for length dependent loss. See Siegman page 414

c = 299792458;

lcr = 532e-9 ; %pump light
%lcr = 1064e-9 ; 

%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
%nktp = 1.8864; %crystal index for green
%dcrys = 10e-3;
%l = 0.05-0.01; %SHG Cavity length - crystal length
%l = (2*0.175)-dcrys;
%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; %meter, this makes FSR closer to the vopo final design document
%L = l+(nktp*dcrys);

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

%gamma = 0.05/2; %arbituary
gamma =0.094; %alog 41622
%gamma = 0.1;
delt_l = -600e-9:1e-10:600e-9;
%delt_l = 0;

%Ein = sqrt(7.4e-3); %OPO threshold power
Pin = .36; 
Ein = sqrt(.7*Pin); %
%Fin = 

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

%E circulated SB
Ecircsb = i*gamma/2.*Ein.*t1./(1-r1.*rloss.*exp(-2.*i.*(pi.*fsb./FSR + 2.*pi.*delt_l./lcr) ));
neg_Ecircsb = i*gamma/2.*Ein.*t1./(1-r1.*rloss.*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-r_all.*exp(-2*i*(pi*fcr/FSR + 2*pi*delt_l/lcr) ).*exp(-abs) );

Etcr = Ecirccr.*t2.*exp( -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.*rloss.*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.*rloss.*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(32)
subplot(3,1,1)
plot(delt_l, PDC, 'r')
title('OPO Refl Calculated (532)');
xlabel('Frequency (Hz)'), ylabel('Power (W)')
grid on
legend('DC Refl power (rloss = 0.99467)')

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
plot(delt_l, Ptrans, 'k') %plot RF power imag part

grid on
%legend('RF power (Im or Q)')
legend('DC trans power')


figure(33)
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


