clear all
close all
clc


%load data from .csv
dataset=readtable('dataset.csv');

rowok=1:1:size(dataset,1);

%%to test only one plot uncomment the following lines
% site='GL';
% rowok=find(strcmp(site,dataset.site));

%%to test only one year uncomment the following lines
yr=2020;
rowok=table2array(dataset(:,4))==yr; 

data=table2array(dataset(rowok,2:end));
er=data(:,13);
ta=data(:,7);

%column number of candidate additional variables in 'data'
addedvar=[9,11]; 

%additional variables provided by the user: uncomment one of the follwoing 
% myi=[]; % empty array to build the model from scratch 
myi=[11,9,4]; % filled with column numbers to fit the model using preselected vatiables  

%open output file
fileID = fopen('ERmodel.txt','w+'); 

%% first step: ER vs Ta

fprintf(fileID,'first step: \n');
fprintf(fileID,'a_0 b_0 varres varexpl AIC h \n');

% exclude rows with NaN values
good=(~isnan(er) & ~isnan(ta));

% exponential fit 
[fit1,gof,output]=fit(ta(good),er(good),'exp1');
a1=[fit1.a, fit1.b];
res=output.residuals;
varexpl=var(a1(1).*exp(ta(good).*a1(2)))/var(er(good));
AIC1=length(er(good))*log(sum(res.^2))+2*length(a1);
h=lillietest(res);

% write fitted values and staistics on output file
fprintf(fileID,[repmat('%g ',1,5),'\n'],a1,varexpl,AIC1,h);

 %% second step: ER vs additional variable #1

fprintf(fileID,'\n');
fprintf(fileID,'second step \n');
fprintf(fileID,'variable b_0 a_0 a_1 pvals varexpl AIC h \n');

n=1; % flag for variable selection

% loop over candidate additional variables
if(isempty(myi))
    
