% T.A.Wright et al (2020) 

% This script plots the pump-power-dependant down-conversion efficiency, 
% shown in Fig.4(c).
% The calculation describing the conversion is described within the
% manuscript; we sum the power at each wavelength using  BS-FWM spectra 
% obtained on an Optical Spectrum Analyser (OSA) and evaluate 
% nu_up = ((P_T - N_T) * omega_Sr)/ (P_Sr * D * omega_T). The spectra were
% saved by the OSA using a csv format. The data files required for this 
% plot are in the form 'wl/AXXXX.CSV', where XXXX denotes numbers between 
% 0000 and 00014. 

% figure object fig1 is shared with plot_fig4a_inset.m, plot_fig4b.m 
% and plot_fig4a.m. 
clearvars -except fig1

% C5 contains five colours in rgb format for plotting.  
C5 = [0.346666666666667,0.536000000000000,0.690666666666667;...
    0.915294117647059,0.281568627450980,0.287843137254902;...
    0.441568627450980,0.749019607843137,0.432156862745098;...
    0.640000000000000,0.640000000000000,0.640000000000000;...
    0.636470588235294,0.375294117647059,0.675294117647059];

% some axes positions (in cm). 
left_a = 0.775;
left_b = 4.775;
bottom = 0.825;
width = 3.15;
height= 3.2;

% array containing pump powers (measured at output of the fibre) [mW]. 
power = [20:20:200];


% loads the data file for the input spectra 
t = csvread('A0014.CSV',318,0);
% input wavelengths in column 1
input_wl = t(:,1);
% input intensity levels in column 2. 
input_lvl = t(:,2);

% The OSA will arbitarily populate some bins with an intensity value of 
% -210, so we remove datapoints corresponding to these values here. 
input_lvl(input_wl == 0) = [];
input_wl(input_wl == 0) = [];
input_wl(input_lvl == -210) = [];
input_lvl(input_lvl == -210) = [];

% remove data values far from input wavelength 
input_lvl(input_wl < 1090) =[];
input_wl(input_wl < 1090) =[];

input_lvl(input_wl > 1094) =[];
input_wl(input_wl > 1094) =[];

% the data file gives intensity in dBm, so when summing it is necessary to
% first convert to power on a linear scale. 
P_int_input = sum(10.^(0.1.*input_lvl));

% load noise data from csv file. 
n = csvread('A0011.CSV',318,0);
% noise wavelengths in column 1
n2wl = n(:,trace(1));
% noise intensity levels in column 2.
n2 = n(:,trace(2));


% The OSA will arbitarily populate some bins with an intensity value of 
% -210, so we remove datapoints corresponding to these values here. 
n2(n2wl == 0) = [];
n2wl(n2wl == 0) = [];
n2wl(n2 == -210) = [];
n2(n2 == -210) = [];

% remove data values far from output wavelength 
n2(n2wl < 1520) =[];
n2wl(n2wl < 1520) =[];
n2(n2wl > 1540) =[];
n2wl(n2wl > 1540) =[];

% the data file gives intensity in dBm, so when summing it is necessary to
% first convert to power on a linear scale. 
Pn = sum(10.^(0.1.*n2));

% This loop steps through the pump powers and opens the data files
% for each, finding the integrated power in each spectra. 
int_power=zeros(1,1);
for k = 1:10
    arry = [9 8 7 6 5 4 3 2 1 0];
    if arry(k) < 10
        fileroot = ['A000' num2str(arry(k)) '.CSV'];
    else
        fileroot = ['A00' num2str(arry(k)) '.CSV'];
    end 
    
    a = csvread(fileroot,318,0);
    % wavelengths in column 1 of data table
    a2wl = a(:,trace(1));
    % intensity levels in column 2 of the data table
    a2 = a(:,trace(2));

    % The OSA will arbitarily populate some bins with an intensity value of 
    % -210, so we remove datapoints corresponding to these values here. 
    a2(a2wl == 0) = [];
    a2wl(a2wl == 0) = [];
    a2wl(a2 == -210) = [];
    a2(a2 == -210) = [];
    
    
    % remove data values far from output wavelength 
    a2(a2wl < 1520) =[];
    a2wl(a2wl < 1520) =[];
    a2(a2wl > 1540) =[];
    a2wl(a2wl > 1540) =[];
    
    % populates int_power as a vector containing the intergrated power
    % around 1092 nm for each input wavelength. The data file gives 
    % intensity in dBm, so when summing it is necessary to first
    % convert to power on a linear scale. 
    int_power(k) = sum(10.^(0.1.*a2));
end

power = [0 power];
int_power = [Pn int_power];

Rp = 80*1e6; % Repetition rate of the laser [Hz] 
Tp = 10e-12; % pulse duration [s]
duty_cycle = Rp * Tp; 
ratio = 1530/1092; % equivalent to omega_Sr/omega_T in Eq.2

% --- calculation of up-conversion efficiency. 

eff = ratio*100.*(int_power-Pn)./(P_int_input*duty_cycle);

% --- Propagation of uncertainties ---
d_Tp = 0.5e-12; % uncertainty in pulse duration [s]
d_Pt = 0.05.*P_int_input; % uncertainty in input power 

% propagating uncertainties in pulse duration and input power. 
Err_Tp = (int_power./Rp).*(-1./(Tp^2 * P_int_input)).*d_Tp;
Err_Pt = (int_power./Rp).*(-1./(Tp * P_int_input^2)).*d_Tp;

% uncertainty in efficiency
Err_eff = ratio*100.*sqrt(Err_Tp.^2 + Err_Pt.^2);

% create axes for fig.4(c). 
ax2 = axes('units','centimeters','position',[left_b bottom width height]); 
hold on
% plot power vs conversion efficiency. 
p1 = plot(power,eff,'marker','v','markersize',2.5,'markerfacecolor',C5(1,:),'markeredgecolor','none','color',C5(1,:),'linestyle','none');
% plot error bars
errorbar(power,eff,Err_eff,'color',C5(1,:),'LineStyle','none','LineWidth',0.75,'capsize',2.5);

% commands to change appearance of figure.
box on;
xlabel('Power [mW]');
ylabel('Efficiency [%]');
ax2.LineWidth = 0.75;
ax2.FontSize = 8;
ax2.FontName = 'CMU Serif';
ax2.XLim = [0 205];
ax2.YLim = [0 5];
ax2.XTick = 0:40:200;
ax2.YTick = 0:1:5;
title('(c)');
