no longer needed
This commit is contained in:
sakul-45
parent
65f5a71e58
commit
05e1ede24f
217
EPR_script.m
217
EPR_script.m
@ -1,217 +0,0 @@
|
|||||||
% 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 = './';
|
|
||||||
% specifiy excel-file for saving parameters
|
|
||||||
table_path = 'example_results.xlsx';
|
|
||||||
|
|
||||||
%% 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 = string(extractBefore(extractAfter(path,asManyOfPattern(wildcardPattern + "/")),'.'));
|
|
||||||
% using Parallel Computing toolbox to speed up (check with "gpuDevice" if you can use this)
|
|
||||||
gpuData = gpuArray(Data);
|
|
||||||
|
|
||||||
%% 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
|
|
||||||
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
|
|
||||||
|
|
||||||
%% saving parameters
|
|
||||||
% concatenate data to existing table
|
|
||||||
table_old = readtable(table_path);
|
|
||||||
table_old.Properties.VariableNames = {'filename', 'date', 'pre-trigger', ...
|
|
||||||
'baseline_points', 'max_area_left', 'max_area_right', 'T_x', 'T_y', 'T_z', ...
|
|
||||||
'sim. g-value', 'sim_D', 'sim_E', 'sim_lw_gauss', 'sim_lw_lorentz'};
|
|
||||||
% new data as table
|
|
||||||
table_new = table(dataname, string(datestr(clock)), pre_trigger, baseline_points, ...
|
|
||||||
max_region(1), max_region(end), Exp.Temperature(1), Exp.Temperature(2), Exp.Temperature(3), ...
|
|
||||||
BestSys.g, BestSys.D(1), BestSys.D(2), BestSys.lw(1), BestSys.lw(2), ...
|
|
||||||
'VariableNames', {'filename', 'date', 'pre-trigger', 'baseline_points', ...
|
|
||||||
'max_area_left', 'max_area_right', 'T_x', 'T_y', 'T_z', 'sim. g-value', ...
|
|
||||||
'sim_D', 'sim_E', 'sim_lw_gauss', 'sim_lw_lorentz'});
|
|
||||||
table_conc = [table_old;table_new];
|
|
||||||
writetable(table_conc,table_path)
|
|
||||||
@ -1,7 +0,0 @@
|
|||||||
function [outputArg1,outputArg2] = double_2D_3D_plot(inputArg1,inputArg2)
|
|
||||||
%DOUBLE_2D_3D_PLOT Summary of this function goes here
|
|
||||||
% Detailed explanation goes here
|
|
||||||
outputArg1 = inputArg1;
|
|
||||||
outputArg2 = inputArg2;
|
|
||||||
end
|
|
||||||
|
|
||||||
Binary file not shown.
Loading…
x
Reference in New Issue
Block a user