for i=addedvar
    
    clear good1 X a2 f2 fit2
    
    good1=~isnan(data(:,i)) & ~isnan(er) & ~isnan(ta); 
    
    % compute the starting values for nonlinear fit
    er2=er(good1)./exp(a1(2).*ta(good1));
    f2=regress(er2,[ones(size(data(good1,i))) data(good1,i)]);
    a2g=[a1(2);f2(1,1);f2(2,1)];
    
    % nonlinear fit
    X=[ta(good1),data(good1,i)];
    modelfun='y~exp(b1*x1)*(b2+b3*x2)';
    fit2=fitnlm(X,er(good1),modelfun,a2g); 
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(er(good1));
    AIC=length(er(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    
     fprintf(fileID,['%s ',repmat('%g ',1,9),'\n'],cell2mat(dataset(:,i+1).Properties.VariableNames),a2,pval,varexpl,AIC,h);
    
    
    if(n==1 && AIC<AIC1) %fitat variable optimizing the AIC compared to step 1
        best=AIC;
        i1=i;
        n=0;
    elseif(n==0 && AIC<best) %variables optimizing the AIC compared to other varibles in the loop
        best=AIC;
        i1=i;
    end
end

end

fprintf(fileID,'\n');
fprintf(fileID,'selected: \n');



if(~isempty(myi))
    
    clear good1 a2 f2 fit2

    x1=data(:,myi(1));
    good1=~isnan(x1) & ~isnan(er) & ~isnan(ta);
    
    % compute the starting values for nonlinear fit
    er2=er(good1)./exp(a1(2).*ta(good1));
    f2=regress(er2,[ones(size(x1(good1))) x1(good1)]);
    a2g=[a1(2);f2(1,1);f2(2,1)];
    
    % nonlinear fit
    modelfun='y~exp(b1*x1)*(b2+b3*x2)';
    fit2=fitnlm([ta(good1),x1(good1)],er(good1),modelfun,a2g);
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(er(good1));
    AIC2=length(er(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,9),'\n'],cell2mat(dataset(:,myi(1)+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
    
elseif(isempty(myi) && n==0) % at least one additional varible optimizing the AIC compared to step 1
    
    clear good1 a2 f2 fit2
    addedvar(addedvar==i1)=[]; % delete the selected variable from the array of candidates

    x1=data(:,i1);
    good1=~isnan(x1) & ~isnan(er) & ~isnan(ta);
    
    % compute the starting values for nonlinear fit
    er2=er(good1)./exp(a1(2).*ta(good1));
    f2=regress(er2,[ones(size(x1(good1))) x1(good1)]);
    a2g=[a1(2);f2(1,1);f2(2,1)];
    
    % nonlinear fit
    modelfun='y~exp(b1*x1)*(b2+b3*x2)';
    fit2=fitnlm([ta(good1),x1(good1)],er(good1),modelfun,a2g);
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(er(good1));
    AIC2=length(er(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,9),'\n'],cell2mat(dataset(:,i1+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
    
else % no additional varible in the loop optimizing the AIC -> test Time_dec & DOY and stop
    
    clear good1 a2 f2 fit2
    i1=5; %Time_Dec
    x1=data(:,i1);
    good1=~isnan(x1) & ~isnan(er) & ~isnan(ta);
    
    % compute the starting values for nonlinear fit
    er2=er(good1)./exp(a1(2).*ta(good1));
    f2=regress(er2,[ones(size(x1(good1))) x1(good1)]);
    a2g=[a1(2);f2(1,1);f2(2,1)]; %guess
    
    % nonlinear fit
    modelfun='y~exp(b1*x1)*(b2+b3*x2)';
    fit2=fitnlm([ta(good1),x1(good1)],er(good1),modelfun,a2g);
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(er(good1));
    AIC2=length(er(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    if(AIC2<AIC1)
        
        % write fitted values and staistics on output file
        fprintf(fileID,['%s ',repmat('%g ',1,9),'\n'],cell2mat(dataset(:,i1+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
       
        i2=4; %DOY
        x2=data(:,i2);
        good2=(~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));

        % compute the starting values for nonlinear fit
        er2=(er(good2)./exp(a2(1).*ta(good2)))-a2(3).*x1(good2);
        f3=regress(er2,[ones(size(x2(good2))) x2(good2)]);
        a3g=[a2(1); f3(1,1); a2(3); f3(2,1)];

        % nonlinear fit
        modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3)';
        fit3=fitnlm([ta(good2),x1(good2),x2(good2)],er(good2),modelfun,a3g);
        a3=fit3.Coefficients.Estimate;
        pval=fit3.Coefficients.pValue;
        res=fit3.Residuals.Raw;
        varexpl=var(fit3.Fitted)/var(er(good2));
        AIC3=length(er(good2))*log(sum(res.^2))+2*length(a3);
        h=lillietest(res);
        
        if(AIC3<AIC2)
            fit2=fit3;
            good1=good2;
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);          
        end
    
    else
        
        clear good1 a2 f2 fit2
        i1=4; %DOY
        x1=data(:,i1);
        good1=~isnan(x1) & ~isnan(er) & ~isnan(ta);

        % compute the starting values for nonlinear fit
        er2=er(good1)./exp(a1(2).*ta(good1));
        f2=regress(er2,[ones(size(x1(good1))) x1(good1)]);
        a2g=[a1(2); f2(1,1);f2(2,1)]; %guess

        % nonlinear fit
        modelfun='y~exp(b1*x1)*(b2+b3*x2)';
        fit2=fitnlm([ta(good1),x1(good1)],er(good1),modelfun,a2g);
        a2=fit2.Coefficients.Estimate;
        pval=fit2.Coefficients.pValue;
        res=fit2.Residuals.Raw;
        varexpl=var(fit2.Fitted)/var(er(good1));
        AIC2=length(er(good1))*log(sum(res.^2))+2*length(a2);
        h=lillietest(res);
        
        if(AIC2<AIC1)
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,9),'\n'],cell2mat(dataset(:,i1+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
        end
       
 
    end    
    
    % plot results: measured vs fitted
    figure
    scatter(er(good1),fit2.Fitted,50,'r','filled')
    hold on
    plot(min(er):max(er),min(er):max(er),'k','LineWidth',1.5)
    set(gca,'FontSize',12)
    xlabel('measured ER (\mu mol/m^2/s)','FontSize',15)
    ylabel('modeled ER (\mu mol/m^2/s)','FontSize',15)
   
    fclose(fileID);
    
    return
end

 %% third step: ER vs additional variable #2
 
fprintf(fileID,'\n');
fprintf(fileID,'third step \n');
fprintf(fileID,'variable b_0 a_0 a_1 a_2 pvals varexpl AIC h \n');

n=1;
if(isempty(myi))
 for i=addedvar
    
    clear good2 X a3 f3 fit3
    good2=(~isnan(data(:,i)) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    %%guess
    er2=(er(good2)./exp(a2(1).*ta(good2)))-a2(3).*x1(good2);
    f3=regress(er2,[ones(size(data(good2,i))) data(good2,i)]);
    a3g=[a2(1); f3(1,1); a2(3); f3(2,1)];

    X=[ta(good2),x1(good2),data(good2,i)];
    
    %%fitnonlin
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3)';
    fit3=fitnlm(X,er(good2),modelfun,a3g);
    
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(er(good2));
    AIC=length(er(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
   
     fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i+1).Properties.VariableNames),a3,pval,varexpl,AIC,h);
    
    
    if(n==1 && AIC<AIC2) 
        best=AIC;
        i2=i;
        n=0;
    elseif(n==0 && AIC<best) 
        best=AIC;
        i2=i;
    end
 end
end
fprintf(fileID,'\n');
fprintf(fileID,'selected: \n');

if(~isempty(myi))

    clear good2 a3 f3 fit3
    
    x2=data(:,myi(2));
    good2=(~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    % compute the starting values for nonlinear fit
    er2=(er(good2)./exp(a2(1).*ta(good2)))-a2(3).*x1(good2);
    f3=regress(er2,[ones(size(x2(good2))) x2(good2)]);
    a3g=[a2(1); f3(1,1); a2(3); f3(2,1)];
    
    %%nonlinear fit
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3)';
    fit3=fitnlm([ta(good2),x1(good2),x2(good2)],er(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(er(good2));
    AIC3=length(er(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,myi(2)+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
    
elseif(isempty(myi) && n==0)
    
    clear good2 a3 f3 fit3
    addedvar(addedvar==i2)=[]; % delete the selected variable from the array of candidates

    x2=data(:,i2);
    good2=(~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    % compute the starting values for nonlinear fit
    er2=(er(good2)./exp(a2(1).*ta(good2)))-a2(3).*x1(good2);
    f3=regress(er2,[ones(size(x2(good2))) x2(good2)]);
    a3g=[a2(1); f3(1,1); a2(3); f3(2,1)];
    
    %%nonlinear fit
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3)';
    fit3=fitnlm([ta(good2),x1(good2),x2(good2)],er(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(er(good2));
    AIC3=length(er(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
    
else
    clear good2 a3 f3 fit3
    i2=5; %Time_Dec
    x2=data(:,i2);
    good2=(~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    % compute the starting values for nonlinear fit
    er2=(er(good2)./exp(a2(1).*ta(good2)))-a2(3).*x1(good2);
    f3=regress(er2,[ones(size(x2(good2))) x2(good2)]);
    a3g=[a2(1); f3(1,1); a2(3); f3(2,1)];
    
    % nonlinear fit
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3)';
    fit3=fitnlm([ta(good2),x1(good2),x2(good2)],er(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(er(good2));
    AIC3=length(er(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    if(AIC3<AIC2)
        % write fitted values and staistics on output file
        fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
        
        i3=4;
        x3=data(:,i3);
        good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));

        %%guess
        er2=(er(good3)./exp(a3(1).*ta(good3)))-a3(3).*x1(good3)-a3(4).*x2(good3);
        f4=regress(er2,[ones(size(x3(good3))) x3(good3)]);
        a4g=[a3(1); f3(1,1); a3(3); a3(4); f3(2,1)];

        %%fitnonlin
        modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
        fit4=fitnlm([ta(good3),x1(good3),x2(good3),x3(good3)],er(good3),modelfun,a4g);

        a4=fit4.Coefficients.Estimate;
        pval=fit4.Coefficients.pValue;
        res=fit4.Residuals.Raw;
        varexpl=var(fit4.Fitted)/var(er(good3));
        AIC4=length(er(good3))*log(sum(res.^2))+2*length(a4);
        h=lillietest(res);
    
        if(AIC4<AIC3)
            fit3=fit4;
            good2=good3;
            fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i3+1).Properties.VariableNames),a4,pval,varexpl,AIC4,h);
        end
         
    else
        clear good2 a3 f3 fit3
        i2=4; %DOY
        x2=data(:,i2);
        good2=(~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));

        % compute the starting values for nonlinear fit
        er2=(er(good2)./exp(a2(1).*ta(good2)))-a2(3).*x1(good2);
        f3=regress(er2,[ones(size(x2(good2))) x2(good2)]);
        a3g=[a2(1); f3(1,1); a2(3); f3(2,1)];

        % nonlinear fit
        modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3)';
        fit3=fitnlm([ta(good2),x1(good2),x2(good2)],er(good2),modelfun,a3g);
        a3=fit3.Coefficients.Estimate;
        pval=fit3.Coefficients.pValue;
        res=fit3.Residuals.Raw;
        varexpl=var(fit3.Fitted)/var(er(good2));
        AIC3=length(er(good2))*log(sum(res.^2))+2*length(a3);
        h=lillietest(res);
        
        if(AIC3<AIC2)
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
        end
        
    end
    
    % plot results: measured vs fitted
    figure
    scatter(er(good2),fit3.Fitted,50,'r','filled')
    hold on
    plot(min(er):max(er),min(er):max(er),'k','LineWidth',1.5)
    set(gca,'FontSize',12)
    xlabel('measured ER (\mu mol/m^2/s)','FontSize',15)
    ylabel('modeled ER (\mu mol/m^2/s)','FontSize',15)
   
    fclose(fileID);
    return
   
end

 %% fourth step: ER vs additional variable #3
 
fprintf(fileID,'\n');
fprintf(fileID,'fourth step \n');
fprintf(fileID,'variable b_0 a_0 a_1 a_2 a_3 pvals varexpl AIC h \n');

n=1;
if(isempty(myi))
 for i=addedvar
    
    clear good3 X a4 f4 fit4
    good3=(~isnan(data(:,i)) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    %%guess
    er2=(er(good3)./exp(a3(1).*ta(good3)))-a3(3).*x1(good3)-a3(4).*x2(good3);
    f4=regress(er2,[ones(size(data(good3,i))) data(good3,i)]);
    a4g=[a3(1); f4(1,1); a3(3); a3(4); f4(2,1)];

    X=[ta(good3),x1(good3),x2(good3),data(good3,i)];
    
    %%fitnonlin
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
    fit4=fitnlm(X,er(good3),modelfun,a4g);  
    a4=fit4.Coefficients.Estimate;
    pval=fit4.Coefficients.pValue;
    res=fit4.Residuals.Raw;
    varexpl=var(fit4.Fitted)/var(er(good3));
    AIC=length(er(good3))*log(sum(res.^2))+2*length(a4);
    h=lillietest(res);
    
    fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i+1).Properties.VariableNames),a4,pval,varexpl,AIC,h);
    
    if(n==1 && AIC<AIC3) 
        best=AIC;
        i3=i;
        n=0;
    elseif(n==0 && AIC<best) 
        best=AIC;
        i3=i;
    end
 end
end
fprintf(fileID,'\n');
fprintf(fileID,'selected: \n');

if(~isempty(myi))

    clear good3 a4 f4 fit4
    
    x3=data(:,myi(3));
    good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    %%guess
    er2=(er(good3)./exp(a3(1).*ta(good3)))-a3(3).*x1(good3)-a3(4).*x2(good3);
    f4=regress(er2,[ones(size(x3(good3))) x3(good3)]);
    a4g=[a3(1); f3(1,1); a3(3); a3(4); f3(2,1)];
    
    %%fitnonlin
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
    fit4=fitnlm([ta(good3),x1(good3),x2(good3),x3(good3)],er(good3),modelfun,a4g);
    a4=fit4.Coefficients.Estimate;
    pval=fit4.Coefficients.pValue;
    res=fit4.Residuals.Raw;
    varexpl=var(fit4.Fitted)/var(er(good3));
    AIC4=length(er(good3))*log(sum(res.^2))+2*length(a4);
    h=lillietest(res);
    
    fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,myi(3)+1).Properties.VariableNames),a4,pval,varexpl,AIC4,h);
    
elseif(isempty(myi) && n==0)
    
    clear good3 a4 f4 fit4
    addedvar(addedvar==i3)=[]; % delete the selected variable from the array of candidates

    x3=data(:,i3);
    good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    %%guess
    er2=(er(good3)./exp(a3(1).*ta(good3)))-a3(3).*x1(good3)-a3(4).*x2(good3);
    f4=regress(er2,[ones(size(x3(good3))) x3(good3)]);
    a4g=[a3(1); f3(1,1); a3(3); a3(4); f3(2,1)];
    
    %%fitnonlin
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
    fit4=fitnlm([ta(good3),x1(good3),x2(good3),x3(good3)],er(good3),modelfun,a4g);
    a4=fit4.Coefficients.Estimate;
    pval=fit4.Coefficients.pValue;
    res=fit4.Residuals.Raw;
    varexpl=var(fit4.Fitted)/var(er(good3));
    AIC4=length(er(good3))*log(sum(res.^2))+2*length(a4);
    h=lillietest(res);
    
    fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i3+1).Properties.VariableNames),a4,pval,varexpl,AIC4,h);
    
else
    
    clear good3 a4 f4 fit4
    i3=5;
    x3=data(:,i3);
    good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
    %%guess
    er2=(er(good3)./exp(a3(1).*ta(good3)))-a3(3).*x1(good3)-a3(4).*x2(good3);
    f4=regress(er2,[ones(size(x3(good3))) x3(good3)]);
    a4g=[a3(1); f3(1,1); a3(3); a3(4); f3(2,1)];
    
    %%fitnonlin
    modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
    fit4=fitnlm([ta(good3),x1(good3),x2(good3),x3(good3)],er(good3),modelfun,a4g);  
    a4=fit4.Coefficients.Estimate;
    pval=fit4.Coefficients.pValue;
    res=fit4.Residuals.Raw;
    varexpl=var(fit4.Fitted)/var(er(good3));
    AIC4=length(er(good3))*log(sum(res.^2))+2*length(a4);
    h=lillietest(res);
    
    if(AIC4<AIC3)
        fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i3+1).Properties.VariableNames),a4,pval,varexpl,AIC4,h);
        
        i4=4;
        x4=data(:,i4);
        good4=(~isnan(x4) & ~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
    
        %%guess
        er2=(er(good4)./exp(a4(1).*ta(good4)))-a4(3).*x1(good4)-a4(4).*x2(good4)-a4(5).*x3(good4);
        f5=regress(er2,[ones(size(x4(good4))) x4(good4)]);
        a5g=[a4(1); f5(1,1); a4(3); a4(4); a4(5); f5(2,1)];
    
        %%fitnonlin
        modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4+b6*x5)';
        fit5=fitnlm([ta(good4),x1(good4),x2(good4),x3(good4),x4(good4)],er(good4),modelfun,a5g);
        a5=fit5.Coefficients.Estimate;
        pval=fit5.Coefficients.pValue;
        res=fit5.Residuals.Raw;
        varexpl=var(fit5.Fitted)/var(er(good4));
        AIC5=length(er(good4))*log(sum(res.^2))+2*length(a5);
        h=lillietest(res);
        
        if(AIC5<AIC4)
            fit4=fit5;
            good3=good4;
            fprintf(fileID,['%s ',repmat('%g ',1,15),'\n'],cell2mat(dataset(:,i4+1).Properties.VariableNames),a5,pval,varexpl,AIC5,h);
        end
    
    else
        
        clear good3 a4 f4 fit4
        i3=4;
        x3=data(:,i3);
        good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));

        %%guess
        er2=(er(good3)./exp(a3(1).*ta(good3)))-a3(3).*x1(good3)-a3(4).*x2(good3);
        f4=regress(er2,[ones(size(x3(good3))) x3(good3)]);
        a4g=[a3(1); f3(1,1); a3(3); a3(4); f3(2,1)];

        %%fitnonlin
        modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
        fit4=fitnlm([ta(good3),x1(good3),x2(good3),x3(good3)],er(good3),modelfun,a4g);
        a4=fit4.Coefficients.Estimate;
        pval=fit4.Coefficients.pValue;
        res=fit4.Residuals.Raw;
        varexpl=var(fit4.Fitted)/var(er(good3));
        AIC4=length(er(good3))*log(sum(res.^2))+2*length(a4);
        h=lillietest(res);
    
        if(AIC4<AIC3)
            fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i3+1).Properties.VariableNames),a4,pval,varexpl,AIC4,h);
        end
    
    end
    % plot results: measured vs fitted
    figure
    scatter(er(good3),fit4.Fitted,50,'r','filled')
    hold on
    plot(min(er):max(er),min(er):max(er),'k','LineWidth',1.5)
    set(gca,'FontSize',12)
    xlabel('measured ER (\mu mol/m^2/s)','FontSize',15)
    ylabel('modeled ER (\mu mol/m^2/s)','FontSize',15)
   
    fclose(fileID);
    return
