% T.A.Wright et al (2020) 

% This script plots the pump-power-dependant conversion efficiency, shown
% in Fig.4(b).
% 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_Sr - N_sr) * omega_T)/ (P_T * D * omega_Sr). 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 00012. 

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

% 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 = 15:15:180;

% load noise data from csv file. 
n = csvread('pwr/A0000.CSV',318,0);
% noise wavelengths in column 1
n2wl = n(:,1);
% noise intensity levels in column 2. 
n2 = n(:,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) = [];

background = mean(n2);

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

% array for intergrated power of converted signal for each input wavelength
int_power=zeros(1,12);

% loads the data file for the input spectra at telecom wavelength,
t = csvread('pwr/A0012.CSV',318,0);
% Telecom wavelengths in column 1
t2wl = t(:,1);
% Telecom intensity levels in column 2. 
t2 = t(:,2);

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

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


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

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

% --- calculation of up-conversion efficiency. 

eff = ratio*100.*(int_power - background)./(Pt*duty_cycle); % Eq.1

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

% propagating uncertainties in pulse duration and telecom power. 
Err_Tp = (int_power./80e6).*(-1./(Tp^2 * Pt)).*d_Tp;
Err_Pt = (int_power./80e6).*(-1./(Tp * Pt^2)).*d_Pt;

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

% create axes for fig.4(b). 
ax3 = axes('units','centimeters','position',[left_a bottom width height]); 
hold on
% plot power vs conversion efficiency. 
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 [%]');
ax3.LineWidth = 0.75;
ax3.FontSize = 8;
ax3.FontName = 'CMU Serif';
ax3.XLim = [0 185];
set(ax3,'XMinorTick','on');
ax3.XAxis.MinorTickValues = 0:25:175;
set(ax3,'YMinorTick','on');
ax3.YAxis.MinorTickValues = 0:10:6;
title('(b)');


