%clear all
%close all

%SQZ measured with PSL LO Jan 18

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

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

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 86%
total_l_90 = esc_l*prop_l*pd_l*hd_l*0.9; %Total loss 78%
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_low = esc_l*prop_l*pd_l*hd_l*1.1; %Total loss 95%

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

threshold_P = 12;  %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);

asqz_lowloss = 1 + total_l_low*4*sqrt(x)./((1-sqrt(x)).^2);
sqz_lowloss = 1 - total_l_low*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_hightheta = sqz*cos(high_theta*1e-3).^2 + asqz*sin(high_theta*1e-3).^2;
asqz_ph_hightheta = asqz*cos(high_theta*1e-3).^2 + sqz*sin(high_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_ph_lowloss = sqz_lowloss*cos(theta*1e-3).^2 + asqz_lowloss*sin(theta*1e-3).^2;
asqz_ph_lowloss = asqz_lowloss*cos(theta*1e-3).^2 + sqz_lowloss*sin(theta*1e-3).^2;


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

sqz_db_hightheta = 10*log10(sqz_ph_hightheta);
asqz_db_hightheta = 10*log10(asqz_ph_hightheta);

sqz_db_lowloss = 10*log10(sqz_ph_lowloss);
asqz_db_lowloss = 10*log10(asqz_ph_lowloss);

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);


%SQZ Measurement
sqz_with_atten = [-3.958, -3.22, -2.56];
asqz_with_atten = [6, 5, 4.6];

sqz_wo_atten = [-3.15 -2.58];
asqz_wo_atten = [5, 4.238];

sqz_3MHz_sat = [-2.53, -2.1];
asqz_3MHz_sat = [5, 4.8];


nlg_with_atten = [2.2811, 2.0283, 1.8165];
nlg_wo_atten = [1.962, 1.7316];
nlg_3MHz_sat = [1.9632, 1.7838];




%pump in measurement
%pump_in = [1.35 0.962]; %SB taken into account, use 48% total
%pump_in = [1.3 0.9]; %SB taken into account, use 50% total
nlg_mea_with_atten = [2.239, 2.0726, 1.8171, 2.2811, 2.0283, 1.8165];
pump_in_with_atten = [(2.6*.78), (2.1*.78), (1.6*.78), 2.67*.78, 2.1*.78, 1.6*.78]; 




%((sqrt((0.00547+6.8e-5)./(0.0015+6.8e-5)) + 1)/2).^2

%% 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(7)
plot1 = plot(nlg,sqz_db_90,'r--',...
    nlg,sqz_db_80,'b--',...
    nlg,sqz_db_70,'k--',...
    nlg,sqz_db,'g--',...
    nlg,sqz_db_hightheta,'g',...
     nlg,asqz_db_90,'r--',...
     nlg,asqz_db_80, 'b--',...
     nlg, asqz_db_70, 'k--',...
     nlg, asqz_db, 'g--',...
     nlg,asqz_db_hightheta,'g','LineWidth' ,1.5)

     
% axis([1, 100, -20, 30])
xlim([1 3.5])
hold on
plot2 = plot(nlg_with_atten,sqz_with_atten,'mo',...
     nlg_wo_atten, sqz_wo_atten, 'kx',...
     nlg_3MHz_sat, sqz_3MHz_sat,'rs',...
     nlg_with_atten,asqz_with_atten,'mo',...
     nlg_wo_atten, asqz_wo_atten, 'kx',...
     nlg_3MHz_sat, asqz_3MHz_sat, 'rs', '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 Mephisto LO Jan18, 2019');
set(hTitle,'FontSize',16)
legend(plot1, {'efficiency = 0.78, theta = 10 mrad', 'efficiency = 0.69, theta = 10 mrad', 'efficiency = 0.61, theta = 10 mrad', 'efficiency = 0.86, theta = 10mrad', 'efficiency = 0.86, theta = 150mrad'}, 'FontSize',14, 'Location', 'Northwest')
%legend(plot2, {'20dB attenuator ON', '20dB attenuator OUT', 'HD saturated'}, 'FontSize',14, 'Location', 'Northeast')

%legend({'efficiency = 0.86, phase noise = 250 mrad'}, '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(9)
semilogy(x*18, nlg, '--', 'LineWidth' ,1.5)
hold on     

semilogy(pump_in_with_atten, nlg_mea_with_atten, 'ro','LineWidth' ,2)

xlabel('Pump power (mW)')
ylabel('Classical Parametric Gain')
xlim([0 5])
ylim([0 5])
%yticks([1 2 3 4 5 6 7 8 9 10])

set(gca,'FontSize',18, 'YMinorTick', 'on');
hTitle = title('nlg measurement Jan18, 2019');
set(hTitle,'FontSize',16)
legend({'Threshold Power = 18'}, 'FontSize',14, 'Location', 'Northwest')

grid on
hold off