end
%% measured vs fitted data

% plot results when DOY has not been tested (see optional): measured vs fitted 
figure
scatter(er(good3),fit4.Fitted,50,'r','filled')
hold on
plot(min(er):max(er),min(er):max(er),'k','LineWidth',1.5)
set(gca,'FontSize',12)
xlabel('measured ER (\mu mol/m^2/s)','FontSize',15)
ylabel('modeled ER (\mu mol/m^2/s)','FontSize',15)

fclose(fileID);
%% optional
% copy and past the fourth step here to test additional variable #4 
% changes 
% good2-> good3=(~isnan(data(:,i)) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
% er2=er(good3)./(a3(1)*a3(2).*ta(good3)./(a3(1)+a3(2).*ta(good3)))-a3(4).*x1(good3)-a3(5).*x1(good3);
% X=[ta(good3),x1(good3),x2(good3),data(good3,i)];
% modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3+b6*x4)';
% a3->a4 f3->f4 fit3->f4
% AIC2->AIC3; AIC3->AIC4

%% recommended: check parameter stability
% good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(er) & ~isnan(ta));
% a4g=a4;
% modelfun='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
% fit4=fitnlm([ta(good3),x1(good3),x2(good3),x3(good3)],er(good3),modelfun,a4g);
% a4=fit4.Coefficients.Estimate;
% pval=fit4.Coefficients.pValue;
% res=fit4.Residuals.Raw;
% varexpl=var(fit4.Fitted)/var(er(good3));
% AIC4=length(er(good3))*log(sum(res.^2))+2*length(a4);
% h=lillietest(res);