clear
close all

c = 299792458;
wl = 425.5e-9;
w = 2*pi*c/wl;
hbar = 1.0545718e-34;

facet_IR_trans = 1;
facet_bl_trans = 1;
facet_ye_trans = 1;

filename = 'dc_raw_dat_tt.xlsx';
sheet = 1;
xlRange = 'B2:J30';

[num,text,raw] = xlsread(filename,sheet,xlRange);

C = {'005' '010' '020' '020_2' '030' '040' '040_2' '050' '050_2' '060' '060_2'...
    '070' '070_2' '080' '080_2' '090' '090_2' '090_3' '090_4' '100' '100_2'...
    '110' '110_2' '110_3' '110_4' '120' '120_2' '120_3' '120_4'};

load('dc_noise.mat');

pump_power = num(:,6);
normed_counts= zeros(size(pump_power));
[s, snc, nnc, Pi_m, in_standdev, tot_count_s, tot_count_n, tot_count_snc, tot_count_nnc] = deal(normed_counts);

load('dc_cdata.mat');

for k = 1:29

combined = sct.(['p' char(C(k))]).combined;

% FOR IR CALIBRATION
CIn_Start = combined(num(k,1),2);
CIn_Stop  = combined(num(k,5),2);
PIn_Start = num(k,7);
PIn_Stop  = num(k,8);

