tr-epr-simulation/EPR_script.m

% high-spin S=1 simulation
% Inputs requested in command line at certain points

clear variables
close all

 % window positions (currently optimised for dual WQHD with main on right)
 % also uncomment all figure(gcf) when working with single monitor
% Give desired position of figure window as number of pixels [pos_x pos_y size_x size_y]:
position = [-1250,50,1200,800];
% Give desired position of figure(2) window (will have two stacked subplots)
% as number of pixels [pos_x pos_y size_x size_y]: 
position2 = [-2000,50,700,800];
% Give desired position of figure(3) window (will contain EPR spectrum and 
% simulation) as number of pixels [pos_x pos_y size_x size_y]: 
position3 = [-1250,50,1200,600];

% specify dir for printing figures
figdir = './';

%% loading Data
path = input('Path to dataset: ','s');
load(path)

whos % what variables have been loaded
params % what information is contained in the structure called 'params'

% get name of dataset
dataname = extractBefore(extractAfter(path,asManyOfPattern(wildcardPattern + "/")),'.');

gpuData = gpuArray(Data); % using Parallel Computing toolbox to speed up

%% Baseline Correcting
% plot the raw data & check the number of points before the signal (pre-trigger)
plot(gpuData)
title('raw data')
set(gcf,'Position',position)

% substract the pre-trigger
pre_trigger = input('Number of pre-trigger points: ');
signal_baseline_time = bsxfun(@minus, gpuData, mean(gpuData(1:pre_trigger,:)));
plot(signal_baseline_time) % plot the corrected data set
title('time corrected data')
set(gcf,'Position',position)
% figure(gcf) % bring figure to foreground

ready = input('Proceed?');

 % plot the transpose and check the number of points to the lower and higher fields of the signal
