#S. L. Chuang, Physics of Optoelectronic Devices table K.3
function generalTernaryInterpolation(b1,b2,x){
    return b1*x+b2*(1-x);
}

#S. L. Chuang, Physics of Optoelectronic Devices table K.3
function generalQuaternaryInterpolationType1(b1,b2,b3,b4,x,y){
    return b1*x*y+b2*x*(1-y)+b3*(1-x)*y+b4*(1-x)*(1-y);
}

#S. L. Chuang, Physics of Photonic Devices table C.3 (newer edition of Physics of Optoelectronic Devices)
function generalQuaternaryInterpolationType2(b1,b2,b3,x,y){
    return b1*(1-x-y)+b1*x*(1-y)+b2*x+b3*y;
}

#Decompose ternary into binaries according to inteprolation formula from S. L. Chuang, Physics of Optoelectronic Devices table K.3
function decomposeTernary(name){
    binaryBasis = cell(2);
    mqw_material_database;
    
    if(name == "Al(x)Ga(1-x)As"){
        binaryBasis{1} = AlAs;
        binaryBasis{2} = GaAs;
    }else if(name == "Ga(x)In(1-x)As"){
        binaryBasis{1} = GaAs;
        binaryBasis{2} = InAs;
    }else if(name == "Al(x)In(1-x)As"){
        binaryBasis{1} = AlAs;
        binaryBasis{2} = InAs;
    }else if(name == "GaAs(x)P(1-x)"){
        binaryBasis{1} = GaAs;
        binaryBasis{2} = GaP;
    }else{
        assert("decomposeTernary: name unknown.",0);
    }
    return binaryBasis;
}

#Decompose ternary into binaries according to inteprolation formula from S. L. Chuang, Physics of Optoelectronic Devices table K.3
function decomposeQuaternaryType1(name){
    binaryBasis = cell(4);
    mqw_material_database;
    
    if(name == "Ga(x)In(1-x)As(y)P(1-y)"){
        binaryBasis{1} = GaAs;
        binaryBasis{2} = GaP;
        binaryBasis{3} = InAs;
        binaryBasis{4} = InP;
    }else{
        assert("decomposeQuaternaryType1: name unknown.",0);
    }
    return binaryBasis;
}

#Decompose ternary into binaries according to inteprolation formula from S. L. Chuang, Physics of Photonics Devices Table C.3
function decomposeQuaternaryType2(name){
    binaryBasis = cell(3);
    mqw_material_database;
    
    if(name == "In(1-x-y)Al(x)Ga(y)As"){
        binaryBasis{1} = InAs;
        binaryBasis{2} = AlAs;
        binaryBasis{3} = GaAs;
    }else{
        assert("decomposeQuaternaryType2: name unknown.",0);
    }
    return binaryBasis;
}

#Implements binary materials from table K.2 in S. L. Chuang, Physics of Optoelectronic Devices
#dEg_dT_300K is user provided temperature dependence where 300K is reference, if dEg_dT_300K=0 band gap is retruned at 300K regardless of T
function buildBinaryMaterialChuang(name,T,dEg_dT_300K){
    mqw_material_database;
    
    AB = struct;
    if(name == "AlAs"){
        AB = AlAs;
    }else if(name == "GaAs"){
        AB = GaAs;
    }else if(name == "InAs"){
        AB = InAs;
    }else if(name == "InP"){
        AB = InP;
    }else if(name == "GaP"){
        AB = GaP;
    }else{
        assert("buildBinaryMaterialChuang: name unknown.",0);
    }
    
    #Band gap custom temperature dependence
    AB.eg = AB.eg + dEg_dT_300K*(T-300);
    
    return AB;
}


