clear
close all

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

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

filename = 'uc_raw_dat_tt.xlsx';
sheet = 1;
xlRange = 'B2:J25';

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

C = {'1','5','10','20','30','40','50','60','70','80','90','90ii','100','110','120','130','140','150','160','160ii','170','180','180ii','190'};
load('uc_noise.mat');
load('uc_cdata.mat');

pump_power = num(:,6);
normed_counts= zeros(size(pump_power));
[s, snc, nnc, Pi_m, ir_standdev, tot_count_s, tot_count_n, tot_count_snc, tot_count_nnc] = deal(normed_counts);
 
for k = 1:24
    
combined = sct.(['p' char(C(k))]).combined;

%% FOR IR CALIBRATION
CIR_Start = combined(num(k,1),4);
CIR_Stop  = combined(num(k,5),4);
PIR_Start = num(k,7);
PIR_Stop  = num(k,8);

In_Pfit = polyfit([CIR_Start CIR_Stop],[PIR_Start PIR_Stop],1);
IR_PWRs = polyval(In_Pfit,combined(num(k,1):num(k,5),4));

if num(k,1) ~= 1
    fill = IR_PWRs(1).*ones([1 (num(k,1)-1)]);
    IR_PWRs = [fill'; IR_PWRs];
else 
end
ir_standdev(k) = std(IR_PWRs(num(k,2):num(k,4))); % uncertainty in IR counts. 

%% Signal and Noise
B_raw_c = combined(num(k,2):num(k,3),2); %(k,2)-(k,3) means includes init_peak
B_no_cor = B_raw_c;
tot_count_snc(k) = 5*sum(B_no_cor); % 5s time bins included;
B_raw_c = B_raw_c./(1- B_raw_c * 43e-9); % DEAD TIME ON SIGNAL
tot_count_s(k) = 5*sum(B_raw_c); % 5s time bins included;
nnc(k) = n(k); %Noise No Dead Time Correction
tot_count_nnc(k) = 25*sum(nnc(k)); % 25s integration time factored back in. 
n(k) = n(k)./(1 - n(k) * 43e-9); % DEAD TIME ON BACKGROUND
tot_count_n(k) = 25*sum(n(k)); % 25s integration time factored back in. 
Pi_m(k) = mean(IR_PWRs(num(k,2):num(k,3)));
s(k) = mean(B_raw_c) * num(k,9); % ND Scaling
snc(k) = mean(B_no_cor)*num(k,9); %ND Scaling SIGNAL NO DEAD TIME CORRECTION

end 

%% FACET LOSS
Pi_m = Pi_m./facet_IR_trans; % FACET LOSS ON IR
s = s./facet_bl_trans; % FACET LOSS ON SIGNAL
n = n./facet_ye_trans; % FACET LOSS ON NOISE 

%% Removal of data files where data files show experimental error. 
pump_power([11 16 18 19 22]) = [];
%normed_counts([11 16 18 19 22]) = [];
s([11 16 18 19 22]) = [];
n([11 16 18 19 22]) = [];
ir_standdev([11 16 18 19 22]) = [];
Pi_m([11 16 18 19 22]) = [];
tot_count_s([11 16 18 19 22]) = [];
tot_count_n([11 16 18 19 22]) = [];
snc([11 16 18 19 22]) = [];
nnc([11 16 18 19 22]) = [];
tot_count_snc([11 16 18 19 22]) = [];
tot_count_nnc([11 16 18 19 22]) = [];
pump_power = pump_power/facet_ye_trans;
n = n';
%% Efficiency 

nopt = 0.2*0.96; 
ndet = 0.86;
T_i = 12.5e-9;
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/8.7749e-9);
sraw = (snc./(Pi_m/8.7749e-9))/160e6;
n_input = (nnc./80e6)./sraw;

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

Err_s = (T_i./T_p).*((hbar*w)./(Pi_m.*ndet.*nopt)).*ds;
Err_n = -1.*(T_i./T_p).*((hbar*w)./(Pi_m.*ndet.*nopt)).*dn;
Err_Pm = ((n-s)./Pi_m.^2).*(T_i./T_p).*((hbar*w)./(ndet.*nopt)).*dPm;
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;
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.*(8.7749e-9)./(nnc.* Pi_m); 
snr_Err_nnc = dnnc.*(snc*8.7749e-9)./((nnc.^2).* Pi_m);
snr_Err_Pm = dPm.*(snc*8.7749e-9)./(nnc.* (Pi_m.^2));
           
srw_Err_snc = (8.7749e-9./(160e6 .* Pi_m)) .* dsnc;
srw_Err_Pm  = ((-snc.*8.7749e-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 Efficiency 
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 40 80 120 160 200];
axes1.YTick = [0 2 4 6 8 10];
axes1.FontSize = 32;
axes1.XLim = [0 200];
axes1.YLim = [0 10.4];
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
C = linspecer(3); 
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 10 20 30 40];
axes2.FontSize = 32;
axes2.XLim = [0 200];
axes2.YLim = [0 45];
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

fig3 = figure; 
set(fig3,'position',[230 250 570 510]);
set(fig3,'resize','off');
scatter(pump_power,n_input,80,C(1,:),'filled','v'); 
axes3 = gca;
axes3.LineWidth = 2;
axes3.XTick = [0 40 80 120 160 200];
axes3.YTick = [0 0.3 0.6 0.9 1.2 1.5 1.8];
axes3.FontSize = 20;
axes3.XLim = [0 200];
hold on; 
errorbar(pump_power,n_input, Err_mu1,'LineStyle','none','LineWidth',1.3,'Color',C(1,:));
xlabel('Pump Power [mW]','FontSize',22);
ylabel('Input Photons Per Pump Pulse','FontSize',22); 
axis square
box on

%% Save temp files 

save('uc_mu1.mat','pump_power','n_input', 'Err_mu1');