plot(signal_baseline_time.')
title('transposed time corrected data')
set(gcf,'Position',position)
% figure(gcf) % bring figure to foreground


 % BASELINE correction
baseline_points = input('Number of baseline points (use smaller value from left and right): ');
l1 = mean(signal_baseline_time(:,1:baseline_points),2); % calculate the mean on the left along the time axis
l2 = mean(signal_baseline_time(:,end-baseline_points:end),2); %calculate the mean on the right along the time axis
baseline_time = (l1 +l2)/2; %take the average

signal_baseline_time_field = bsxfun(@minus, signal_baseline_time, baseline_time); % subtract the background in the time-domain 

 % plot the corrected data set
plot(signal_baseline_time_field.')
title('transposed fully corrected data')
set(gcf,'Position',position)
% figure(gcf) % bring figure to foreground

clear ready
ready = input('Proceed?');

 % plot the transpose to find the region of maximum signal. Use this below
plot(signal_baseline_time_field)
title('fully corrected data')
set(gcf,'Position',position)
% figure(gcf) % bring figure to foreground

 % contour plot: The index gives the number of contours
% contourf(signal_baseline_field_time,6)

 % NORMALISING
max_region = input('Region of the maximum signal as [x1:x2]: ');
% take the mean over the maxium region. You can decide how wide it is
signal_baseline_time_field_mean = (mean(signal_baseline_time_field(max_region,:)));
% normalise the amplitude to 1
signal_baseline_time_field_mean_norm = signal_baseline_time_field_mean/max(signal_baseline_time_field_mean);

%% Creating figure with two subplots
figure(2)
set(gcf,'PaperUnits','centimeters')
set(gcf,'Position',position2)
set(gcf,'InvertHardcopy','off','Color',[1 1 1])
set(0,'DefaultAxesFontSize', 12,'DefaultAxesLineWidth',2)

cont_or_surf = input('Should lower subplot be contour(1) or surface(2) plot? (1/2): ');
subplot(2,1,2)
if cont_or_surf == 1
    % contour plot: add the time and field axes
    contourf(0.1*params.Field_Vector, TimeBase*1e6 ,signal_baseline_time_field,'LineColor','none')
elseif cont_or_surf == 2
    % surface plot: add the time and field axes
    surf(0.1*params.Field_Vector, TimeBase*1e6 ,signal_baseline_time_field)
    colormap default
    shading interp
end

xlabel('Magnetic Field / mT')
ylabel('Time / \mus')

subplot(2,1,1)
 % plot the spectrum
plot(0.1*params.Field_Vector,signal_baseline_time_field_mean_norm,'LineWidth',2)

xlabel('Magnetic Field / mT')
axis('tight')
box off

%% Simulation section

Exp.mwFreq = params.mwFreq; % GHz
Exp.nPoints = length(params.Field_Vector);
Exp.CenterSweep = 0.1*[params.Field_Center params.Field_Sweep]; % mT (converted from Gauss)
Exp.Harmonic = 0; % zeroth harmonic

init_proceed = 'n';
while init_proceed == 'n'
    % populations of the triplet sub-levels
    % these need to be varied manually to get the right shape
    Exp.Temperature = input('Input population of triplett sublevels as [T_x T_y T_z]: ');
    % initial simulation settings
    Sys.S = 1; % Total Spin
    Sys.g = input('g value: '); % needs to be optimised
    Sys.D = input('D and E value as [D E]: '); % mT, The D and E values need to be optimised
    Sys.lw = input('Isotropic line broadening at FWHM as [Gaussian Lorentzian]: '); % mT, linewidth needs to be optimised
    
    [bfield,spec] = pepper(Sys,Exp); % perform a simulation with the parameters above
    spec_norm = spec/max(spec); % normalize the simulation
    
    figure(3)
    set (gcf,'PaperUnits','centimeters')
    set (gcf,'Position',position3) % set the position, size and shape of the plot
    set (gcf,'InvertHardcopy','off','Color',[1 1 1])
    set(0,'DefaultAxesFontSize', 16,'DefaultAxesLineWidth',1.5)
    plot(0.1*params.Field_Vector,signal_baseline_time_field_mean_norm,'r', bfield,spec_norm,'b','LineWidth',1);
    axis('tight')
    legend('experimental','simulation')
    legend boxoff
    xlabel('Magnetic Field / mT')
    ylabel('EPR signal / A. U.')
    set(gca,'Box','Off', 'XMinorTick','On', 'YMinorTick','On', 'TickDir','Out', 'YColor','k')
    pause(2);
    init_proceed = input('Spectrum shape manually fitted? [y/n]: ','s');
end

 % variation settings for simulation
Vary.g = 0.01; 
Vary.D = [10 10];
Vary.lw = [1 0];
 % further setup
FitOpt.Method = 'simplex fcn';
FitOpt.Scaling = 'lsq';

 % When you have got a good fit by eye, use esfit to optimise
simu_proceed = 'n';
while simu_proceed == 'n'
    % fitting routine
    [BestSys,BestSpc] = esfit('pepper',signal_baseline_time_field_mean_norm,Sys,Vary,Exp,[],FitOpt);
    % plot best fit
    figure(3)
    plot(0.1*params.Field_Vector,signal_baseline_time_field_mean_norm,'r',...
         0.1*params.Field_Vector,BestSpc,'b','LineWidth',1);
    axis('tight')
    legend('experimental','simulation')
    legend boxoff
    xlabel('Magnetic Field / mT')
    ylabel('EPR signal / A. U.')
    set(gca,'Box','Off', 'XMinorTick','On', 'YMinorTick','On', 'TickDir','Out', 'YColor','k')
    
    simu_proceed = input('Did the simulation converge? [y/n]: ','s');
    if simu_proceed == 'n'
        simu_val = input('Do you want to repeat the simulation with new best values? [y/n]: ','s');
        if simu_val == 'y'
            Sys.g = BestSys.g;
            Sys.D = BestSys.D;
            Sys.lw = BestSys.lw;
        end
    end
end

%% printing figures and saving parameters
printing = input('Do you want to print figure(3)? [y/n]: ','s');
if printing == 'y'
    figure(3)
    set(gcf,'Units','Inches');
    pos = get(gcf,'Position');
    set(gcf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
    print(gcf,strcat(figdir,dataname),'-dpdf','-r0');
end