function [kappa, kappastd] = kappacalc(variables,error)

%This function was written by Kelly Karl January 2022

%This function calculates the area under the curves, the error, and the
%value for kappa and its error based on the paper "Quantitative Assessment of Ligand Bias from Bias Plots:
%The Bias Coefficient Kappa" by Karl and Hristova.

%To use enter the fit paramters in the following manner ([Etopligx, KAligx,
%nligx, Etopligy, KAligy, nligy, Etoprefx, KArefx, nrefx, Etoprefy, KArefy,
%nrefy](errors in same order)]

%This section takes the values that the usere has entered and assigns
%them to different variable names
EtopA=variables(1);
KA=variables(2);
nA=variables(3);
EtopB=variables(4);
KB=variables(5);
nB=variables(6);

EtopAref=variables(7);
KAref=variables(8);
nAref=variables(9);
EtopBref=variables(10);
KBref=variables(11);
nBref=variables(12);

%This section takes the error the user has entered and assigns it
deltaEtopA=error(1);
deltaKA=error(2);
deltanA=error(3);
deltaEtopB=error(4);
deltaKB=error(5);
deltanB=error(6);

deltaEtopAref=error(7);
deltaKAref=error(8);
deltanAref=error(9);
deltaEtopBref=error(10);
deltaKBref=error(11);
deltanBref=error(12);

%Defines the x-value to terminate the area calculation
n=min(EtopA,EtopAref);




%Calculate average area under both curves
x=0:.001:n-.001;
l=length(x);
yval=[]; yvalref=[];
for i=1:l
    z=x(i);
    y=EtopB*((z*KA/(EtopA-z))^(nB/nA))/(KB+((z*KA/(EtopA-z))^(nB/nA)));
    yref=EtopBref*((z*KAref/(EtopAref-z))^(nBref/nAref))/(KBref+((z*KAref/(EtopAref-z))^(nBref/nAref)));
    yval(end+1)=y;
    yvalref(end+1)=yref;

end

yval;
areaMC=[];
tarea=[]; tarearef=[];
for j=1:l-1
   minnie=min(yval(j),yval(j+1));
   minnieref=min(yvalref(j),yvalref(j+1));
   area=0.001*minnie+0.5*0.001*abs(yval(j+1)-yval(j));
   arearef=0.001*minnieref+0.5*0.001*abs(yvalref(j+1)-yvalref(j));
   tarea(end+1)=area;
   tarearef(end+1)=arearef;
end


totarea=sum(tarea);
totarearef=sum(tarearef);
rng('default');
kappa=[];
%Calculate the contribution of each error point to the curve
EtopAer=EtopA-deltaEtopA;
EtopAer2=EtopA+deltaEtopA;
KAer=KA-deltaKA;
KAer2=KA+deltaKA;
EtopBer=EtopB-deltaEtopB;
EtopBer2=EtopB+deltaEtopB;
KBer=KB-deltaKB
KBer2=KB+deltaKB;
nAer=nA-deltanA;
nAer2=nA+deltanA;
nBer=nB-deltanB;
nBer2=nB+deltanB;
m=min(EtopA,EtopAref);

