function [effluxReflux,mx] = refluxEfflux(tef,pairs,ind,dataDay)

%% INPUTS:
%1) tef ? tef output structure from calcTEF.m

%2) paris ? pairs of tef transects to calculate reflux/efflux between (1st column is 
%oceanward transect, 2nd column is glacierward transect.)

%3) ind ? indices to calculate reflux/efflux if not using full time series (otherwise can
%cause errors at beginning of run when tef time series aren't equilibrated)

%4) dataDay ? day to pull model data from

%OUTPUTS
%1) effluxReflux ? snapshot of effluxReflux from dataDay
%2) mx ? Full time series of efflux/reflux for all indices specified by ind

%%
[M,N] = size(pairs);
for i = 1:M
% for i = 1;
    
    mx(i).Q1 = abs(tef(pairs(i,1)).Qin);
    mx(i).Q2 = abs(tef(pairs(i,1)).Qout);
    mx(i).Q3 = abs(tef(pairs(i,2)).Qin);
    mx(i).Q4 = abs(tef(pairs(i,2)).Qout);
    
    mx(i).S1 = tef(pairs(i,1)).Sin;
    mx(i).S2 = tef(pairs(i,1)).Sout;
    mx(i).S3 = tef(pairs(i,2)).Sin;
    mx(i).S4 = tef(pairs(i,2)).Sout;
    
    mx(i).Qs1 = abs(tef(pairs(i,1)).Qs_in);
    mx(i).Qs2 = abs(tef(pairs(i,1)).Qs_out);
    mx(i).Qs3 = abs(tef(pairs(i,2)).Qs_in);
    mx(i).Qs4 = abs(tef(pairs(i,2)).Qs_out);
   %% 
   if prod(size(ind)) == 0
       ind = 36:length(mx(i).Q1)-36; %default to entire time series (except edges effects of filter)
   end
       
    for ii = ind
    
        if mx(i).Qs1(ii) + mx(i).Qs4(ii) ~= mx(i).Qs2(ii) + mx(i).Qs3(ii)
                cmp1 = 1; cmp2 = 1; cmp3 = 1; cmp4 = 1;
                
                %pre-allocate adjusted values as original values as default
                mx(i).S1adj(ii) = mx(i).S1(ii);
                mx(i).S2adj(ii) = mx(i).S2(ii);
                mx(i).S3adj(ii) = mx(i).S3(ii);
                mx(i).S4adj(ii) = mx(i).S4(ii);
                
                mx(i).Q1adj(ii) = mx(i).Q1(ii);
                mx(i).Q2adj(ii) = mx(i).Q2(ii);
                mx(i).Q3adj(ii) = mx(i).Q3(ii);
                mx(i).Q4adj(ii) = mx(i).Q4(ii);
                
                mx(i).S1min(ii,:) = NaN(1,4);
                mx(i).S2min(ii,:) = NaN(1,4);
                mx(i).S3min(ii,:) = NaN(1,4);
                
                mx(i).S2max(ii,:) = NaN(1,4);
                mx(i).S3max(ii,:) = NaN(1,4);
                mx(i).S4max(ii,:) = NaN(1,4);
                
                %find all min and max values for each salinity inequality. Adjust
                %salinities as necessary. Adjust associated volume fluxes so that
                %Qs = Qadj+Sadj
                mn = mean([mx(i).S4adj(ii) mx(i).S2adj(ii)]);
                df2 = diff([mx(i).S2adj(ii) mn]);
                df4 = diff([mx(i).S4adj(ii) mn]);
                mx(i).S2min(ii,1) = mx(i).S2adj(ii) + df2;
                mx(i).S4max(ii,1) = mx(i).S4adj(ii) + df4;
                if mx(i).S4adj(ii) > mx(i).S2adj(ii)
                    cmp1 = 0;
                    mx(i).S2adj(ii) = mx(i).S2adj(ii) + df2 + 0.0001;
                    mx(i).S4adj(ii) = mx(i).S4adj(ii) + df4 - 0.0001;
                    mx(i).Q4adj(ii) = mx(i).Qs4(ii)/mx(i).S4adj(ii);  %enforce Qs = Q*S
                    mx(i).Q2adj(ii) = mx(i).Qs2(ii)/mx(i).S2adj(ii);
                end
                
                mn = mean([mx(i).S2adj(ii) mx(i).S1adj(ii)]);
                df1 = diff([mx(i).S1adj(ii) mn]);
                df2 = diff([mx(i).S2adj(ii) mn]);
                mx(i).S2max(ii,2) = mx(i).S2adj(ii) + df2;
                mx(i).S1min(ii,2) = mx(i).S1adj(ii) + df1;   
                if mx(i).S2adj(ii) > mx(i).S1adj(ii)
                    cmp2 = 0;
                    mx(i).S1adj(ii) = mx(i).S1adj(ii) + df1 + 0.0001;
                    mx(i).S2adj(ii) = mx(i).S2adj(ii) + df2 - 0.0001;
                    mx(i).Q2adj(ii) = mx(i).Qs2(ii)/mx(i).S2adj(ii);
                    mx(i).Q1adj(ii) = mx(i).Qs1(ii)/mx(i).S1adj(ii);
                end
                
                mn = mean([mx(i).S4adj(ii) mx(i).S3adj(ii)]);
                df4 = diff([mx(i).S4adj(ii) mn]);
                df3 = diff([mx(i).S3adj(ii) mn]);
                mx(i).S3min(ii,3) = mx(i).S3adj(ii) + df3;
                mx(i).S4max(ii,3) = mx(i).S4adj(ii) + df4;    
                if mx(i).S4adj(ii) > mx(i).S3adj(ii)
                    cmp3 = 0;  
                    mx(i).S3adj(ii) = mx(i).S3adj(ii) + df3 + 0.0001;
                    mx(i).S4adj(ii) = mx(i).S4adj(ii) + df4 - 0.0001;
                    mx(i).Q3adj(ii) = mx(i).Qs3(ii)/mx(i).S3adj(ii); 
                    mx(i).Q4adj(ii) = mx(i).Qs4(ii)/mx(i).S4adj(ii);
                end
                
                mn = mean([mx(i).S3adj(ii) mx(i).S1adj(ii)]);
                df1 = diff([mx(i).S1adj(ii) mn]);
                df3 = diff([mx(i).S3adj(ii) mn]);
                mx(i).S1min(ii,4) = mx(i).S1adj(ii) + df1;
                mx(i).S3max(ii,4) = mx(i).S3adj(ii) + df3;
                if mx(i).S3adj(ii) > mx(i).S1adj(ii) 
                    cmp4 = 0;
                    mx(i).S1adj(ii) = mx(i).S1adj(ii) + df1 + 0.0001;
                    mx(i).S3adj(ii) = mx(i).S3adj(ii) + df3 - 0.0001;
                    mx(i).Q3adj(ii) = mx(i).Qs3(ii)/mx(i).S3adj(ii); 
                    mx(i).Q1adj(ii) = mx(i).Qs1(ii)/mx(i).S1adj(ii);         
                end
                
                %check Q1+Q4=Q2+Q3. Further adjust volume fluxes if
                %necessary to preserve volume conservation
                
                if mx(i).Q1adj(ii) + mx(i).Q4adj(ii) ~= mx(i).Q2adj(ii) + mx(i).Q3adj(ii)
                    
                    mn = mean([mx(i).Q1adj(ii) + mx(i).Q4adj(ii), mx(i).Q2adj(ii) + mx(i).Q3adj(ii)]);
                    dfa = diff([mx(i).Q1adj(ii) + mx(i).Q4adj(ii) mn]);
                    dfb = diff([mx(i).Q2adj(ii) + mx(i).Q3adj(ii) mn]);

                    mx(i).Q1adj(ii) = mx(i).Q1adj(ii) + dfa / 2;
                    mx(i).Q4adj(ii) = mx(i).Q4adj(ii) + dfa / 2;
                    mx(i).Q2adj(ii) = mx(i).Q2adj(ii) + dfb / 2;
                    mx(i).Q3adj(ii) = mx(i).Q3adj(ii) + dfb / 2;
                end
                
                %adjust salt fluxes to conserve mass and volume (Q1S1 + Q4S4 = Q2S2 + Q3S3) by
                %adjusting salinity within each appropriate range. Try to split the
                %difference bettwen adjustments, but if adjusted salinity falls
                %outside of appropriate range, then difference is made up
                %from accompanying salinity
                
                if mx(i).Qs1(ii) + mx(i).Qs4(ii) ~= mx(i).Qs2(ii) + mx(i).Qs3(ii)
                    mn = mean([mx(i).S1adj(ii)*mx(i).Q1adj(ii)+mx(i).S4adj(ii)*mx(i).Q4adj(ii), ...
                        mx(i).S2adj(ii)*mx(i).Q2adj(ii)+mx(i).S3adj(ii)*mx(i).Q3adj(ii)]);
                    dfa = diff([mx(i).S1adj(ii)*mx(i).Q1adj(ii)+mx(i).S4adj(ii)*mx(i).Q4adj(ii) mn]);
                    dfb = diff([mx(i).S2adj(ii)*mx(i).Q2adj(ii)+mx(i).S3adj(ii)*mx(i).Q3adj(ii) mn]);

                    %start with right side of equation
                    s1 = mx(i).S1adj(ii) + (dfa/2)/mx(i).Q1adj(ii);
                    s4 = mx(i).S4adj(ii) + (dfa/2)/mx(i).Q4adj(ii);
                    
                    if s1 > nanmax(mx(i).S1min(ii,:)) && s4 < nanmin(mx(i).S4max(ii,:))
                        mx(i).S1adj(ii) = s1;
                        mx(i).S4adj(ii) = s4;
                    elseif s4 > nanmin(mx(i).S4max(ii,:))
                        mx(i).S4adj(ii) = nanmin(mx(i).S4max(ii,:)) - 0.0001;
                        mx(i).S1adj(ii) = (abs(mx(i).S4adj(ii)*mx(i).Q4adj(ii) - mn))/mx(i).Q1adj(ii);
                    elseif s1 < nanmax(mx(i).S1min(ii,:))
                        mx(i).S1adj(ii) = nanmax(mx(i).S1min(ii,:)) + 0.0001;
                        mx(i).S4adj(ii) = (abs(mx(i).S1adj(ii)*mx(i).Q1adj(ii) - mn))/mx(i).Q4adj(ii);
                    else 
                        error('Debug: Using mean value will not work to balance equation')
                    end
                    
                    s2 = mx(i).S2adj(ii) + (dfb/2)/mx(i).Q2adj(ii);
                    s3 = mx(i).S3adj(ii) + (dfb/2)/mx(i).Q3adj(ii);
                    
                    if s2 > nanmax(mx(i).S2min(ii,:)) && s2 < nanmin(mx(i).S2max(ii,:)) && s3 > nanmax(mx(i).S3min(ii,:)) && s3 < nanmin(mx(i).S3max(ii,:))
                        mx(i).S2adj(ii) = s2;
                        mx(i).S3adj(ii) = s3;
                    elseif s2 < nanmax(mx(i).S2min(ii,:))
                        mx(i).S2adj(ii) = nanmax(mx(i).S2min(ii,:)) + 0.0001;
                        mx(i).S3adj(ii) = (abs(mx(i).S2adj(ii)*mx(i).Q2adj(ii) - mn))/mx(i).Q3adj(ii);
                    elseif s2 > nanmin(mx(i).S2max(ii,:))
                        mx(i).S2adj(ii) = nanmin(mx(i).S2max(ii,:)) - 0.0001;
                        mx(i).S3adj(ii) = (abs(mx(i).S2adj(ii)*mx(i).Q2adj(ii) - mn))/mx(i).Q3adj(ii);
                    elseif s3 < nanmax(mx(i).S3min(ii,:))
                        mx(i).S3adj(ii) = nanmax(mx(i).S3min(ii,:)) + 0.0001;
                        mx(i).S2adj(ii) = (abs(mx(i).S3adj(ii)*mx(i).Q3adj(ii) - mn))/mx(i).Q2adj(ii);
                    elseif s3 > nanmin(mx(i).S3max(ii,:))
                        mx(i).S3adj(ii) = nanmin(mx(i).S3max(ii,:)) - 0.0001;
                        mx(i).S2adj(ii) = (abs(mx(i).S3adj(ii)*mx(i).Q3adj(ii) - mn))/mx(i).Q2adj(ii);
                    else
                        error('Debug: Using mean value will not work to balance equation')
                    end
                    
                    %calculate adjusted salt fluxes
                    mx(i).Qs1adj(ii) = mx(i).Q1adj(ii) * mx(i).S1adj(ii);
                    mx(i).Qs2adj(ii) = mx(i).Q2adj(ii) * mx(i).S2adj(ii);
                    mx(i).Qs3adj(ii) = mx(i).Q3adj(ii) * mx(i).S3adj(ii);
                    mx(i).Qs4adj(ii) = mx(i).Q4adj(ii) * mx(i).S4adj(ii);
                    
                    %test
                    mx(i).LHS(ii) = mx(i).Qs1adj(ii) + mx(i).Qs4adj(ii);
                    mx(i).RHS(ii) = mx(i).Qs2adj(ii) + mx(i).Qs3adj(ii);
                    
                end
                
            mx(i).adjPercQs(ii,:) = [abs(mx(i).Qs1adj(ii) - mx(i).Qs1(ii))/mx(i).Qs1(ii) ... 
                abs(mx(i).Qs2adj(ii) - mx(i).Qs2(ii))/mx(i).Qs2(ii) ...
                abs(mx(i).Qs3adj(ii) - mx(i).Qs3(ii))/mx(i).Qs3(ii) ...
                abs(mx(i).Qs4adj(ii) - mx(i).Qs4(ii))/mx(i).Qs4(ii)];
            
            mx(i).adjPercQ(ii,:) = [abs(mx(i).Q1adj(ii) - mx(i).Q1(ii))/mx(i).Q1(ii) ... 
                abs(mx(i).Q2adj(ii) - mx(i).Q2(ii))/mx(i).Q2(ii) ...
                abs(mx(i).Q3adj(ii) - mx(i).Q3(ii))/mx(i).Q3(ii) ...
                abs(mx(i).Q4adj(ii) - mx(i).Q4(ii))/mx(i).Q4(ii)];
            
            mx(i).adjPercS(ii,:) = [abs(mx(i).S1adj(ii) - mx(i).S1(ii))/mx(i).S1(ii) ... 
                abs(mx(i).S2adj(ii) - mx(i).S2(ii))/mx(i).S2(ii) ...
                abs(mx(i).S3adj(ii) - mx(i).S3(ii))/mx(i).S3(ii) ...
                abs(mx(i).S4adj(ii) - mx(i).S4(ii))/mx(i).S4(ii)];
        end
  
        mx(i).a21(ii) = (mx(i).Q2adj(ii)./mx(i).Q1adj(ii)).*((mx(i).S2adj(ii)-mx(i).S4adj(ii))./(mx(i).S1adj(ii)-mx(i).S4adj(ii)));
        mx(i).a34(ii) = (mx(i).Q3adj(ii)./mx(i).Q4adj(ii)).*((mx(i).S1adj(ii)-mx(i).S3adj(ii))./(mx(i).S1adj(ii)-mx(i).S4adj(ii)));
        mx(i).a31(ii) = (mx(i).Q3adj(ii)./mx(i).Q1adj(ii)).*((mx(i).S3adj(ii)-mx(i).S4adj(ii))./(mx(i).S1adj(ii)-mx(i).S4adj(ii)));
        mx(i).a24(ii) = (mx(i).Q2adj(ii)./mx(i).Q4adj(ii)).*((mx(i).S1adj(ii)-mx(i).S2adj(ii))./(mx(i).S1adj(ii)-mx(i).S4adj(ii)));
        
        %compute absolute errors
        mx(i).Q1err(ii) = abs(mx(i).Q1adj(ii)-mx(i).Q1(ii));
        mx(i).Q2err(ii) = abs(mx(i).Q2adj(ii)-mx(i).Q2(ii));
        mx(i).Q3err(ii) = abs(mx(i).Q3adj(ii)-mx(i).Q3(ii));
        mx(i).Q4err(ii) = abs(mx(i).Q4adj(ii)-mx(i).Q4(ii));
        
        mx(i).S1err(ii) = abs(mx(i).S1adj(ii)-mx(i).S1(ii));
        mx(i).S2err(ii) = abs(mx(i).S2adj(ii)-mx(i).S2(ii));
        mx(i).S3err(ii) = abs(mx(i).S3adj(ii)-mx(i).S3(ii));
        mx(i).S4err(ii) = abs(mx(i).S4adj(ii)-mx(i).S4(ii));
        
        %compute total errors [ alpha = (Qa/Qb)*((Sa-Sb)/(Sc-Sd)) ]
        %a21
        Qa = mx(i).Q2adj(ii); 
        QaErr = mx(i).Q2err(ii);
        Qb = mx(i).Q1adj(ii);
        QbErr = mx(i).Q1err(ii);
        Sa = mx(i).S2adj(ii);
        SaErr = mx(i).S2err(ii);
        Sb = mx(i).S4adj(ii);
        SbErr = mx(i).S4err(ii);
        Sc = mx(i).S1adj(ii);
        ScErr = mx(i).S1err(ii);
        Sd = mx(i).S4err(ii);
        SdErr = mx(i).S4err(ii);
        alpha = mx(i).a21(ii);
        alphaType = 'a21';
        [mx(i).a21err(ii)] = calcError(Qa,QaErr,Qb,QbErr,Sa,SaErr,Sb,SbErr,Sc,ScErr,Sd,SdErr,alpha,alphaType);
        
        %a31
        Qa = mx(i).Q3adj(ii); 
        QaErr = mx(i).Q3err(ii);
        Qb = mx(i).Q1adj(ii);
        QbErr = mx(i).Q1err(ii);
        Sa = mx(i).S3adj(ii);
        SaErr = mx(i).S3err(ii);
        Sb = mx(i).S4adj(ii);
        SbErr = mx(i).S4err(ii);
        Sc = mx(i).S1adj(ii);
        ScErr = mx(i).S1err(ii);
        Sd = mx(i).S4err(ii);
        SdErr = mx(i).S4err(ii);
        alpha = mx(i).a31(ii);
        alphaType = 'a31';
        [mx(i).a31err(ii)] = calcError(Qa,QaErr,Qb,QbErr,Sa,SaErr,Sb,SbErr,Sc,ScErr,Sd,SdErr,alpha,alphaType);
        
        %a24
        Qa = mx(i).Q2adj(ii); 
        QaErr = mx(i).Q2err(ii);
        Qb = mx(i).Q4adj(ii);
        QbErr = mx(i).Q4err(ii);
        Sa = mx(i).S1adj(ii);
        SaErr = mx(i).S1err(ii);
        Sb = mx(i).S2adj(ii);
        SbErr = mx(i).S2err(ii);
        Sc = mx(i).S1adj(ii);
        ScErr = mx(i).S1err(ii);
        Sd = mx(i).S4err(ii);
        SdErr = mx(i).S4err(ii);
        alpha = mx(i).a24(ii);
        alphaType = 'a24';
        [mx(i).a24err(ii)] = calcError(Qa,QaErr,Qb,QbErr,Sa,SaErr,Sb,SbErr,Sc,ScErr,Sd,SdErr,alpha,alphaType);
        
        %a34
        Qa = mx(i).Q3adj(ii); 
        QaErr = mx(i).Q3err(ii);
        Qb = mx(i).Q4adj(ii);
        QbErr = mx(i).Q4err(ii);
        Sa = mx(i).S1adj(ii);
        SaErr = mx(i).S1err(ii);
        Sb = mx(i).S3adj(ii);
        SbErr = mx(i).S3err(ii);
        Sc = mx(i).S1adj(ii);
        ScErr = mx(i).S1err(ii);
        Sd = mx(i).S4err(ii);
        SdErr = mx(i).S4err(ii);
        alpha = mx(i).a34(ii);
        alphaType = 'a34';
        [mx(i).a34err(ii)] = calcError(Qa,QaErr,Qb,QbErr,Sa,SaErr,Sb,SbErr,Sc,ScErr,Sd,SdErr,alpha,alphaType);
    
    end
       
        %output synthesis for data day
        effluxReflux(i).a21 = nanmean(mx(i).a21(dataDay*24:dataDay*24+23));
        effluxReflux(i).a34 = nanmean(mx(i).a34(dataDay*24:dataDay*24+23));
        effluxReflux(i).a31 = nanmean(mx(i).a31(dataDay*24:dataDay*24+23));
        effluxReflux(i).a24 = nanmean(mx(i).a24(dataDay*24:dataDay*24+23));
        effluxReflux(i).a21err = nanmean(mx(i).a21err(dataDay*24:dataDay*24+23));
        effluxReflux(i).a34err = nanmean(mx(i).a34err(dataDay*24:dataDay*24+23));
        effluxReflux(i).a31err = nanmean(mx(i).a31err(dataDay*24:dataDay*24+23));
        effluxReflux(i).a24err = nanmean(mx(i).a24err(dataDay*24:dataDay*24+23));
        effluxReflux(i).Q1 = nanmean(mx(i).Q1adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Q2 = nanmean(mx(i).Q2adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Q3 = nanmean(mx(i).Q3adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Q4 = nanmean(mx(i).Q4adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).S1 = nanmean(mx(i).S1adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).S2 = nanmean(mx(i).S2adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).S3 = nanmean(mx(i).S3adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).S4 = nanmean(mx(i).S4adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Qs1 = nanmean(mx(i).Qs1adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Qs2 = nanmean(mx(i).Qs2adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Qs3 = nanmean(mx(i).Qs3adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).Qs4 = nanmean(mx(i).Qs4adj(dataDay*24:dataDay*24+23));
        effluxReflux(i).attributes = ...
['========================================          ' ; 
'                 |     a24       |                ' ; 
'     Q4---->     |-------------->|   Q2---->      ' ; 
'              S4 |    \     /    | S2             ' ; 
'land  - - - - - -| a34 |   | a21 | - - - - - ocean' ; 
'                 |    /     \    |                ' ; 
'     <----Q3     |<--------------|   <----Q1      ' ; 
'              S3 |      a31      | S1             ' ; 
'    =========================================     '] ;

end


function [error] = calcError(Qa,QaErr,Qb,QbErr,Sa,SaErr,Sb,SbErr,Sc,ScErr,Sd,SdErr,alpha,alphaType)
%inputs given for equation: alpha = (Qa/Qb)*((Sa-Sb)/(Sc-Sd)). Either Sc or
%Sd will be repeated, but the calcError takes care of this issue.
%Input/output errors are absolute.

    if alphaType(end) == '4' %a24 and a34
        % Form1: F = (a/b)*((c-d)/(c-e))
            a = Qa;
            b = Qb;
            c = Sa;
            d = Sb;
            e = Sd;
        
            aErr = QaErr;
            bErr = QbErr;
            cErr = SaErr;
            dErr = SbErr;
            eErr = SdErr;
    
            %partial derivatives
            dFda = (1/b)*((c-d)/(c-e));
            dFdb = -(a/(b^2))*((c-d)/(c-e));
            dFdc = (a/b)*(-(c-d)/((c-e)^2) + 1/(c-e));
            dFdd = -(a/b)*(1/(c-e));
            dFde = -(a/b)*((c-d)/((c-e)^2));

    elseif alphaType(end) == '1' %a21 and a31
            %Form 2: F = (a/b)*((c-d)/(e-d))
            a = Qa;
            b = Qb;
            c = Sa;
            d = Sb;
            e = Sc;
        
            aErr = QaErr;
            bErr = QbErr;
            cErr = SaErr;
            dErr = SbErr;
            eErr = ScErr;
    
            %partial derivatives
            dFda = (1/b)*((c-d)/(e-d));
            dFdb = -(a/(b^2))*((c-d)/(e-d));
            dFdc = (a/b)*(1/(e-d));
            dFdd = (a/b)*((c-d)/((e-d)^2) - 1/(e-d));
            dFde = -(a/b)*((c-d)/((e-d)^2));
    end
            
            %relative errors
            Ra = aErr/abs(a);
            Rb = bErr/abs(b);
            Rc = cErr/abs(c);
            Rd = dErr/abs(d);
            Re = eErr/abs(e);
        
            %absolute error for alpha coefficient
            error = alpha*sqrt((dFda*Ra)^2 + (dFdb*Rb)^2 + (dFdc*Rc)^2 + (dFdd*Rd)^2 + (dFde*Re)^2);
end
    
  
end

