
clear
close all


% load MMS data
% trex = readtable('MMS_NaOH_stdAdd_20180530.txt');




Istray = 50; 

MMS.Dark = trex{:,6}';
% MMS.I_blank = trex{1,10:end}';
% MMS.I_blank_D = MMS.I_blank - MMS.Dark(1); 
MMS.I = trex{:,10:end}';
MMS.I_D = MMS.I - MMS.Dark - Istray;

% use example output from 
trex = readtable('MSC_examp_reply_blank_samp.txt');





MMS.Araw = log10(MMS.I_D(:,1) ./ MMS.I_D(:, 2:end));


% Lamda calculation
% UVVIS1. Pretty sure this is the correct UVVIS 
C0 = +2.455646892E+02;
C1 = +2.136483472E+00;
C2 = +1.005266419E-04;
C3 = -7.432017061E-07;
C4 = +0.0;

X = [1:256]'; 
MMS.Lambda = C0 + C1.*X + C2.*X.^2 + C3.*X.^3;
% Baseline correction
ibase = inrange(MMS.Lambda, [750 760]);
MMS.Araw_blc = MMS.Araw - mean(MMS.Araw(ibase,:));

n = 2; % order polynomial
iuse = inrange(MMS.Lambda, [573 583]);
npixel = 6;


MMS.I_D_mono = NaN(size(MMS.I_D));
MMS.x = MMS.I_D_mono;




for i = npixel+1:length(MMS.Lambda)-npixel

    
    % get 5 points around pixel;    
    iuse = i-npixel:i+npixel;
    
    for j = 1:size(MMS.Araw,2)
    
        % extract I and lambda
        I = MMS.I_D(iuse,j);
%         I = MMS.Araw(iuse,j);
        wl = MMS.Lambda(iuse);
        
        c = polyfit(wl, I, n); % fit polynomial curve to intensity
        
        % X grid for calculations
        X = wl(1):0.01:wl(end);
        
%         % get minima assuming gaussian distribution around pixel lambda
        mu = wl(npixel+1);
        sigma = sqrt(4);
        ndist = get_normdist(X, mu, sigma);
        [f, S, Cnt] = LMFsolve(@(x) I(npixel+1) - trapz(X,polyval(c,X) .* ndist .* x), 1, 'default');
        
        % assume all wavelenghts contribute evenly to the average
%         [f, S, Cnt] = LMFsolve(@(x) I(npixel+1) - mean(polyval(c,X)).* x, 1, 'default');
        
        
        MMS.I_D_mono(i,j) = MMS.I_D(i,j) * f;
%         MMS.I_D_mono(i,j) = MMS.Araw(i,j) * f;
        MMS.x(i,j) = f;
        
        
    end
end

% calculate absorbance

MMS.I_blank = MMS.I_D_mono(:,1);
MMS.I_samp = MMS.I_D_mono(:,2:end);

MMS.A = log10(MMS.I_blank ./ MMS.I_samp);

% baseline correction
ibase = inrange(MMS.Lambda, [750 760]);
MMS.A_blc = MMS.A - mean(MMS.A(ibase,:));
% MMS.Araw_blc = MMS.Araw - mean(MMS.Araw(ibase,:));

figure
hold on; box on;
plot(MMS.Lambda, MMS.A_blc);
plot(MMS.Lambda, MMS.Araw_blc);


figure
plot(MMS.Lambda, MMS.A_blc - MMS.Araw_blc);
% plot(MMS.Lambda, MMS.I_D_mono(:,1:end-1) - MMS.Araw);
ylabel('A_{corrected} - A{raw}');
xlim([400 800]);


% interpolate onto 1 nm grid
X_ = 400:1:800;
MMS.A_blcint = NaN(length(X_), size(MMS.A_blc,2));
for i = 1:size(MMS.A_blc,2)
    
    MMS.A_blcint(:,i) = interp1(MMS.Lambda, MMS.A_blc(:,i), X_, 'spline');
    
    
end


% 
% I = MMS.I_D(iuse,1); % 156 is 578 nm, so 5 points centering around 578
% wl = MMS.Lambda(iuse);
% 
% c = polyfit(wl, I, n); %fit polynomial curve
% 
% X = wl(1):0.01:wl(end);
% mu =wl(3);
% sigma = sqrt(3);
% 
% ndist = get_normdist(X, mu, sigma);
% % 
% % I_ = polyval(c,X) .* ndist .* f;
% [x_, S, Cnt] = LMFsolve(@(x) I(3) - trapz(X,polyval(c,X) .* ndist .* x), 1, 'default');
% 
% 
% I_real = I(3).*x_;


% [x, S, Cnt] = LMFsolve(@fun2, x0, 'default'); % use matalb exchange function














