%% very simplified model 10
%an input output model, does not care about the states
% like #7 but uses fsolve instead of solve
function [y,sym] = very_Simple_model10(u,par)

%% part1 calculate x3 and x4
z = u(1);
P1 = u(2);
Prt = u(3);
Prb = u(4);


%% riser part of model

%alternative, probably better solution of this problem, it has an
%analytical solution, nice
rho_g2 = Prt * par.M_G / (par.R * par.T2);
rho_Mix2_av = (Prb - Prt)/(par.g*par.L2);

V2_l = (rho_Mix2_av-rho_g2)/(par.rho_L-rho_g2)*par.V2;
V2_g = par.V2 - V2_l;

x3_ss = rho_g2*V2_g;
x4_ss = par.rho_L*V2_l;



%% pipeline part of model
%use fsolve for this part of the system
rho_g1 = par.M_G * P1 / (par.R * par.T2);

h = newtonRaphson(0.008,par,u,x4_ss,rho_g2);
%h1=0.009398682123451;

%%

Ag0 = par.A2 * (max((par.hc - h),0)/par.hc)^2;

Al0 = par.A2 - Ag0;

alpha_l2_t = 2 * x4_ss / (par.V2*par.rho_L) - Al0/par.A1;

alpha_l2_t_m = (alpha_l2_t* par.rho_L) / ((alpha_l2_t * par.rho_L)+(1-alpha_l2_t)*rho_g2);

rho_2t = alpha_l2_t * par.rho_L + (1-alpha_l2_t) * rho_g2;

w = par.K_pc * z * sqrt ( rho_2t * (Prt - par.P0-2e3));

wg1 = par.K_g * Ag0 * sqrt(rho_g1 *(P1 - Prb));
wl1 = par.K_o * Al0 * sqrt(par.rho_L * (P1 - Prb + par.rho_L*par.g*h));

wldiff = alpha_l2_t_m*w - wl1;
wgdiff = (1-alpha_l2_t_m)*w - wg1;

wldiffRatio = wldiff/wl1;
wgdiffRatio = wgdiff/wg1;

wl2 = alpha_l2_t_m*w;
wg2 = (1-alpha_l2_t_m)*w;


Ql = wl1/par.rho_L;
Qg = wg1/par.Rog_n;

y(1) = Qg * 60 * 1000;
y(2) = Ql * 60 * 1000;  

sym(1) =  x3_ss;
sym(2) =  x4_ss;
sym(3) =  h;
sym(4) =  Ag0;
sym(5) =  Al0;
end







function x = newtonRaphson(x0,par,u,x4_ss,rho_g2)
    xk = x0;
    xprev = xk-1; % to make sure it iterates at least one time
    i=0;
    [eq,diff] = equation_Jacobi(xk,par,u,x4_ss,rho_g2);
     eq = real(eq);
     diff = real(diff);
    while ( abs(xk-xprev)>1e-16 && i<1e4 && diff>1e-4)
        xprev = xk;
        temp = eq/diff;
        xk = xprev + 0.1*temp;
        [eq,diff] = equation_Jacobi(xk,par,u,x4_ss,rho_g2);
        eq = real(eq);
        diff = real(diff);
        i=i+1;
    end
    if(i==1e4)
        x=x0;
    end
    x=xk; 
end


function [eq,diff] = equation_Jacobi(h,par,u,x4_ss,rho_g2)
    z = u(1);
    P1 = u(2);
    Prt = u(3);
    Prb = u(4);
    
    %% write equations here
    %% calc differentials of somethings with respect to h
    Ag = par.A1 * (max((par.hc - max(h,0)),0)/par.hc)^2;
    Al= par.A1 - Ag;  
    
    if h>0
        diffAg_h = 2 * par.A1 * (max((par.hc - h),0)/par.hc) * (1/par.hc);
    else
        h = 0;
        diffAg_h = 0;
    end
    diffAl_h = -diffAg_h;
    
    %% midflows-outflow-innflows, differentiated with respect to h
    rho_g1 = par.M_G * P1 / (par.R * par.T2);

    wG = par.K_g * Ag * sqrt(rho_g1 *(max(P1-Prb,0)));
    wL = par.K_o * Al * sqrt(par.rho_L*(max((P1 - Prb + par.rho_L*par.g*h),0)));
    
    diffwG_h = diffAg_h * par.K_g * sqrt(rho_g1 *(max(P1-Prb,0)));
    diffwL_h = diffAl_h * par.K_o * sqrt(par.rho_L * (max((P1 - Prb + par.rho_L*par.g*h),0))) + (1/2) * par.K_o * Al /(sqrt(par.rho_L * (max((P1 - Prb + par.rho_L*par.g*h),0)))) * par.rho_L; 
    
    %% densities and mass fractions differentiated with respect to h
    alpha_l2_t = 2 * x4_ss / (par.V2*par.rho_L) - Al/par.A1;
    alpha_l2_t_m = (alpha_l2_t* par.rho_L) / ((alpha_l2_t * par.rho_L)+(1-alpha_l2_t)*rho_g2);
    rho_2t = alpha_l2_t * par.rho_L + (1-alpha_l2_t) * rho_g2;
    
    diffalpha_l2_t_h = -diffAl_h/par.A1;
    diffalpha_l2_t_m_h = diffalpha_l2_t_h * ( par.rho_L/((alpha_l2_t * par.rho_L)+(1-alpha_l2_t)*rho_g2) - (par.rho_L - rho_g2)/((alpha_l2_t * par.rho_L)+(1-alpha_l2_t)*rho_g2)^2);
    diffrho_2t_h = diffalpha_l2_t_h * (par.rho_L - rho_g2);
    
    %% total flow differentiated with respect to h
    w = par.K_pc * z * sqrt ( rho_2t * (Prt - par.P0 - 2e3));
    
    diffw_h = par.K_pc * z / sqrt(rho_2t * (Prt - par.P0 - 2e3)) * diffrho_2t_h*(Prt - par.P0 - 2e3);
    
    %% differentiate my equations which are set to zero with respect to h
    eq2 = w*alpha_l2_t_m - wL;
    eq3 = w*(1-alpha_l2_t_m) - wG;  
    
    diffeq2_h = diffw_h*alpha_l2_t_m + diffalpha_l2_t_m_h*w - diffwL_h;
    diffeq3_h = diffw_h*(1-alpha_l2_t_m) - diffalpha_l2_t_m_h*w - diffwG_h;
    
    eq = eq2-eq3;
    diff = diffeq2_h-diffeq3_h;
end