#Implements ternary portion of table K.3 from S. L. Chuang, Physics of Optoelectronic Devices
#dEg_dT_300K is user provided temperature dependence where 300K is reference, if dEg_dT_300K=0 band gap is retruned at 300K regardless of T
function buildTernaryMaterialChuang(name,T,dEg_dT_300K,x){
       
    binaryBasis = decomposeTernary(name);
    b1 = binaryBasis{1};
    b2 = binaryBasis{2};
    
    #Parameters from the general interpolation formula
    ABC = struct;
    ABC.lc = generalTernaryInterpolation(b1.lc,b2.lc,x);
    ABC.eg = generalTernaryInterpolation(b1.eg,b2.eg,x);
    ABC.vb = generalTernaryInterpolation(b1.vb,b2.vb,x);    
    ABC.gamma1 = generalTernaryInterpolation(b1.gamma1,b2.gamma1,x);
    ABC.gamma2 = generalTernaryInterpolation(b1.gamma2,b2.gamma2,x);
    ABC.gamma3 = generalTernaryInterpolation(b1.gamma3,b2.gamma3,x);
    ABC.ep = generalTernaryInterpolation(b1.ep,b2.ep,x);
    ABC.c11 = generalTernaryInterpolation(b1.c11,b2.c11,x);
    ABC.c12 = generalTernaryInterpolation(b1.c12,b2.c12,x);
    ABC.ac = generalTernaryInterpolation(b1.ac,b2.ac,x);
    ABC.av = generalTernaryInterpolation(b1.av,b2.av,x);
    ABC.b = generalTernaryInterpolation(b1.b,b2.b,x);
    ABC.me = generalTernaryInterpolation(b1.me,b2.me,x);
    ABC.mhh = generalTernaryInterpolation(b1.mhh,b2.mhh,x);
    ABC.mlh = generalTernaryInterpolation(b1.mlh,b2.mlh,x);
    
    #Override parameters with specific formulas
    if(name == "Al(x)Ga(1-x)As"){
        if (x<0.4){
            ABC.eg = 1.424+1.247*x;
        }
        ABC.me = 0.067 + 0.083*x;
        ABC.mhh = 0.5 + 0.29*x;
        ABC.mlh = 0.087 + 0.063*x;
    }else if(name == "Ga(x)In(1-x)As"){
        ABC.eg = 0.324+0.7*x+0.4*x^2;
        ABC.me = 0.025*(1-x) + 0.071*x - 0.0163*x*(1-x);
    }else if(name == "Al(x)In(1-x)As"){
        if ((x > 0.44) and (x < 0.54)){
            ABC.eg = 0.357+2.29*x;
        }else{
            ABC.eg = 0.36+2.35*x+0.24*x^2;
        }
    }
    
    #Band gap custom temperature dependence
    ABC.eg = ABC.eg + dEg_dT_300K*(T-300);
    
    return ABC;
}

