pkg load optim;

# model -  subsequent binding
# H  + G <--> HG   K1
# HG + G <--> HG2  K2

# global parameters
global K1;
global K2;


# ----------------------------------------
# equilibrium

function err = err_equi(x)
    global K1;
    global K2;
    
    global cH;
    global cG;
    global cHG;
    global cHG2;

    global cG0;
    global cH0;
    
    cHG = x(1);
    if( cHG <= 0.0 )
        cHG = 1e-90;
    end
    cHG2 = x(2);
    if( cHG2 <= 0.0 )
        cHG2 = 1e-90;
    end    
# cH0 = cH + cHG +   cHG2
# cG0 = cG + cHG + 2*cHG2

    cH = cH0 - cHG - cHG2;
    if( cH <= 0.0 )
        cH = 1e-90;
    end
    
    cG = cG0 - cHG - 2*cHG2;
    if( cG <= 0.0 )
        cG = 1e-90;
    end

# H  + G <--> HG   K1 = cHG/(cH*cG)
# HG + G <--> HG2  K2 = cHG2/(cHG*cG)
 
    err(1) = log(K1) - log(cHG) + log(cH) + log(cG);
    err(2) = log(K2) - log(cHG2) + log(cHG) + log(cG);
    
    #err
end


function solve_equilibrium()
    
    global cG0;
    global cH0;
    global cH;
    global cG;
    global cHG;
    global cHG2;
        
    # boundaries
    lb(1) = 0;
    ub(1) = cH0/2;
    x0(1) = cH0/4;
    
    lb(2) = 0;
    ub(2) = cG0/3;
    x0(2) = cG0/8;    
      
    x =  lsqnonlin(@err_equi,x0,lb,ub,optimset('MaxIter',1000));
    cHG  = x(1);
    cHG2 = x(2);
    cH = cH0 - cHG - cHG2;
    cG = cG0 - cHG - 2*cHG2;
end

K1 = 2*1000;
K2 = 1000/2;

cH0 = 0.001;

for cG0=0.0001:0.0002:0.006
    solve_equilibrium();
    printf("%20.8e %20.8e %20.8e %20.8e %20.8e\n",cG0,cH,cG,cHG,cHG2);
end