From 4a83b50ca0a9b9746bcbc3bf772cf8c65aa79d51 Mon Sep 17 00:00:00 2001
From: sakul-45 <81963567+sakul-45@users.noreply.github.com>
Date: Sun, 25 Apr 2021 00:06:09 +0200
Subject: [PATCH] Update EPR_script.m

Added:
 - parallel computing to baseline correction
 - main simultation loop with decisions at its end
 - printing decision for fig(3)
 - extracting name of dataset for further use (printing figures, fir values in excel, ...)
---
 EPR_script.m | 85 ++++++++++++++++++++++++++++++++--------------------
 1 file changed, 53 insertions(+), 32 deletions(-)

diff --git a/EPR_script.m b/EPR_script.m
index ea76a8a..5b7c4f9 100644
--- a/EPR_script.m
+++ b/EPR_script.m
@@ -5,11 +5,18 @@ 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');
@@ -18,15 +25,20 @@ 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(Data)
+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, Data, mean(Data(1:pre_trigger,:)));
+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)
@@ -34,9 +46,8 @@ set(gcf,'Position',position)
 
 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')
+plot(signal_baseline_time.')
 title('transposed time corrected data')
 set(gcf,'Position',position)
 % figure(gcf) % bring figure to foreground
@@ -51,7 +62,7 @@ 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')
+plot(signal_baseline_time_field.')
 title('transposed fully corrected data')
 set(gcf,'Position',position)
 % figure(gcf) % bring figure to foreground
@@ -128,7 +139,7 @@ while init_proceed == 'n'
     
     figure(3)
     set (gcf,'PaperUnits','centimeters')
-    set (gcf,'Position',position) % set the position, size and shape of the plot
+    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);
@@ -138,41 +149,51 @@ while init_proceed == 'n'
     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
 
-return
  % 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
-% esfit('pepper',signal_baseline_field_time_mean_norm,Sys,Vary,Exp,[],FitOpt);%fitting route
-% return
+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
 
-[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',position) % 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')
-
-% set(gcf,'Units','Inches');
-% pos = get(gcf,'Position');
-% set(gcf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
-% print(gcf,'..\Abbildungen\Regression5','-dpdf','-r0');
\ No newline at end of file
+%% 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