% T.A.Wright et al (2020) 

% This script contains all of the commands used to load, process and plot 
% the data used in Fig.3. Fig.3(b) is plotted first, then Fig.3(a).

% In Fig.3(b) we plot distributions of pulse energy in FWM idler pulses for
% different seed powers. The raw data takes the form of .dat files, which 
% contain oscilliscope traces obtained when monitoring a silicon photodiode. 

% In Fig.3(a) we plot the seeded and unseeded FWM spectra obtained on an
% Optical Spectrum Analyser (OSA). These spectra were saved by the OSA
% using a csv format. 

close all
clear

%% Fig.3(b): Normalised Peak Voltage (Pulse Energy) histograms 

% Fig.3(b) data file input names: 

% In Fig.3(b) we plot distributions of pulse energy in FWM idler pulses 
% for different seed powers. Lines show the least-square Gaussian fits.

% To plot the normalised peak-voltage distributions it is necessary to load 
% the many oscillisope traces used to record the photodiode signal. These 
% were saved from the oscilliscope the with file names of the form 
% C3Trace00xxx.dat.
%
% C3Trace00236.dat to C3Trace00256.dat correspond to the data when the seed
% power was 00 mW. 

% C3Trace00258.dat to C3Trace00279.dat correspond to the data when the seed
% power was 05 mW. 

% C3Trace00281.dat to C3Trace00302.dat correspond to the data when the seed
% power was 10 mW. 

% C3Trace00304.dat to C3Trace00326.dat correspond to the data when the seed
% power was 20 mW. 

% C3Trace00330.dat to C3Trace00354.dat correspond to the data when the seed
% power was 30 mW. 

% These are stored in the variable file_nums2: 
file_nums2 = [236 256; 258 279; 281 302; 304 326; 330 354];

%% Fig.3(b): defining bins for histograms 

% When we produce histograms of normalised peak voltage we use the 
% histcounts function. [N,edges] = histcounts(X,defined_edges) sorts X into bins with the bin 
% edges specified by the vector, defined_edges.

% We normalise our peak voltages to the mean value. So the the
% defined_edges vector specified below means our bins run from 0 signal to
% twice the value of the mean peak height. 

defined_edges = 0:0.01:2;

%% Fig.3(b): Some variables for plotting 

% C5 contains five colours in rgb format. 
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];

% these variables contain positions for both axes within figure (in cm)
left_a = 0.775-0.35;
left_b = 4.775;
bottom = 0.825;
width = 3.15;
l_width = 3.5;
height= 3.2;

% create figure object 
fig3 = figure('color','white','units','centimeters','position',[5.5,9,8,4.5]);

% create an axes object for figure 3(b)
ax2 = axes('units','centimeters','position',[left_b bottom width height]);
hold on;

%% Fig.3(b): Importing oscilliscope traces and data analysis

