function linearization(net_data,model)
    """
    extracts the senstivity coefficients from the non-convex Weymouth equation
    """
    con_nl = model[:w_eq]
    var_nl = [model[:ğ�›‘];model[:Ï†];model[:Îº]]
    # aux functions
    raw_index(var::MOI.VariableIndex) = var.value
    raw_index(con::NonlinearConstraintIndex) = con.value
    # extract variables
    var = all_variables(model)
    Nv = length(var)
    # extract Jacobian structure
    d = NLPEvaluator(model)
    MOI.initialize(d, [:Jac])
    jac_str = MOI.jacobian_structure(d)
    # extract operating point
    opertaing_point = value.(var)
    # evaluate Weymouth eq. Jacobian at the operating point
    V = zeros(length(jac_str))
    MOI.eval_constraint_jacobian(d, V, opertaing_point)
    # prepare the result
    I = [js[1] for js in jac_str] # rows, equations
    J = [js[2] for js in jac_str] # cols, variables
    jac = sparse(I, J, V)
    rows = raw_index.(index.(con_nl))
    cols = raw_index.(index.(var_nl))
    # Jacobian
    Jac = Matrix(jac[rows, cols])
    Jac_Ï€ = Jac[net_data[:E],net_data[:N]]
    Jac_Ï† = Jac[net_data[:E],length(net_data[:N]) .+ net_data[:E]]
    Jac_Îº = Jac[net_data[:E],(length(net_data[:N]) .+ length(net_data[:E])) .+ net_data[:E]]
    # operating point
    ğ�›‘Ì‡ = opertaing_point[raw_index.(index.(var_nl))][net_data[:N]]
    Ï†Ì‡ = opertaing_point[raw_index.(index.(var_nl))][length(net_data[:N]) .+ net_data[:E]]
    ÎºÌ‡ = opertaing_point[raw_index.(index.(var_nl))][(length(net_data[:N]) .+ length(net_data[:E])) .+ net_data[:E]]
    # linearization coefficients
    Ï‚1 = inv(Jac_Ï†) * (Jac_Ï€ * ğ�›‘Ì‡ + Jac_Îº * ÎºÌ‡ + Jac_Ï† * Ï†Ì‡)
    Ï‚2 = -inv(Jac_Ï†) * Jac_Ï€
    Ï‚3 = -inv(Jac_Ï†) * Jac_Îº
    # pressure-related
    Ï‚Ì‚2 = net_data[:A]*Ï‚2
    Ï‚Ì‚3 = net_data[:B] + net_data[:A]*Ï‚3
    Ï‚Ì†2 = remove_col_and_row(Ï‚Ì‚2,net_data[:ref])
    Ï‚Ì†2 = full_matrix(inv(Ï‚Ì†2),net_data[:ref])
    # # flow-related
    Ï‚Ì€2 = Ï‚2*Ï‚Ì†2
    Ï‚Ì€3 = Ï‚2*Ï‚Ì†2*Ï‚Ì‚3 - Ï‚3
    # save linearization results
    lin_res = Dict( :jac => Jac,
                    :ğ�›‘Ì‡ => ğ�›‘Ì‡,
                    :Ï†Ì‡ => Ï†Ì‡,
                    :ÎºÌ‡ => ÎºÌ‡,
                    :Ï‚1 => round.(Ï‚1, digits = 10),
                    :Ï‚2 => round.(Ï‚2, digits = 10),
                    :Ï‚3 => round.(Ï‚3, digits = 10),
                    :Ï‚Ì‚2 => round.(Ï‚Ì‚2, digits = 10),
                    :Ï‚Ì‚3 => round.(Ï‚Ì‚3, digits = 10),
                    :Ï‚Ì†2 => round.(Ï‚Ì†2, digits = 10),
                    :Ï‚Ì€2 => round.(Ï‚Ì€2, digits = 10),
                    :Ï‚Ì€3 => round.(Ï‚Ì€3, digits = 10)
                    )
    maximum(-inv(Jac_Ï†) * Jac_Ï€) >= 1e6 ? @info("most likely you are at Bifurcation point") : NaN
    maximum(-inv(Jac_Ï†) * Jac_Îº) >= 1e6 ? @info("most likely you are at Bifurcation point") : NaN
    # solve linearized gas network optimization
    sol_linearized = linearized_opt(net_data,lin_res)
    # check lienarization quality
    maximum(sol_non_convex[:ğ�›‘] .- sol_linearized[:ğ�›‘]) <= 1 ? @info("linearization successful") : @warn("linearization fails")
    return  lin_res
end