#Implements quaternary portion of table K.3 from S. L. Chuang, Physics of Optoelectronic Devices
#dEg_dT_300K is user provided temperature dependence where 300K is reference, if dEg_dT_300K=0 band gap is retruned at 300K regardless of T
function buildQuaternaryMaterialChuang(name,T,dEg_dT_300K,x,y){
    
    ABCD = struct;
    if(name == "Ga(x)In(1-x)As(y)P(1-y)"){
        binaryBasis = decomposeQuaternaryType1(name);
        b1 = binaryBasis{1};
        b2 = binaryBasis{2};
        b3 = binaryBasis{3};
        b4 = binaryBasis{4};
        
        #Parameters from the general interpolation formula        
        ABCD.lc = generalQuaternaryInterpolationType1(b1.lc,b2.lc,b3.lc,b4.lc,x,y);
        ABCD.eg = generalQuaternaryInterpolationType1(b1.eg,b2.eg,b3.eg,b4.eg,x,y);
        ABCD.vb = generalQuaternaryInterpolationType1(b1.vb,b2.vb,b3.vb,b4.vb,x,y);
        ABCD.gamma1 = generalQuaternaryInterpolationType1(b1.gamma1,b2.gamma1,b3.gamma1,b4.gamma1,x,y);
        ABCD.gamma2 = generalQuaternaryInterpolationType1(b1.gamma2,b2.gamma2,b3.gamma2,b4.gamma2,x,y);
        ABCD.gamma3 = generalQuaternaryInterpolationType1(b1.gamma3,b2.gamma3,b3.gamma3,b4.gamma3,x,y);
        ABCD.ep = generalQuaternaryInterpolationType1(b1.ep,b2.ep,b3.ep,b4.ep,x,y);
        ABCD.c11 = generalQuaternaryInterpolationType1(b1.c11,b2.c11,b3.c11,b4.c11,x,y);
        ABCD.c12 = generalQuaternaryInterpolationType1(b1.c12,b2.c12,b3.c12,b4.c12,x,y);
        ABCD.ac = generalQuaternaryInterpolationType1(b1.ac,b2.ac,b3.ac,b4.ac,x,y);
        ABCD.av = generalQuaternaryInterpolationType1(b1.av,b2.av,b3.av,b4.av,x,y);
        ABCD.b = generalQuaternaryInterpolationType1(b1.b,b2.b,b3.b,b4.b,x,y);
        ABCD.me = generalQuaternaryInterpolationType1(b1.me,b2.me,b3.me,b4.me,x,y);
        ABCD.mhh = generalQuaternaryInterpolationType1(b1.mhh,b2.mhh,b3.mhh,b4.mhh,x,y);
        ABCD.mlh = generalQuaternaryInterpolationType1(b1.mlh,b2.mlh,b3.mlh,b4.mlh,x,y);
    }else if(name == "In(1-x-y)Al(x)Ga(y)As"){
        binaryBasis = decomposeQuaternaryType2(name);
        b1 = binaryBasis{1};
        b2 = binaryBasis{2};
        b3 = binaryBasis{3};
        
        #Parameters from the general interpolation formula
        ABCD.lc = generalQuaternaryInterpolationType2(b1.lc,b2.lc,b3.lc,x,y);
        ABCD.eg = generalQuaternaryInterpolationType2(b1.eg,b2.eg,b3.eg,x,y);
        ABCD.vb = generalQuaternaryInterpolationType2(b1.vb,b2.vb,b3.vb,x,y);
        ABCD.gamma1 = generalQuaternaryInterpolationType2(b1.gamma1,b2.gamma1,b3.gamma1,x,y);
        ABCD.gamma2 = generalQuaternaryInterpolationType2(b1.gamma2,b2.gamma2,b3.gamma2,x,y);
        ABCD.gamma3 = generalQuaternaryInterpolationType2(b1.gamma3,b2.gamma3,b3.gamma3,x,y);
        ABCD.ep = generalQuaternaryInterpolationType2(b1.ep,b2.ep,b3.ep,x,y);
        ABCD.c11 = generalQuaternaryInterpolationType2(b1.c11,b2.c11,b3.c11,x,y);
        ABCD.c12 = generalQuaternaryInterpolationType2(b1.c12,b2.c12,b3.c12,x,y);
        ABCD.ac = generalQuaternaryInterpolationType2(b1.ac,b2.ac,b3.ac,x,y);
        ABCD.av = generalQuaternaryInterpolationType2(b1.av,b2.av,b3.av,x,y);
        ABCD.b = generalQuaternaryInterpolationType2(b1.b,b2.b,b3.b,x,y);
        ABCD.me = generalQuaternaryInterpolationType2(b1.me,b2.me,b3.me,x,y);
        ABCD.mhh = generalQuaternaryInterpolationType2(b1.mhh,b2.mhh,b3.mhh,x,y);
        ABCD.mlh = generalQuaternaryInterpolationType2(b1.mlh,b2.mlh,b3.mlh,x,y);
    }else{
        assert("buildQuaternaryMaterialChuang: name unknown.",0);
    }
    
    #Override parameters with specific formulas
    if(name == "Ga(x)In(1-x)As(y)P(1-y)"){
        if(x==0.1894*y/(0.4184-0.013*y)){#lattice matched to InP
            ABCD.eg = 1.35-0.775*y+0.149*y^2; #this is actually for 298K, but we take it as 300K.
            ABCD.me = 0.08-0.039*y;
            ABCD.mhh = 0.46;
            ABCD.mlh = 0.12-0.099*y+0.03*y^2;
        }else{
            ABCD.eg = 1.35+0.668*x-1.068*y+0.758*x^2+0.078*y^2-0.069*x*y-0.322*x^2*y+0.03*x*y^2;
            ABCD.me = 0.08-0.116*y+0.026*x-0.059*x*y+(0.064-0.02*y)*x^2+(0.06+0.032*x)*y^2; #second edition textbook
            #ABCD.me = 0.08-0.116*x+0.026*y-0.059*x*y+(0.064-0.02*x)*y^2+(0.06+0.032*y)*x^2; #first edition textbook
        }
        ABCD.lc = (5.8688-0.4176*x+0.1896*y+0.0125*x*y)*1e-10;
    }else if(name == "In(1-x-y)Al(x)Ga(y)As"){
        if(y==0.468-0.983*x){#lattice matched to InP
            z=x/0.48;
            ABCD.eg = 0.76+0.49*z+0.2*z^2;
            ABCD.me = 0.0427+0.0328*z;
        }else{
            ABCD.eg = 0.36+2.093*x+0.629*y+0.577*x^2+0.436*y^2+1.013*x*y-2.0*x*y*(1-x-y);
        }        
    }
    
    #Band gap custom temperature dependence
    ABCD.eg = ABCD.eg + dEg_dT_300K*(T-300);
    
    return ABCD;
}