% This outer loop just runs through repeating the same steps of importing, 
% processing, and plotting data for seed powers of 0, 5, 10, 20 and 30 mW. 
for seed = 1:length(file_nums2)

    % The cumulative selection of normalised peak heights are stored in 
    % the array, tot_pks. We reuse this for every seed power, so it is 
    % reset here. 
    tot_pks =[];
    
    %Loop for reading in oscilliscope traces from data files, finding peaks...
    % in those traces and then normalising those peaks to the mean value. 
    for index = file_nums2(seed,1):file_nums2(seed,2)
        if index < 10
            structroot1 = sprintf('C3Trace0000%d.dat',index);
        elseif index < 100 && index > 9 
            structroot1 = sprintf('C3Trace000%d.dat',index);
        else
            structroot1 = sprintf('C3Trace00%d.dat',index);
        end
        a = load(structroot1);

        mk = 0.2*max(a(:,2));
        [pks,locs] = findpeaks(a(:,2),'MinPeakProminence',mk);
        
        % find mean peak height 
        avg_height = mean(pks);
        % normalise current selection of peaks to the mean height * 
        nrm_peaks = pks./avg_height;
        % Cumulative selection of normalised peak heights are stored here
        tot_pks = [tot_pks; nrm_peaks];
        
        % * Note that normalising relative to the mean peak height for each
        % individual recorded trace is useful for washing out how the the 
        % FWM sideband might be effected by experimental parameters that 
        % might be changing over longer time scales
    end
    
    %[N,edges] = histcounts(X,defined_edges) sorts X into bins with the bin 
    % edges specified by the vector, defined_edges.
    [n, edges] = histcounts(tot_pks,defined_edges);
    edges = edges(2:end) - (edges(2)-edges(1))/2;
    % convert n in to a percentage, as oppose to a raw number
    n = 100.*n./sum(n);
    % apply a gaussian fit to the data 
    f = fit(edges',n','gauss1');
    % remove unoccupied bins (for the purposes of plotting)
    edges(n==0) =[];
    n(n==0) = [];
    % scatter plot of the percenatge of peaks in each bin, 
    % plotted on axes 2 (Fig.3(b)).
    sc1 = scatter(ax2,edges,n,0.5,'filled','v');
    % recolour the plot according to the colour plot defined by the
    % variable C5
    sc1.MarkerFaceColor = C5(seed,:);
    sc1.MarkerEdgeColor = C5(seed,:);
    % hold command to enable multiple things to be plotted on the same axes
    hold on 
    % Plot the fitted gaussian distribution 
    plo = plot(f);
    % change things such as line width and colour of the line
    plo.LineWidth = 0.75;
    plo.Color = C5(seed,:);
    set(get(get(sc1,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
    % std(tot_pks) if uncommented this prints standard deviation to the command window
end
%% Some commands to change the appearance of Fig.3(b)
% places a legend on axes 2 (Fig.3(b)).
leg = legendflex({'0 mW','5 mW','10 mW','20 mW','30 mW'},'xscale',0.3,'fontsize',7,'box','off','fontname','CMU Serif', 'anchor', [3 3], 'buffer', [-0.75 -5]);

title('(b)');    
ylab = ylabel('Pulses Per Bin [%]');
xlabel('Norm. Peak Voltage');
set(ax2,'FontSize',8);
set(ax2,'LineWidth',0.75);
set(ax2,'FontName','CMU Serif');
xlim([0.45 1.55]);
ylim([0 12]);
ax2.XTick = 0.5:0.25:1.5;
ax2.YTick = 0:3:12;
box on

set(ylab, 'position', [0.3651,6,-1]);  % shift the y label to the left by 0.1


%% Fig.3(a): The seeded and unseeded FWM spectra 

% clear some of the workspace used to plot the pulse-energy histograms. 
clearvars -except fig2 ax2 left_a bottom width height l_width

% define colours used for plotting. 
p_colour1 = [0.346666666666667,0.536000000000000,0.690666666666667];
p_colour2 = [0.915294117647059,0.281568627450980,0.287843137254902];

% define the offset in dB used to seperate the spectra in Fig.3(a)
sep = 90.*[0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12];

% The spectra we plot in Fig.3(a) are saved from the OSA in csv files. 
% This loop reads those .csv files and stores them in the workspace as 
% an array of structs

for index = 90:113 

    if index == 0
        fname = 'W0000.CSV';
    elseif index >= 1 && index < 10
        fname = strcat('W000',num2str(index),'.CSV');
    elseif index <100 && index >= 10
        fname = strcat('W00',num2str(index),'.CSV');
    else 
        fname = strcat('W0',num2str(index),'.CSV');
    end
    
    % nosa_read is a short function which reads the wavelength and
    % intensity columns out of the csv file. The syntax for selecting
    % wavelengths for the nth spectra is: wavelength = A(n).WavelengthOut. 
    % To select the intensity for the nth spectra, enter:
    % intensity = A(n).LevelOut. 
    A(index+1) = nosa_read(fname);
end

% As A has a lot of empty rows it is tidier to just define a new variable 
% to keep only the bits of A which actually contain any data. AP is this
% new variable. 
AP = A(91:114);

% create an axes for plotting Fig.3(a)
ax1 = axes('units','centimeters','position',[left_a bottom l_width height]);
hold on
box on

% This loop removes all data for wavelengths below 600 nm. 
for kcount = 1:24
    AP(kcount).LevelOut(AP(kcount).WavelengthOut<600) = [],
    AP(kcount).WavelengthOut(AP(kcount).WavelengthOut<600) = [],
end


% Plots the unseeded FWM spectra (seperated by the offsets defined 
% in the variable 'sep')
ax1p2 = plot(ax1,... 
    AP(6).WavelengthOut,AP(6).LevelOut+sep(1),...
    AP(10).WavelengthOut,AP(10).LevelOut+sep(2),...
    AP(14).WavelengthOut,AP(14).LevelOut+sep(3),...
    AP(18).WavelengthOut,AP(18).LevelOut+sep(4),...
    AP(22).WavelengthOut,AP(22).LevelOut+sep(5),...
    'LineWidth',0.75,'Color',p_colour2);

% Plots the seeded FWM spectra (seperated by the offsets defined 
% in the variable 'sep')
ax1p1 = plot(ax1,... 
    AP(5).WavelengthOut,AP(5).LevelOut+sep(1),...
    AP(9).WavelengthOut,AP(9).LevelOut+sep(2),...
    AP(13).WavelengthOut,AP(13).LevelOut+sep(3),...
    AP(17).WavelengthOut,AP(17).LevelOut+sep(4),...
    AP(21).WavelengthOut,AP(21).LevelOut+sep(5),...
    'LineWidth',0.75,'Color',p_colour1);

% This loop adds the labels indicating the pump power to the plot. 
sepcount = 0;
for kcount = 300:100:700
    sepcount = sepcount-1;
    powerstring = num2str(kcount);
    fullstring = [powerstring ' mW'];
    text(1110, 5 +sepcount*90,fullstring,'FontName','CMU Serif','FontSize',7);
end

% The following commands are used to change the appearance of Fig.3(a),
% and to add axis labels. 
title('(a)');
lim = ylim(ax1);    ax1.YLim = [lim(1) 0];
xlabel(ax1,'Wavelength [nm]');
ylab1 = ylabel(ax1,'Intensity [100 dB/div.]');
ax1.LineWidth = 0.75;
ax1.FontSize = 8;
ax1.FontName = 'CMU Serif';
ax1.XLim = [590 1325];
ax1.YLim = [-485 15];
ax1.YTickLabel = [];