In_Pfit = polyfit([151 CIn_Start CIn_Stop],[0 PIn_Start PIn_Stop],1);
In_PWRs = polyval(In_Pfit,combined(num(k,1):num(k,5),2));
if num(k,1) ~= 1
    fill = In_PWRs(1).*ones([1 (num(k,1)-1)]);
    In_PWRs = [fill'; In_PWRs];
else 
end
in_standdev(k) = std(In_PWRs(num(k,1):num(k,5))); % uncertainty in IR counts. 
sig_raw_c = combined(num(k,2):num(k,3),4); %(k,2)-(k,3) means includes init_peak
sig_no_cor = sig_raw_c;
tot_count_snc(k) = 5.*sum(sig_no_cor);
sig_raw_c = sig_raw_c./(1- sig_raw_c * 10e-6); % Dead time on signal 
tot_count_s(k) = 5.*sum(sig_raw_c); % 5s time bins included; 
nnc(k) = n(k);
tot_count_nnc(k) = 25*sum(nnc(k));
n(k) = n(k)./(1 - n(k) * 10e-6); % Dead time on background
tot_count_n(k) = 25*sum(n(k)); % 25s integration time factored back in. 
B_pnorm_c = (sig_raw_c - n(k))./(In_PWRs(num(k,2):num(k,3))/(0.7661)); % 2-3 means includes peak.
Pi_m(k) = mean(In_PWRs(num(k,2):num(k,3)));
m_B_pnorm_c = mean(B_pnorm_c);
m_B_pnorm_c = m_B_pnorm_c * num(k,9); % ND Scaling
normed_counts(k) = m_B_pnorm_c;
s(k) = mean(sig_raw_c) * num(k,9); % ND Scaling SIGNAL
snc(k) = mean(sig_no_cor)*num(k,9); %ND Scaling SIGNAL no dead time correction

end 

Pi_m = Pi_m./facet_bl_trans; % Facet loss on IR
s = s./facet_IR_trans; % Facet loss on signal
n = n./facet_IR_trans;
pump_power = pump_power/facet_ye_trans;
n = n';

% Efficiency 
nopt = 0.45*0.9*0.9*0.93*0.97; 
ndet = 0.095;
T_i = 300e-12;
T_p = 2.9250e-11;
Pi_m = Pi_m.*1e-9;
next2 = 100.*((s - n)./Pi_m).*((hbar *w)./(nopt*ndet))*((T_i)/(T_p));

% SNR
SNR = (snc./nnc)./(Pi_m/0.7661e-9);
sraw = (snc./(Pi_m/0.7661e-9))/160e6;
n_input = (nnc./80e6)./sraw;

% Errors
ds = sqrt(tot_count_s)./5;
dn = sqrt(tot_count_n)./25;
dPm = in_standdev.*1e-9;
dTi = 10e-12;
dTp = 2.5e-12;
dnd = 0.005;
dno = 0.005; 

Err_s = (T_i./T_p).*((hbar*w)./(Pi_m.*ndet.*nopt)).*ds; %Sig
Err_n = -1.*(T_i./T_p).*((hbar*w)./(Pi_m.*ndet.*nopt)).*dn; % Noise
Err_Pm = ((n-s)./Pi_m.^2).*(T_i./T_p).*((hbar*w)./(ndet.*nopt)).*dPm; % Power
Err_Ti = ((s-n)./Pi_m).*(1/T_p).*((hbar*w)./(ndet.*nopt)).*dTi;
Err_Tp = ((s-n)./Pi_m).*(T_i/T_p.^2).*((hbar*w)./(ndet.*nopt)).*dTp;
Err_nd = ((n-s)./Pi_m).*(T_i/T_p).*((hbar*w)./((ndet.^2).*nopt)).*dnd;% Detector
Err_no = ((n-s)./Pi_m).*(T_i/T_p).*((hbar*w)./(ndet.*(nopt.^2))).*dno;%  

Err_next = 100.*sqrt(Err_s.^2 + Err_n.^2 + Err_Pm.^2 + Err_Ti.^2 + ...
           Err_Tp.^2 +Err_nd.^2 + Err_no.^2);

dsnc = sqrt(tot_count_snc)./5;
dnnc = sqrt(tot_count_nnc)./25;
      
snr_Err_snc = dsnc.*(0.7661e-9)./(nnc.* Pi_m); 
snr_Err_nnc = dnnc.*(snc*0.7661e-9)./((nnc.^2).* Pi_m);
snr_Err_Pm = dPm.*(snc*0.7661e-9)./(nnc.* (Pi_m.^2));

srw_Err_snc = (0.7661e-9./(160e6 .* Pi_m)) .* dsnc;
srw_Err_Pm  = ((-snc.*0.7661e-9)./(160e6 .*(Pi_m.^2))).*dPm;

Err_snr = sqrt(snr_Err_snc.^2 + snr_Err_nnc.^2 + snr_Err_Pm.^2);
Err_sraw= sqrt(srw_Err_snc.^2 + srw_Err_Pm.^2);

mu1_Err_nnc = (1./(sraw.*80e6)) .* dnnc;
mu1_Err_srw = (nnc./((sraw.^2).*80e6)) .* Err_sraw;

Err_mu1 = sqrt(mu1_Err_nnc.^2 + mu1_Err_srw.^2);

% Plot EFF

C = linspecer(3); 
fig1 = figure; 
set(fig1,'position',[230 250 570 510]);
set(fig1,'resize','off');
scatter(pump_power,next2,150,C(1,:),'filled','v'); 
axes1 = gca;
axes1.LineWidth = 4.2;
axes1.XTick = [0 30 60 90 120 150];
axes1.YTick = [0 0.2 0.4 0.6 0.8 1.0 1.2];
axes1.XLim = [0 125];
axes1.YLim = [0 1.25];
axes1.FontSize = 32;
hold on; 
errorbar(pump_power,next2, Err_next,'LineStyle','none','LineWidth',2.2,'Color',C(1,:));
xlabel('Pump Power [mW]','FontSize',32);
ylabel('Efficiency [%]','FontSize',32); 
axis square
box on
a123 = [pump_power]\next2;
P = [a123;0];
pump_power_zero = [0; pump_power];
yf = polyval(P,pump_power_zero);
plot(pump_power_zero, yf,'color',C(1,:),'linewidth',3);

% Plot SNR

fig2 = figure; 
set(fig2,'position',[230 250 570 510]);
set(fig2,'resize','off');
scatter(pump_power,SNR,150,C(1,:),'filled','v'); 
axes2 = gca;
axes2.LineWidth = 4.2;
axes2.XTick = [0 40 80 120 160 200];
axes2.YTick = [0 20 40 60 80 100 120];
axes2.XLim = [0 125];
axes2.YLim = [0 130];
axes2.FontSize = 32;
hold on; 
errorbar(pump_power,SNR, Err_snr,'LineStyle','none','LineWidth',2.2,'Color',C(1,:));
xlabel('Pump Power [mW]','FontSize',32);
ylabel('SNR','FontSize',32); 
axis square
box on

% Plot mu_1

% struct1 values are mu_1 values obtained from upconversion matlab file 
struct1 = load('uc_mu1.mat'); 

fig3 = figure; 
set(fig3,'position',[230 250 570 510]);
set(fig3,'resize','off');
DFG = scatter(struct1.pump_power,struct1.n_input,150,C(1,:),'filled','d'); hold on
SFG = scatter(pump_power,n_input,150,C(2,:),'filled','v'); 
axes3 = gca;
axes3.LineWidth = 4.2;
axes3.XTick = [0 40 80 120 160 200];
axes3.FontSize = 32;
axes3.XLim = [0 200];
axes3.YLim = [0 1];

hold on; 
errorbar(pump_power,n_input, Err_mu1,'LineStyle','none','LineWidth',2,'Color',C(2,:));
errorbar(struct1.pump_power,struct1.n_input, struct1.Err_mu1,'LineStyle','none','LineWidth',2,'Color',C(1,:));
set(axes3,'yscale','log')
xlabel('Pump Power [mW]','FontSize',32);
ylabel('Input Photons Per Pulse','FontSize',32); 
axis square
box on
axes3.Layer = 'top';
[fig3, ~] = legend('SFG','DFG','location','northeast');
legend('boxoff');