% ------info ---------------------------------------
% Routine provides fit of raw ETR vs. E data 
% The format of the file needs to be as follows:
%  data on sheet1, row 1 = header; 
%  columns should alternate PAR, ETR for each curve
% ---- END info -------------------------------------

clear all;

w=cd;
fprintf('Specify the file containing PAM data \n');
[fname, file_path]=uigetfile('*.*');
cd (file_path)
fid         = fopen(fname);
[A, sheet]  = xlsfinfo(fname);

type = input('\n Specify your fit preference \n 1: rectangular hyperbola (red) \n 2: neg. exponential (green) \n 3: hyperbolic tan (blue) \n 4: tanh with photoinbition (black) \n 5: All \n:');
% MODEL 1. rect. hyp =  (Pmax*E)/ (E + (Pmax/alpha))
% MODEL 2. neg exp   =  Pmax(1- exp(-E*(alpha/Pmax));
% MODEL 3. hyp. tan  =  Pmax*tanh((E*alpha)/Pmax);
% MODEL 4. photoinh  =  Pmax*tanh((E*alpha)/Pmax)* 1/(1+ beta*E);

% read .xls file: format should be par-etr par-etr etc for successive samples
    [data, header]  = xlsread(fname,char(sheet(1)));
    status=fclose(fid);
% strip header row leaving numeric data
    ndata = data (1:end,:);

% create stucture for PAM light curve data
odd    = 1:2:length(ndata)-1;
even   = 2:2:length(ndata);
pam.id = header(1,odd);
pam.par= ndata(:,odd);
pam.etr= ndata(:,even);

if (type==5),
    loop = type-1;
else
    loop = 1;
end

for iterations = 1:loop,
    if (type==5 | iterations>1),
        type=iterations;
    end

    %fit curves
    for i=1:1:(length(ndata)/2),
        
       % o2(:,i)=.25*0.5*abs(i)*(pam.etr(:,i)./(.84*.5));

        x=pam.par(:,i);
        y=pam.etr(:,i);

        %provide initial guesses for alpha, beta and pmax;
        beta(1) = max(y);
        temp    = polyfit(x(1:5),y(1:5),1);
        beta(2) = temp(1);
        beta(3) = 0.0;
        if (type==1)
            fun = inline ('(beta(1).*x)./(x+(beta(1)./beta(2)))','beta','x');
            color= 'r-';
        elseif (type==2)
            fun = inline ('beta(1)*(1-exp((-x.*beta(2))./beta(1)))','beta','x');
            color= 'g-';
        elseif (type==3)
            fun = inline ('beta(1).* tanh((x.*beta(2))./beta(1))', 'beta', 'x');
            color = 'b-';
        elseif (type==4)
            fun = inline ('beta(1).*(tanh((x.*beta(2))./beta(1)).*(1./(1+(beta(3).*x))))', 'beta', 'x');
            color= 'k-';
        end
        %estimate model parameters with nlinfit
        [param, r, J] = nlinfit(x, y, fun, beta);
        ypred = nlpredci(fun,x,param,r,J);
        sse = norm(r).^2;
        sst = norm(y - mean(y)).^2;
        rsquare = 1 - sse/sst; 

        model(i).pmax  = param(1);
        model(i).alpha = param(2);
        model(i).ik    = param(1)/param(2);
        model(i).beta  = param(3);
        model(i).ypred = ypred;
        model(i).rsq   = rsquare;

% PLOTTING--- up to 23 profiles---------------------------------------
        figure (1);

% Change the subplot total to fit your data!
        subplot (2,3,i);
     
        plot (pam.par(:,i),pam.etr(:,i),'k.'); hold on;
        plot (pam.par(:,i),model(i).ypred,color);
        axis([0 1200 0 50]);
        title (char(pam.id(i))); xlim ([0 1200]);
    end
    
    clear x y r J color param ypred fun sse sst rsquare
end

cd (w);
 stats(:,1)=[model.pmax]';
 stats(:,2)=[model.alpha]';
 stats(:,3)=[model.ik]';
 stats(:,4)=[model.beta]';
 stats(:,5)= type;
 stats(:,6)=[model.rsq]';
 

clear alpha b s file_path fid fname wd i  j s sampleid status param temp loop iterations type