%clear all
%close all

%SQZ measured with PSL LO Jan 10

%% Losses
% https://git.ligo.org/haocun.yu/lho_squeezing/wikis/squeezer-budget

theta = 10; %50; %rms phase fluctuations (mrad)

esc_l = 0.9863; %OPO escape efficiency (LHO)
prop_l = 0.98*0.96; %Propagation loss in HD path
% prop_l = 0.5;
pd_l = 0.99; %PD quantum efficiency

vis = 0.97; %homodyne fringe visibility
hd_l = vis^2; %effective homodyne loss

total_l = esc_l*prop_l*pd_l*hd_l; %Total loss
total_l_90 = esc_l*prop_l*pd_l*hd_l*0.9; %Total loss
total_l_80 = esc_l*prop_l*pd_l*hd_l*0.8; %Total loss
total_l_70 = esc_l*prop_l*pd_l*hd_l*0.7; %Total loss

%total_l = esc_l*0.975; %ANU
%total_l = 0.95
% total_l = 1;

threshold_P = 17; %threshold power in mW

%%
%Fractional Power to thershold
x = linspace(0,1,5000);

%Classical parametric gain
nlg = ((1+sqrt(x))./(1-x)).^2;
nlg_deamp = ((1-sqrt(x))./(1-x)).^2;

%Sqz/ anti-sqz (output variance)
asqz = 1 + total_l*4*sqrt(x)./((1-sqrt(x)).^2);
sqz = 1 - total_l*4*sqrt(x)./((1+sqrt(x)).^2);

asqz_90 = 1 + total_l_90*4*sqrt(x)./((1-sqrt(x)).^2);
sqz_90 = 1 - total_l_90*4*sqrt(x)./((1+sqrt(x)).^2);

asqz_80 = 1 + total_l_80*4*sqrt(x)./((1-sqrt(x)).^2);
sqz_80 = 1 - total_l_80*4*sqrt(x)./((1+sqrt(x)).^2);

asqz_70 = 1 + total_l_70*4*sqrt(x)./((1-sqrt(x)).^2);
sqz_70 = 1 - total_l_70*4*sqrt(x)./((1+sqrt(x)).^2);

%Sqz/ anti-sqz with phase fluctuations
sqz_ph = sqz*cos(theta*1e-3).^2 + asqz*sin(theta*1e-3).^2;
asqz_ph = asqz*cos(theta*1e-3).^2 + sqz*sin(theta*1e-3).^2;

sqz_ph_90 = sqz_90*cos(theta*1e-3).^2 + asqz_90*sin(theta*1e-3).^2;
asqz_ph_90 = asqz_90*cos(theta*1e-3).^2 + sqz_90*sin(theta*1e-3).^2;

sqz_ph_80 = sqz_80*cos(theta*1e-3).^2 + asqz_80*sin(theta*1e-3).^2;
asqz_ph_80 = asqz_80*cos(theta*1e-3).^2 + sqz_80*sin(theta*1e-3).^2;

sqz_ph_70 = sqz_70*cos(theta*1e-3).^2 + asqz_70*sin(theta*1e-3).^2;
asqz_ph_70 = asqz_70*cos(theta*1e-3).^2 + sqz_70*sin(theta*1e-3).^2;

%Sqz/anti-sqz (in dB)
sqz_db = 10*log10(sqz_ph);
asqz_db = 10*log10(asqz_ph);

sqz_db_90 = 10*log10(sqz_ph_90);
asqz_db_90 = 10*log10(asqz_ph_90);

sqz_db_80 = 10*log10(sqz_ph_80);
asqz_db_80 = 10*log10(asqz_ph_80);

sqz_db_70 = 10*log10(sqz_ph_70);
asqz_db_70 = 10*log10(asqz_ph_70);

%Measurement
sqz_mea = [-2.5, -2.34, -2.3];
Jan14_sqz_mea = [-2.4, -2];
asqz_mea = [4.4, 4, 3.7 ];
Jan14_asqz_mea = [4.5, 3.5]
nlg_mea = [2.1, 1.98, 1.784];
Jan14_nlg_mea = [2.1, 1.9];

%% Plots
% Pump power vs. gain

% figure(1)
% semilogy(x,nlg,...
% x,nlg_deamp)
% hold on
% xlabel('Fractional power to threshold (P/P_{threshold})')
% ylabel('Classical Parametric Gain')
% grid on
% hold off

%% P/Pth vs. SQZ/ASQZ
% figure(2)
% plot(x,sqz_db,'b-',...
% x,asqz_db,'r-')
% hold on
% xlabel('Fractional power to threshold (P/P_{threshold})')
% ylabel('Sqz/ Anti-Sqz (dB rel to SN)')
% grid on
% hold off

%% NLG vs. SQZ/ASQZ
figure(3)
plot(nlg,sqz_db_90,'r--',...
nlg,sqz_db_80,'b--',...
nlg,sqz_db_70,'k--',...
nlg,sqz_db,'g--',...
nlg,asqz_db_90,'r--',...
nlg,asqz_db_80, 'b--',...
nlg, asqz_db_70, 'k--',...
nlg, asqz_db, 'g--','LineWidth' ,1.5)

% axis([1, 100, -20, 30])
xlim([1 5])
hold on
plot(nlg_mea,sqz_mea,'mo',...
nlg_mea,asqz_mea,'mo',...
Jan14_nlg_mea, Jan14_sqz_mea, 'mx',...
Jan14_nlg_mea, Jan14_asqz_mea, 'mx', 'LineWidth' ,2)
xlabel('Classical Parametric Gain')
ylabel('Sqz/ Anti-Sqz (dB rel to SN)')
%yticks([-5 -2.5 0 2.5 5 7.5])
set(gca,'FontSize',18);
hTitle = title('SQZ measurement with PSL LO Jan 9 (o), Jan14 (x), 2019');
set(hTitle,'FontSize',16)
legend({'efficiency = 0.78, phase noise = 10 mrad', 'efficiency = 0.69, phase noise = 10 mrad', 'efficiency = 0.61, phase noise = 10 mrad', 'efficiency = 0.86 (expected)'}, 'FontSize',14, 'Location', 'Northwest')

grid on
grid minor
hold off

%% Pump power vs. SQZ/ASQZ
% figure(4)
% plot(x*threshold_P,sqz_db,'b-',...
% x*threshold_P,asqz_db,'r-')
% hold on
% xlabel('Pump power (mW)')
% ylabel('Sqz/ Anti-Sqz (dB rel to SN)')
% xlim([0 20])
% grid on
% hold off

%% Threshold power vs. nlg

% figure(5)
% semilogy(x*threshold_P,nlg,...
% x*threshold_P,nlg_deamp)
% hold on
% xlabel('Pump power (mW)')
% ylabel('Classical Parametric Gain')
% xlim([0 280])
% grid on
% hold off