% T.A.Wright et al (2020) 

% This script plots the wavelength-dependant conversion efficiency, shown
% in Fig.4(a) (inset).
% 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 
% 0001 and 00025. 

% figure object fig1 is shared with plot_fig4a_inset.m, plot_fig4b.m 
% and plot_fig4c.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];


% axes positions 
left_inset = 4.3;
bottom_inset = 1.75+4.5-0.1;
width_inset = 2;
height_inset = 1.75;

% Calculates input wavelength in nm from frequency of input laser. 
input_wavelength = 1e9.*3e8./((195.5:0.025:196.1).*1e12);

% read in noise spectra
n = csvread('wl/A0024.CSV',318,0); 
% noise wavelengths in column 3
n2wl = n(:,3);
% noise intensity levels in column 4. 
n2 = n(:,4);

% 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,length(input_wavelength));

% input spectra at telecom wavelength, recorded for purposes of measuring
% power. We show a normalised efficiency so this isn't actually needed. 
t = csvread('wl/A0025.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 input wavelengths and opens the data files
% for each, finding the integrated power in each spectra. 
for k = 1:length(input_wavelength)
    
    file_num = 0:24;
    
    if file_num(k) < 10
        fileroot = ['wl/A000' num2str(file_num(k)) '.CSV'];
    else
        fileroot = ['wl/A00' num2str(file_num(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) = [];
    
    % ppulates int_power as a vector containing the intergrated power
    % around 1092 nm for each input wavelength. 
    int_power(k) = sum(10.^(0.1.*a2));
end

Rp = 80e6; % Repetition rate of the laser [Hz]
Tp = 9.5e-12; % Pump-pulse duration [s]
duty_cycle = Rp * Tp; 

% Calculates efficiency according to Eq.1 in the manuscript. 
% We neglect the omegas here as we normalise to the maximum value anyway. 
eff = 100.*int_power./(Pt*duty_cycle);  

% axes object for the plot inset in Fig.4(a). 
ax1 = axes('units','centimeters','position',[left_inset bottom_inset width_inset height_inset]); 

hold on
% Plots normalised conversion efficiency as a function of input wavelength.
plot(input_wavelength,eff./max(eff),'marker','v','markersize',2.5,'markerfacecolor',C5(3,:),'markeredgecolor','none','color',C5(1,:),'linestyle','none');

% Plot commands to add labels and change the appearance of the figure. 
ax1.LineWidth = 0.75;
ax1.FontSize = 7;
ax1.FontName = 'CMU Serif';
ax1_pos = ax1.Position; 
ax2 = axes('units','centimeters','Position',ax1_pos,...
    'XAxisLocation','top',...
    'YAxisLocation','right',...
    'Color','none');
ax2.LineWidth = 0.75;
ax2.FontSize = 7;
ax2.FontName = 'CMU Serif';
ax2.XTick = 1091:2:1093;
set(ax2,'YAxisLocation','right');
ax2.YLim = [0 1.0];
ax1.YLim = [0 1.0];
ax2.YTick = [];
ax2.YTickLabel = [];
ax1.YTickLabel = [];
set(ax2,'XAxisLocation','top');
set(ax2,'TickLength',[0.025, 0.025])
set(ax1,'TickLength',[0.025, 0.025])
ax1.XTick = [1529:3:1535];
ax1.XLim = [1528 1535];
ax2.XLim = -1./(1/977.2 - 1/777.5 - 1./([1528 1535])); 