EtopAerref=EtopAref-deltaEtopAref;
EtopAer2ref=EtopAref+deltaEtopAref;
KAerref=KAref-deltaKAref;
KAer2ref=KAref+deltaKAref;
EtopBerref=EtopBref-deltaEtopBref;
EtopBer2ref=EtopBref+deltaEtopBref;
KBerref=KBref-deltaKBref;
KBer2ref=KBref+deltaKBref;
nAerref=nAref-deltanAref;
nAer2ref=nAref+deltanAref;
nBerref=nBref-deltanBref;
nBer2ref=nBref+deltanBref;
m2=min(EtopAer, EtopAerref); 


    q2=0:0.001:m2-.001;
    q=0:0.001:m-.001;
    l2=length(q);
    l3=length(q2);
   
    yvaler1=[]; yvaler2=[]; yvaler3=[]; yvaler4=[]; yvaler5=[]; yvaler6=[]; yvaler7=[]; yvaler8=[]; yvaler9=[]; yvaler10=[]; yvaler11=[]; yvaler12=[];
    yvaler1ref=[]; yvaler2ref=[]; yvaler3ref=[]; yvaler4ref=[]; yvaler5ref=[]; yvaler6ref=[]; yvaler7ref=[]; yvaler8ref=[]; yvaler9ref=[]; yvaler10ref=[]; yvaler11ref=[]; yvaler12ref=[];
    
    if EtopAer<EtopAerref
        
        for k=1:l3
        z=q2(k);
        yer1=EtopB*((z*KA/(EtopAer-z))^(nB/nA))/(KB+((z*KA/(EtopAer-z))^(nB/nA)));
        yvaler1(end+1)=yer1;
        end
        
        for k=1:l3
        z=q2(k);
        yer1ref=EtopBref*((z*KAref/(EtopAerref-z))^(nBref/nAref))/(KBref+((z*KAref/(EtopAerref-z))^(nBref/nAref)));
        yvaler1ref(end+1)=yer1ref;
        end
        
    end
    

 
    
    for c=1:l2
    z=q(c);
    yer2=EtopB*((z*KA/(EtopAer2-z))^(nB/nA))/(KB+((z*KA/(EtopAer2-z))^(nB/nA)));
    yer3=EtopB*((z*KAer/(EtopA-z))^(nB/nA))/(KB+((z*KAer/(EtopA-z))^(nB/nA)));
    yer4=EtopB*((z*KAer2/(EtopA-z))^(nB/nA))/(KB+((z*KAer2/(EtopA-z))^(nB/nA)));
    yer5=EtopBer*((z*KA/(EtopA-z))^(nB/nA))/(KB+((z*KA/(EtopA-z))^(nB/nA)));
    yer6=EtopBer2*((z*KA/(EtopA-z))^(nB/nA))/(KB+((z*KA/(EtopA-z))^(nB/nA)));
    yer7=EtopB*((z*KA/(EtopA-z))^(nB/nA))/(KBer+((z*KA/(EtopA-z))^(nB/nA)));
    yer8=EtopB*((z*KA/(EtopA-z))^(nB/nA))/(KBer2+((z*KA/(EtopA-z))^(nB/nA)));
    yer9=EtopB*((z*KA/(EtopA-z))^(nB/nAer))/(KB+((z*KA/(EtopA-z))^(nB/nAer)));
    yer10=EtopB*((z*KA/(EtopA-z))^(nB/nAer2))/(KB+((z*KA/(EtopA-z))^(nB/nAer2)));
    yer11=EtopB*((z*KA/(EtopA-z))^(nBer/nA))/(KB+((z*KA/(EtopA-z))^(nBer/nA)));
    yer12=EtopB*((z*KA/(EtopA-z))^(nBer2/nA))/(KB+((z*KA/(EtopA-z))^(nBer2/nA)));

    
    yvaler2(end+1)=yer2;
    yvaler3(end+1)=yer3;
    yvaler4(end+1)=yer4;
    yvaler5(end+1)=yer5;
    yvaler6(end+1)=yer6;
    yvaler7(end+1)=yer7;
    yvaler8(end+1)=yer8;
    yvaler9(end+1)=yer9;
    yvaler10(end+1)=yer10;
    yvaler11(end+1)=yer11;
    yvaler12(end+1)=yer12;
    
    
    yer2ref=EtopBref*((z*KAref/(EtopAer2ref-z))^(nBref/nAref))/(KBref+((z*KAref/(EtopAer2ref-z))^(nBref/nAref)));
    yer3ref=EtopBref*((z*KAerref/(EtopAref-z))^(nBref/nAref))/(KBref+((z*KAerref/(EtopAref-z))^(nBref/nAref)));
    yer4ref=EtopBref*((z*KAer2ref/(EtopAref-z))^(nBref/nAref))/(KBref+((z*KAer2ref/(EtopAref-z))^(nBref/nAref)));
    yer5ref=EtopBerref*((z*KAref/(EtopAref-z))^(nBref/nAref))/(KBref+((z*KAref/(EtopAref-z))^(nBref/nAref)));
    yer6ref=EtopBer2ref*((z*KAref/(EtopAref-z))^(nBref/nAref))/(KBref+((z*KAref/(EtopAref-z))^(nBref/nAref)));
    yer7ref=EtopBref*((z*KAref/(EtopAref-z))^(nBref/nAref))/(KBerref+((z*KAref/(EtopAref-z))^(nBref/nAref)));
    yer8ref=EtopBref*((z*KAref/(EtopAref-z))^(nBref/nAref))/(KBer2ref+((z*KAref/(EtopAref-z))^(nBref/nAref)));
    yer9ref=EtopBref*((z*KAref/(EtopAref-z))^(nBref/nAerref))/(KBref+((z*KAref/(EtopAref-z))^(nBref/nAerref)));
    yer10ref=EtopBref*((z*KAref/(EtopAref-z))^(nBref/nAer2ref))/(KBref+((z*KAref/(EtopAref-z))^(nBref/nAer2ref)));
    yer11ref=EtopBref*((z*KAref/(EtopAref-z))^(nBerref/nAref))/(KBref+((z*KAref/(EtopAref-z))^(nBerref/nAref)));
    yer12ref=EtopBref*((z*KAref/(EtopAref-z))^(nBer2ref/nAref))/(KBref+((z*KAref/(EtopAref-z))^(nBer2ref/nAref)));
    
    
    yvaler2ref(end+1)=yer2ref;
    yvaler3ref(end+1)=yer3ref;
    yvaler4ref(end+1)=yer4ref;
    yvaler5ref(end+1)=yer5ref;
    yvaler6ref(end+1)=yer6ref;
    yvaler7ref(end+1)=yer7ref;
    yvaler8ref(end+1)=yer8ref;
    yvaler9ref(end+1)=yer9ref;
    yvaler10ref(end+1)=yer10ref;
    yvaler11ref(end+1)=yer11ref;
    yvaler12ref(end+1)=yer12ref;
    
    
    end
   
    
    tareaer1=[];tareaer2=[];tareaer3=[];tareaer4=[];tareaer5=[];tareaer6=[];tareaer7=[];tareaer8=[];tareaer9=[];tareaer10=[];tareaer11=[];tareaer12=[];
    tareaer1ref=[];tareaer2ref=[];tareaer3ref=[];tareaer4ref=[];tareaer5ref=[];tareaer6ref=[];tareaer7ref=[];tareaer8ref=[];tareaer9ref=[];tareaer10ref=[];tareaer11ref=[];tareaer12ref=[];
    if EtopAer<EtopAerref
        
        for m=1:l3-1
            minnieer1=min(yvaler1(m),yvaler1(m+1));
            areaer1=0.001*minnieer1+0.5*0.001*abs(yvaler1(m+1)-yvaler1(m));
            tareaer1(end+1)=areaer1;
        end
        for m=1:l3-1
            minnieer1ref=min(yvaler1ref(m),yvaler1ref(m+1));
            areaer1ref=0.001*minnieer1ref+0.5*0.001*abs(yvaler1ref(m+1)-yvaler1ref(m));
            tareaer1ref(end+1)=areaer1ref;
        end
    end
    
    
   
    
    for m=1:l2-1
        minnieer2=min(yvaler2(m),yvaler2(m+1));
        minnieer3=min(yvaler3(m),yvaler3(m+1));
        minnieer4=min(yvaler4(m),yvaler4(m+1));
        minnieer5=min(yvaler5(m),yvaler5(m+1));
        minnieer6=min(yvaler6(m),yvaler6(m+1));
        minnieer7=min(yvaler7(m),yvaler7(m+1));
        minnieer8=min(yvaler8(m),yvaler8(m+1));
        minnieer9=min(yvaler9(m),yvaler9(m+1));
        minnieer10=min(yvaler10(m),yvaler10(m+1));
        minnieer11=min(yvaler11(m),yvaler11(m+1));
        minnieer12=min(yvaler12(m),yvaler12(m+1));
      
        areaer2=0.001*minnieer2+0.5*0.001*abs(yvaler2(m+1)-yvaler2(m));
        areaer3=0.001*minnieer3+0.5*0.001*abs(yvaler3(m+1)-yvaler3(m));
        areaer4=0.001*minnieer4+0.5*0.001*abs(yvaler4(m+1)-yvaler4(m));
        areaer5=0.001*minnieer5+0.5*0.001*abs(yvaler5(m+1)-yvaler5(m));
        areaer6=0.001*minnieer6+0.5*0.001*abs(yvaler6(m+1)-yvaler6(m));
        areaer7=0.001*minnieer7+0.5*0.001*abs(yvaler7(m+1)-yvaler7(m));
        areaer8=0.001*minnieer8+0.5*0.001*abs(yvaler8(m+1)-yvaler8(m));
        areaer9=0.001*minnieer9+0.5*0.001*abs(yvaler9(m+1)-yvaler9(m));
        areaer10=0.001*minnieer10+0.5*0.001*abs(yvaler10(m+1)-yvaler10(m));
        areaer11=0.001*minnieer11+0.5*0.001*abs(yvaler11(m+1)-yvaler11(m));
        areaer12=0.001*minnieer12+0.5*0.001*abs(yvaler12(m+1)-yvaler12(m));
       
        tareaer2(end+1)=areaer2;
        tareaer3(end+1)=areaer3;
        tareaer4(end+1)=areaer4;
        tareaer5(end+1)=areaer5;
        tareaer6(end+1)=areaer6;
        tareaer7(end+1)=areaer7;
        tareaer8(end+1)=areaer8;
        tareaer9(end+1)=areaer9;
        tareaer10(end+1)=areaer10;
        tareaer11(end+1)=areaer11;
        tareaer12(end+1)=areaer12;
        
        
        minnieer2ref=min(yvaler2ref(m),yvaler2ref(m+1));
        minnieer3ref=min(yvaler3ref(m),yvaler3ref(m+1));
        minnieer4ref=min(yvaler4ref(m),yvaler4ref(m+1));
        minnieer5ref=min(yvaler5ref(m),yvaler5ref(m+1));
        minnieer6ref=min(yvaler6ref(m),yvaler6ref(m+1));
        minnieer7ref=min(yvaler7ref(m),yvaler7ref(m+1));
        minnieer8ref=min(yvaler8ref(m),yvaler8ref(m+1));
        minnieer9ref=min(yvaler9ref(m),yvaler9ref(m+1));
        minnieer10ref=min(yvaler10ref(m),yvaler10ref(m+1));
        minnieer11ref=min(yvaler11ref(m),yvaler11ref(m+1));
        minnieer12ref=min(yvaler12ref(m),yvaler12ref(m+1));
        
       
        areaer2ref=0.001*minnieer2ref+0.5*0.001*abs(yvaler2ref(m+1)-yvaler2ref(m));
        areaer3ref=0.001*minnieer3ref+0.5*0.001*abs(yvaler3ref(m+1)-yvaler3ref(m));
        areaer4ref=0.001*minnieer4ref+0.5*0.001*abs(yvaler4ref(m+1)-yvaler4ref(m));
        areaer5ref=0.001*minnieer5ref+0.5*0.001*abs(yvaler5ref(m+1)-yvaler5ref(m));
        areaer6ref=0.001*minnieer6ref+0.5*0.001*abs(yvaler6ref(m+1)-yvaler6ref(m));
        areaer7ref=0.001*minnieer7ref+0.5*0.001*abs(yvaler7ref(m+1)-yvaler7ref(m));
        areaer8ref=0.001*minnieer8ref+0.5*0.001*abs(yvaler8ref(m+1)-yvaler8ref(m));
        areaer9ref=0.001*minnieer9ref+0.5*0.001*abs(yvaler9ref(m+1)-yvaler9ref(m));
        areaer10ref=0.001*minnieer10ref+0.5*0.001*abs(yvaler10ref(m+1)-yvaler10ref(m));
        areaer11ref=0.001*minnieer11ref+0.5*0.001*abs(yvaler11ref(m+1)-yvaler11ref(m));
        areaer12ref=0.001*minnieer12ref+0.5*0.001*abs(yvaler12ref(m+1)-yvaler12ref(m));
        
       
        tareaer2ref(end+1)=areaer2ref;
        tareaer3ref(end+1)=areaer3ref;
        tareaer4ref(end+1)=areaer4ref;
        tareaer5ref(end+1)=areaer5ref;
        tareaer6ref(end+1)=areaer6ref;
        tareaer7ref(end+1)=areaer7ref;
        tareaer8ref(end+1)=areaer8ref;
        tareaer9ref(end+1)=areaer9ref;
        tareaer10ref(end+1)=areaer10ref;
        tareaer11ref(end+1)=areaer11ref;
        tareaer12ref(end+1)=areaer12ref;
    end
  
    totareaer1=sum(tareaer1);
    totareaer2=sum(tareaer2);
    totareaer3=sum(tareaer3);
    totareaer4=sum(tareaer4);
    totareaer5=sum(tareaer5);
    totareaer6=sum(tareaer6);
    totareaer7=sum(tareaer7);
    totareaer8=sum(tareaer8);
    totareaer9=sum(tareaer9);
    totareaer10=sum(tareaer10);
    totareaer11=sum(tareaer11);
    totareaer12=sum(tareaer12);
    
    totareaer1ref=sum(tareaer1ref);
    totareaer2ref=sum(tareaer2ref);
    totareaer3ref=sum(tareaer3ref);
    totareaer4ref=sum(tareaer4ref);
    totareaer5ref=sum(tareaer5ref);
    totareaer6ref=sum(tareaer6ref);
    totareaer7ref=sum(tareaer7ref);
    totareaer8ref=sum(tareaer8ref);
    totareaer9ref=sum(tareaer9ref);
    totareaer10ref=sum(tareaer10ref);
    totareaer11ref=sum(tareaer11ref);
    totareaer12ref=sum(tareaer12ref);
   
    
