%etsp_rates_2011_indvT0.m
%10 july 2014
%asantoro@ucsb.edu

%calculates the nitrification rates for six stations in the etsp in 2011.  

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

%modified 14 apr 2011 to allow for input of individual T0 values when there
%is a lot of variability between T0 values.

%modified 21 aug 2013 (has it really been 3 years?!?!) to calculate NO2
%uptake rate from k and starting [NO2-]

%removed above addition 10 july 2014 for etsp data.  added some code so
%that the script can figure out the depths and they don't have to be
%specified.

%input and output file names need to be adjusted for each station in the following lines: 27, 31, 112, 115.

clear all

%load data on initial conditions: 1-depth, 2-AP15NH4 (or NO2), 3-No_15 (NO3-),
%4-No_14 (NO3-)

initials = load('initials_st11b_NO3add.csv');

%load isotope data.  columns: 1-depth, 2-bottle (=treatment, must start at 1 and be sequential), 3-time(h), 4-15R-NO3

data = load('data_st11b_NO3add.csv');


botnum = max(data(:,2)); %number of bottles/treatments

depths = unique(initials(:,1)); %vector of the depths

depnum =  length(depths);%number of depths

%set curve fitting parameters
x0 = [0.0001 0.001]; %initial guess for k and F 
lb = 0;  %sets lower bound of zero, i.e. k and F are positive
ub = [1]; %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
etsprates = zeros(depnum,4);

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

for i=1:depnum  %for each depth
    
    dep = depths (i); %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
    
    
        
    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
        c = initials (b+(i-1)*botnum,2:4); %specify coefficients for rate fitting
        
                
        %use least squares curve fitting to calculate a, the vector of unknown
        %coeffs.  a(1) = k, a(2) = Fin.  see also comments within function 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 are blue stars
        hold on
        plot(time,fit,'r-') %fit is red line
        
        rt = a(2)*1000*24; %convert rate to nM per day
        k_uptake = 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)); %standard deviation
        
        Rinv = inv(R); %inverse of the Hessian approximates the covariance matrix
        
        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 of rate data
        etsprates (b+(i-1)*botnum,1) = dep;   %first column of rate matrix is depth
        etsprates (b+(i-1)*botnum,2) = b;     %second column is bottle #
        etsprates (b+(i-1)*botnum,3) = rt;    %third column is rate in nM per day
        etsprates (b+(i-1)*botnum,4) = SE_nM;    %fourth column is SE in nM per day        
        
    end
    
end


%output graph
print -painters -dpng etsp_st11b_NO3add.png 

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





