%metzyme_rates.m
%22 nov 2013
%asantoro@umces.edu

%additional comments and clarifications added for BCO-DMO submission, 25
%mar 2016

%calculates the nitrification rates for three stations in the equatorial
%pacific during a cruise in october 2011. 15NH4 was added to to duplicate bottles,
%plus a no addition control.

%this script must be run separately for each station. the input files are
%hard-coded and must be edited manual as indicated below between lines marked with '*****'.

%this script calls the function nitr_rate_rev.m and requires the MATLAB optimization
%toolbox (for lsqcurvefit.m).

%bottle numbers:
%1-15NH4+ (bottle A)
%2-15NH4+ (bottle B)
%3-no addition (bottle C)
%=====================================

clear all

%load data on initial conditions, columns:
%1-depth (m)
%2-AP15NH4, starting atom fraction 15N (15N/(15N+14N) calculated from mass balance of tracer addition and starting [NH4+]
%3-No_15 (NO3-), the initial concentration (uM) of 15N-NO3-
%4-No_14 (NO3-), the initial concentration (uM) of 14N-NO3-

%********************
%select input files here  (initials and isotope data) by uncommenting the
%appropriate files

initials = load('initials_mzst1.csv');
%initials = load('initials_mzst3.csv');
%initials = load('initials_mzst5.csv');


%load isotope data.  columns: 1-depth, 2-bottle# (see above), 3-time(h), 4-15R-NO3

data = load('data_mzst1.csv');
% data = load('data_mzst3.csv');
 %data = load('data_mzst5.csv');
%***********************

 
depnum = length(initials); %number of depths
botnum = 3; %number of treatments

%set curve fitting parameters
x0 = [0.00001 0.0001]; %initial guess for k and F 
lb = 0;  %sets lower bound of zero, i.e. k and F are positive
ub = []; %upper bound set to one, change to empty matrix if there is a fitting error


%set up figure for data and fitting results
figure

%setup matrix for data
metzymerates = zeros(depnum,5);


%break data into matrices by depth and bottle, then fit rate function

for i=1:depnum  %for each depth
    
    dep = initials (i,1); %temporary depth
    eval(['data_' num2str(dep) '= data (find (data(:,1) == dep),:);']) %break into matrix by depth
    tempdata = eval(['data_' num2str(dep)]); %set tempdata to current depth
    c = initials (i,2:4); %specify coefficients for rate fitting 
    
        
    for b=1:botnum %for each bottle
        
        bot = tempdata (find (tempdata(:,2) == b),:); %find the data for bottle = botnum
        time = bot(:,3); %time data for fitting
        ydata = bot(:,4); %15R data for fitting
                
        %use least squares curve fitting to calculate a, the vector of unknown
        %coeffs.  a(1) = k, a(2) = F.  see nitr_rate.m
        [a,resnorm,r,exitflag,output,lambda,J] = lsqcurvefit(@(a,time) nitr_rate_rev(a,time,c), x0,time,ydata,lb,ub);
        fit = nitr_rate_rev(a,time,c);
        
        subplot(depnum,botnum,b+(i-1)*botnum) %plot
        plot(time,ydata,'*') %real data
        hold on
        plot(time,fit,'r-') %fit 
        
        rt = a(2)*1000*24; %convert rate to nM per day
        k_removal = a(1) *24; % convert rate constant to per day
        
       
        % calculate standard error of fit coefficient
        % approximate the Hessian using the Jacobian, and use the inverse of that
        % as an approximation to the covariance matrix:
        % 
        [Q,R] = qr(J,0); %decomposition of the jacobian  
        mse = sum(abs(r).^2)/(size(J,1)-size(J,2));
        Rinv = inv(R);
        Sigma = Rinv*Rinv'*mse;
        SE_uM = sqrt(diag(Sigma)); %SE in uM per h
        SE_nM = SE_uM(2)*1000*24;     %SE in nM per day 
        % 
        % where r is the residual vector and J is the Jacobian, both outputs of
        % LSQCURVEFIT. 
        
        %build matrix
        metzymerates (b+(i-1)*botnum,1) = dep;   %first column of rate matrix is depth
        metzymerates (b+(i-1)*botnum,2) = b;     %second column is bottle #
        metzymerates (b+(i-1)*botnum,3) = rt;    %third column is rate in nM per day
        metzymerates (b+(i-1)*botnum,4) = SE_nM; %fourth column is SE of the rate in nM per day
        metzymerates (b+(i-1)*botnum,5) = a(1);  %fifth column is the removal rate constant in units of per hour
        
    end
    
end

%output graph
print -painters -dpng metzyme_rates_st1.png 

%output rates to csv file
dlmwrite('metzyme_rates_st1.txt', metzymerates, 'delimiter', '\t')