%% FIGURE 4 A 
clear;
clc;

% Load Data
DataDir = 'Data_Raw\';
FileName = 'shfsg_20220503_002.mat';
data = load(strcat(DataDir, FileName));
XData{1} = 1e-9 .* data.smSHFSG_savePULSED.PULSED.Data.xAxis; % time data in s
YData{1} = data.smSHFSG_savePULSED.PULSED.Data.Counts(1,:); % count rate in CPS
rawdata = [ XData{1}; YData{1} ];

% Remove Background Drift
ConstBackGround = mean(YData{1});
f = fit(XData{1}', YData{1}', 'poly9');
toSubstract = f(XData{1}');
YDataCorr{1} = (YData{1}' - toSubstract + ConstBackGround)';
YDataCorrNorm{1} = YDataCorr{1} / max(YDataCorr{1}(10:end));

% Fitting Oscillation
startparam = [0.03, 240e3, 3*pi/4, 0.982, 0.000250, 10];
func = @(param, x) ( param(1)*cos(2*pi*param(2)*x+param(3)).*exp(-(x./param(5)).^param(6)) + param(4) );
[beta, R, J, CovB] = nlinfit(XData{1}, YDataCorrNorm{1}, func, startparam);
fitX = linspace(XData{1}(1), XData{1}(end), 10000);
fitY = func(beta, fitX);
fitValues = [beta' sqrt(diag(CovB))];

% Fourier Signal
DataStruct.XData{1} = XData{1};
DataStruct.YData{1} = YDataCorrNorm{1};
DataFourier = FourierAnalysisSHFSG(DataStruct, [200e3, 280e3], 0, [], 0, 1, 0);

% Save Results
dataX       = XData{1};
dataY       = YDataCorrNorm{1};
fourierX    = 1e-3 * DataFourier.XData{1};
fourierY    = abs(DataFourier.YData{1});%/max(abs(DataFourier.YData{1}));

save('Data_Processed_Ready_To_Plot\FIG_4_A.mat', 'dataX', 'dataY', 'fitX', 'fitY', 'fourierX', 'fourierY', 'fitValues');

%% FIGURE 4 B 
clear;
clc;

% Load Data
DataDir = 'Data_Raw\';
FileName = 'shfsg_20220530_009.mat';
data = load(strcat(DataDir, FileName));
XData{1} = 1e-9 .* data.smSHFSG_savePULSED.PULSED.Data.xAxis; % time data in s
YData{1} = data.smSHFSG_savePULSED.PULSED.Data.Counts(1,:); % count rate in CPS
rawdata = [ XData{1}; YData{1} ];

% Remove Background Drift
ConstBackGround = mean(YData{1});
f = fit(XData{1}', YData{1}', 'poly9');
toSubstract = f(XData{1}');
YDataCorr{1} = (YData{1}' - toSubstract + ConstBackGround)';
YDataCorrNorm{1} = YDataCorr{1} / max(YDataCorr{1}(10:end));

% Fitting Oscillation
startparam = [0.03, 240e3, 3*pi/4, 0.982, 0.000250];
func = @(param, x) ( param(1)*cos(2*pi*param(2)*x+param(3)).*exp(-x./param(5)) + param(4) );
[beta, R, J, CovB] = nlinfit(XData{1}, YDataCorrNorm{1}, func, startparam);
fitX = linspace(XData{1}(1), XData{1}(end), 10000);
fitY = func(beta, fitX);
fitValues = [beta' sqrt(diag(CovB))];

% Fourier Signal
DataStruct.XData{1} = XData{1};
DataStruct.YData{1} = YDataCorrNorm{1};
DataFourier = FourierAnalysisSHFSG(DataStruct, [200e3, 280e3], 0, [], 0, 1, 0);

% Save Results
dataX       = XData{1};
dataY       = YDataCorrNorm{1};
fourierX    = 1e-3 * DataFourier.XData{1};
fourierY    = abs(DataFourier.YData{1});%/max(abs(DataFourier.YData{1}));

save('Data_Processed_Ready_To_Plot\FIG_4_B.mat', 'dataX', 'dataY', 'fitX', 'fitY', 'fourierX', 'fourierY', 'fitValues');