#Implements InP lattice matched quaternary portion of table K.3 from S. L. Chuang, Physics of Optoelectronic Devices
#dEg_dT_300K is user provided temperature dependence where 300K is reference, if dEg_dT_300K=0 band gap is retruned at 300K regardless of T
function buildQuaternaryMaterialInPMatchedChuang(name,T,dEg_dT_300K,x,y){
    
    ABCD = struct;
    if(name == "Ga(x)In(1-x)As(y)P(1-y)"){
        binaryBasis = decomposeQuaternaryType1(name);
        b1 = binaryBasis{1};
        b2 = binaryBasis{2};
        b3 = binaryBasis{3};
        b4 = binaryBasis{4};
        
        #Parameters from the general interpolation formula        
        ABCD.lc = generalQuaternaryInterpolationType1(b1.lc,b2.lc,b3.lc,b4.lc,x,y);
        ABCD.eg = generalQuaternaryInterpolationType1(b1.eg,b2.eg,b3.eg,b4.eg,x,y);
        ABCD.vb = generalQuaternaryInterpolationType1(b1.vb,b2.vb,b3.vb,b4.vb,x,y);
        ABCD.gamma1 = generalQuaternaryInterpolationType1(b1.gamma1,b2.gamma1,b3.gamma1,b4.gamma1,x,y);
        ABCD.gamma2 = generalQuaternaryInterpolationType1(b1.gamma2,b2.gamma2,b3.gamma2,b4.gamma2,x,y);
        ABCD.gamma3 = generalQuaternaryInterpolationType1(b1.gamma3,b2.gamma3,b3.gamma3,b4.gamma3,x,y);
        ABCD.ep = generalQuaternaryInterpolationType1(b1.ep,b2.ep,b3.ep,b4.ep,x,y);
        ABCD.c11 = generalQuaternaryInterpolationType1(b1.c11,b2.c11,b3.c11,b4.c11,x,y);
        ABCD.c12 = generalQuaternaryInterpolationType1(b1.c12,b2.c12,b3.c12,b4.c12,x,y);
        ABCD.ac = generalQuaternaryInterpolationType1(b1.ac,b2.ac,b3.ac,b4.ac,x,y);
        ABCD.av = generalQuaternaryInterpolationType1(b1.av,b2.av,b3.av,b4.av,x,y);
        ABCD.b = generalQuaternaryInterpolationType1(b1.b,b2.b,b3.b,b4.b,x,y);
        ABCD.me = generalQuaternaryInterpolationType1(b1.me,b2.me,b3.me,b4.me,x,y);
        ABCD.mhh = generalQuaternaryInterpolationType1(b1.mhh,b2.mhh,b3.mhh,b4.mhh,x,y);
        ABCD.mlh = generalQuaternaryInterpolationType1(b1.mlh,b2.mlh,b3.mlh,b4.mlh,x,y);
    }else if(name == "In(1-x-y)Al(x)Ga(y)As"){
        binaryBasis = decomposeQuaternaryType2(name);
        b1 = binaryBasis{1};
        b2 = binaryBasis{2};
        b3 = binaryBasis{3};
        
        #Parameters from the general interpolation formula
        ABCD.lc = generalQuaternaryInterpolationType2(b1.lc,b2.lc,b3.lc,x,y);
        ABCD.eg = generalQuaternaryInterpolationType2(b1.eg,b2.eg,b3.eg,x,y);
        ABCD.vb = generalQuaternaryInterpolationType2(b1.vb,b2.vb,b3.vb,x,y);
        ABCD.gamma1 = generalQuaternaryInterpolationType2(b1.gamma1,b2.gamma1,b3.gamma1,x,y);
        ABCD.gamma2 = generalQuaternaryInterpolationType2(b1.gamma2,b2.gamma2,b3.gamma2,x,y);
        ABCD.gamma3 = generalQuaternaryInterpolationType2(b1.gamma3,b2.gamma3,b3.gamma3,x,y);
        ABCD.ep = generalQuaternaryInterpolationType2(b1.ep,b2.ep,b3.ep,x,y);
        ABCD.c11 = generalQuaternaryInterpolationType2(b1.c11,b2.c11,b3.c11,x,y);
        ABCD.c12 = generalQuaternaryInterpolationType2(b1.c12,b2.c12,b3.c12,x,y);
        ABCD.ac = generalQuaternaryInterpolationType2(b1.ac,b2.ac,b3.ac,x,y);
        ABCD.av = generalQuaternaryInterpolationType2(b1.av,b2.av,b3.av,x,y);
        ABCD.b = generalQuaternaryInterpolationType2(b1.b,b2.b,b3.b,x,y);
        ABCD.me = generalQuaternaryInterpolationType2(b1.me,b2.me,b3.me,x,y);
        ABCD.mhh = generalQuaternaryInterpolationType2(b1.mhh,b2.mhh,b3.mhh,x,y);
        ABCD.mlh = generalQuaternaryInterpolationType2(b1.mlh,b2.mlh,b3.mlh,x,y);
    }else{
        assert("buildQuaternaryMaterialChuang: name unknown.",0);
    }
    
    #Override parameters with specific formulas
    if(name == "Ga(x)In(1-x)As(y)P(1-y)"){
        ABCD.eg = 1.35-0.775*y+0.149*y^2; #this is actually for 298K, but we take it as 300K.
        ABCD.me = 0.08-0.039*y;
        ABCD.mhh = 0.46;
        ABCD.mlh = 0.12-0.099*y+0.03*y^2;
    }else if(name == "In(1-x-y)Al(x)Ga(y)As"){
        z=x/0.48;
        ABCD.eg = 0.76+0.49*z+0.2*z^2;
        ABCD.me = 0.0427+0.0328*z;
    }
    ABCD.lc = 5.8688e-10;
    
    #Band gap custom temperature dependence
    ABCD.eg = ABCD.eg + dEg_dT_300K*(T-300);
    
    return ABCD;
}

#Bulk effective DOS
function effectiveDOS3D(m,T){
    kb=1.380649e-23;
    return 2*(2*pi*m*9.1e-31*kb*T/h^2)^1.5*1e6; #m^{-3}
}

#Bulk intrinsic density
function ni3D(me,mh,T,Eg){
    kb=1.380649e-23;
    Nc = effectiveDOS3D(me,T);
    Nv = effectiveDOS3D(mh,T);
    return sqrt(Nc*Nv)*exp(-Eg*e/2/kb/T);
}

#Bulk intrinsic work function
function wf3D(me,mh,T,Eg,vb){
    kb=1.380649e-23;
    return vb-Eg/2.0-kb*T/2.0/e*log(effectiveDOS3D(me,T)/effectiveDOS3D(mh,T));
}