%calculate the mean area under each curve

arealigmean = (totareaer1+totareaer2+totareaer3+totareaer4+totareaer5+totareaer6+totareaer7+totareaer8+totareaer9+totareaer10+totareaer11+totareaer12)/12

arearefmean = (totareaer1ref+totareaer2ref+totareaer3ref+totareaer4ref+totareaer5ref+totareaer6ref+totareaer7ref+totareaer8ref+totareaer9ref+totareaer10ref+totareaer11ref+totareaer12ref)/12

%calculate the error for each area under the curve
arealiger = ((totareaer1-arealigmean)^2 +(totareaer2-arealigmean)^2 +(totareaer3-arealigmean)^2 +(totareaer4-arealigmean)^2 +(totareaer5-arealigmean)^2 +(totareaer6-arealigmean)^2 +(totareaer7-arealigmean)^2 +(totareaer8-arealigmean)^2 +(totareaer9-arealigmean)^2 +(totareaer10-arealigmean)^2 +(totareaer11-arealigmean)^2 +(totareaer12-arealigmean)^2)^0.5/2

arearefer = ((totareaer1ref-arearefmean)^2 +(totareaer2ref-arearefmean)^2 +(totareaer3ref-arearefmean)^2 +(totareaer4ref-arearefmean)^2 +(totareaer5ref-arearefmean)^2 +(totareaer6ref-arearefmean)^2 +(totareaer7ref-arearefmean)^2 +(totareaer8ref-arearefmean)^2 +(totareaer9ref-arearefmean)^2 +(totareaer10ref-arearefmean)^2 +(totareaer11ref-arearefmean)^2 +(totareaer12ref-arearefmean)^2)^0.5/2

%calculate kappa
kappa = arealigmean/arearefmean
%calculate kappa errors from errors above
kappaer = ((arearefer/arearefmean)^2+(arealiger/arealigmean)^2)^0.5*kappa

if l>0
    clear kappa
end
end
 