clear all; clc; clear global;
addpath("C:\Program Files (x86)\REFPROP")
Substance = 'CO2';

% lamda = 0.2903;
% PR = 3.1207;

% Initial guess for the parameters: [lamda, PR]
initial_guess = [0.35, 2.9207];

lamda = initial_guess(1);
PR = initial_guess(2);

% Bounds for the parameters: [lamda, PR]
lb = [0, 1.75];               % Lower bounds
ub = [1, 3.3];         % Upper bounds

% Define optimization options       |       for algorithm: sqp or
% interior point or active set
options = optimoptions('fmincon', 'Display', 'iter-detailed', 'Algorithm', 'interior-point','ConstraintTolerance',1e-4);


% Run the optimization
[optimal_params, neg_n_th] = fmincon(@(params) -calculate_n_th(params), initial_guess, [], [], [], [], lb, ub, @(params) nonlinear_constraints(params), options);

% Display results
fprintf('Initial values:\nlamda = %.4f\nPR = %.4f\n\n', ...
        lamda, PR);

fprintf('Optimal parameters:\nlamda = %.4f\nPR = %.4f\n\nMaximized n_th = %.4f\n', ...
        optimal_params(1), optimal_params(2), -neg_n_th);


[n_th, results] = calculate_n_th(optimal_params);

% Save each field of the struct as a separate variable in the workspace
fields = fieldnames(results);   % Get all field names in the struct
for ii = 1:numel(fields)
    varName = fields{ii};        % Current field name
    assignin('base', varName, results.(varName)); % Assign to base workspace
end


fprintf('OPTIMIZATION RESULTS:\nW_t = %.4f MW\nW_c = %.4f MW\nW_rc = %.4f MW\n\nQ_c = %.4f MW\nQ_LTR = %.4f MW\nQ_HTR = %.4f MW\n', ...
        W_t/10^6, W_c/10^6, W_rc/10^6, Q_c/10^6, Q_ltr/10^6, Q_htr/10^6);
    
%%%%%%%%%%%%%%%%%%%%%% DATA IN TABLE %%%%%%%%%%%%%%%%%%%%%%%
num_points = 15; % Define the number of points
point_numbers = {"1" "2is" "2" "3" "4" "5" "6" "7is" "7" "8" "9" "10" "11" "12is" "12"}';  % Point numbers
pressure_MPa = ([P_1 P_2 P_2 P_3 P_4 P_5 P_6 P_7 P_7 P_8 P_9 P_10 P_11 P_12 P_12]./1000)';
pressure_bar = pressure_MPa*10;
temperature_K = ([T_1 T_2_is T_2 T_3 T_4 T_5 T_6 T_7_is T_7 T_8 T_9 T_10 T_11 T_12_is T_12])';
enthalpy = ([h_1 h_2_is h_2 h_3 h_4 h_5 h_6 h_7_is h_7 h_8 h_9 h_10 h_11 h_12_is h_12]./1000)'; %in J/kg
entropy = ([S_1 S_2_is S_2 S_3 S_4 S_5 S_6 S_7_is S_7 S_8 S_9 S_10 S_11 S_12_is S_12]./1000)'; %in kJ/kgK
mass_flow = ([m_1 m_2_is m_2 m_3 m_4 m_5 m_6 m_7_is m_7 m_8 m_9 m_10 m_11 m_12_is m_12])'; %in kg/s


% Create a table
Data = table(point_numbers, pressure_MPa, temperature_K, enthalpy, entropy, mass_flow, 'VariableNames', {'Point Number', 'Pressure (MPa)', 'Temperature (K)', 'Enthalpy (kJ/kg)', 'Entropy (kJ/kgK)', 'Mass Flow (kg/s)'});

% Display the table
disp(Data);

figure;
entropy = ([S_1 S_2_is S_2 S_3 S_4 S_6 S_7_is S_7 S_8 S_9 S_10 S_1 S_12_is S_12]./1000)'; %in kJ/kgK
temperature_K = ([T_1 T_2_is T_2 T_3 T_4 T_6 T_7_is T_7 T_8 T_9 T_10 T_1 T_12_is T_12])';

plot(entropy, temperature_K, '-*')
hold on;

% Adding labels for each point
for i = 1:length(temperature_K)
    text(entropy(i), temperature_K(i), num2str(i), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right')
end

ylabel('Temperature (K)')
xlabel('Entropy (kJ/kgK)')
title('Temperature vs. Entropy Diagram')
grid on;
hold off;


EQUATIONS FOR CHECK
h_9 - (1 - lamda) * h_8 - lamda * h_12
h_3 - (1 - lamda) * (h_8 - h_7) - h_4
h_10 - h_9 - h_2 + h_3



function [n_th, results] = calculate_n_th(params)

lamda = params(1);
PR = params(2);
    

    m_t = 275; %initial guess
    % Struct initialization in order to save all the variables in Workspace
    results = struct;

    % Calculating intermediate variables

    addpath("C:\Program Files (x86)\REFPROP")
    Substance = 'CO2';

    %n_m_T = 0.97; %mechanical efficiency of turbine
    n_is_T = 0.92; %isentropic efficiency of turbine  %%CUSTOM
    n_is_Comp = 0.85; %isentropic efficiency of compressor
    n_is_ReComp = 0.85; %isentropic efficiency of re-compressor
    e_cooler = 0.95; % efficiency of cooler - see below
    e_ltr = 0.93; % efficiency of LTR
    e_htr = 0.95; % efficiency of HTR

    T_min = 35; %Temperature of supercritical CO2 at cooler in C --------- Compressor Inlet Temperature (CIT)
    T_min = T_min + 273.15; %in Kelvin (ready for refpropm)

    T_seawater = 32; % Temperature of sea water in inlet
    T_seawater = T_seawater + 273.15; %in Kelvin (ready for refpropm)
    if T_min < T_seawater + 1
        T_min = T_seawater + 1;
    end

    T_max = 600; %Temperature of supercritical CO2 at turbine entrance in C -------- Turbine Inlet Temperature (TIT)
    T_max = T_max + 273.15; %in Kelvin (ready for refpropm)

    P_min = 74; %Vaccum pressure in compressor in bar SET AS THE MINIMUM CRITICAL POINT PRESSURE (74 bar)
    P_min = P_min*100; %in kPa (ready for refpropm)
    
    P_mc = 250; %CO2 pressure at Main Compressor outlet in bar ---------- Compressor Outlet Pressure (COP)
    P_mc = P_mc*100; %in kPa (ready for refpropm)

    %PR = P_mc / P_min; % MUST be above 1 
    %PR = 1.715;
    P_lowest = P_mc / PR;   %-------- lowest pressure point of the cycle - after the Main Compressor

    %k = 1.0229; % Helper value to calculate h_8

    %lamda = 0.416; % lamda = m_b / m_t; % flow split ratio -> flow directed to re-compressor over the total flow

    Q_in = 60 * 10^6; % reactor power in watts - INITIAL ESTIMATION
    W_net_target = 29.13 * 10^6; % Net power in watts - INITIAL TARGET

    dP_hx = 0; % Pressure difference in heat exchanger in kPa
    dP_HTR_h = 100; % Pressure difference in HTR hot side in kPa
    dP_HTR_c = 150; % Pressure difference in HTR cold side in kPa
    dP_LTR_h = 100; % Pressure difference in LTR hot side in kPa
    dP_LTR_c = 150; % Pressure difference in LTR cold side in kPa
    dP_cooler = 0; % Pressure difference in cooler in kPa

    Cp_g = 1050; %specific heat capacity of gases, stable and irrelevant of Temperasture in J/(kg*K)
    U_htr = 2040; %HTR coeff of heat transport in W/(m^2*K)
    U_ltr = 2150; %LTR coeff of heat transport in W/(m^2*K)



%%%%%%%%%%%%%%%%%%%%%%% Pressure CALCULATIONS PER POINT %%%%%%%%%%%%%%%%%%%%%%%

P_7 = P_mc; % Pmax
P_6 = P_lowest; % after the Main Compressor - derived from the PR
P_5 = P_6 + dP_cooler;
P_4 = P_5;
P_3 = P_4 + dP_LTR_h;
P_2 = P_3 + dP_HTR_h;
P_8 = P_7 - dP_LTR_c;
P_9 = P_8;
P_10 = P_9 - dP_HTR_c;
P_1 = P_10 - dP_hx;
P_11 = P_4;
P_12 = P_9;


%%%%%%%%%%%%%%%%%%%%%%% Temperature, Enthalpy, Entropy CALCULATIONS PER POINT %%%%%%%%%%%%%%%%%%%%%%%

T_1 = T_max;
[h_1,S_1] = refpropm('HS','T',T_1,'P',P_1,Substance); %in J/kg and J/kgK

S_2_is = S_1;
[h_2_is,T_2_is] = refpropm('HT','P',P_2,'S',S_2_is,Substance); %considering isentropic apotonosi s2'=s1

h_2 = h_1 - n_is_T * (h_1 - h_2_is);
[T_2,S_2] = refpropm('TS','P',P_2,'H',h_2,Substance); %in K and J/kgK

T_6 = T_min; % after the cooler
[h_6,S_6] = refpropm('HS','T',T_6,'P',P_6,Substance); %in J/kg and J/kgK

S_7_is = S_6;
[h_7_is,T_7_is] = refpropm('HT','P',P_7,'S',S_7_is,Substance); %considering isentropic apotonosi s7'= s6

h_7 = ((h_7_is - h_6) / n_is_Comp) + h_6;
[T_7,S_7] = refpropm('TS','P',P_7,'H',h_7,Substance); %in K and J/kgK

% h_12 - h_11 calculation based on wanted NET power
% ((W_net_target/m_t) - (h_1 - h_2) + (1 - lamda)*(h_7 - h_6)) / (-lamda) = (h_12 - h_11) 
%h12_h11 = ((W_net_target/m_t) - (h_1 - h_2) + (1 - lamda)*(h_7 - h_6)) / (-lamda);

    T_4 = 325; % initial guess in Kelvin
    [h_4,S_4] = refpropm('HS','T',T_4,'P',P_4,Substance); %in J/kg and J/kgK
    h_4_old = h_4;
 
flag_repeat = 0;
ii = 1;

% Initial guess for the unknowns: [h_3, h_8, h_9, h_10, h_12, h_11, W_net_target, m_t]
initial_guess = [6.67*10^5 6.1*10^5 6.05*10^5 8.76*10^5 5.93*10^5 5.08*10^5 m_t];

for hh = 1:1
%while flag_repeat == 0



    %T_11 = T_5;
    %h_11 = h_5;
    %S_11 = S_5;


% Define the system of equations

% Known values
knowns.e_ltr = e_ltr;
knowns.h_7 = h_7;
knowns.lamda = lamda;
%knowns.h_4 = h_4;
knowns.e_htr = e_htr;
knowns.h_2 = h_2;
%knowns.n_is_ReComp = n_is_ReComp;
%knowns.h_12_is = h_12_is;
knowns.W_net_target = W_net_target;
knowns.h_1 = h_1;
knowns.h_6 = h_6;
knowns.Q_in = Q_in;
%knowns.m_t = m_t;

% Use fsolve to solve the system of equations
options = optimoptions('fsolve', 'Display', 'none');
solution = fsolve(@(vars) system_of_equations(vars, knowns), initial_guess, options);
h_3 = solution(1);
h_8 = solution(2);
h_9 = solution(3);
h_10 = solution(4);
h_12 = solution(5);
h_11 = solution(6);
%W_net_target = solution(7);
m_t = solution(7);



% Guess for the unknowns AFTER FIRST ITERATION: [h_3, h_8, h_9, h_10, h_12, h_11, W_net_target, m_t]
initial_guess = [h_3 h_8 h_9 h_10 h_12 h_11 m_t];


%m_t = Q_in/(h_1 - h_10);

% 
% 
% system = @(x) [
%     e_ltr * (x(1) - h_7) - (x(2) - h_7);                             % Eq1
%     x(1) - h_4 - (1 - lamda) * (x(2) - h_7);                        % Eq2
%     e_htr * (h_2 - x(3)) - (x(4) - x(3));                           % Eq3
%     h_2 - x(1) - (x(4) - x(3));                                     % Eq4
%     x(5) - (x(3) - (1 - lamda) * x(2)) / lamda;                     % Eq5
%     x(6) + (x(5) * n_is_ReComp - h_12_is) / (n_is_ReComp - 1);      % Eq6
% ];
% 
% % Initial guess for fsolve
% x0 = [6.67 6.1 6.05 8.76 5.93 h_4]*10^5; % Generic initial guess
% 
% % Solve the system using fsolve
% options = optimoptions('fsolve', 'Display', 'none'); % Optional: Display iterations
% solution = fsolve(system, x0, options);
% 
% % Display the result
% 
% h_3 = solution(1);
% h_8 = solution(2);
% h_9 = solution(3);
% h_10 = solution(4);
% h_12 = solution(5);
% h_11 = solution(6);

diff_Temp4_old = 0;
    %h_4 = h_11;
    % [T_4,S_4] = refpropm('TS','H',h_4,'P',P_4,Substance); %in J/kg and J/kgK

    % Initialize variables
    tolerance = 10; % Convergence threshold
    max_iterations = 10000; % Maximum number of iterations
    diff_Temp4 = abs(h_11 - h_4_old);
    diff_h12_h11 = abs(h_12 - h_11);
% Check for convergence
    if  abs(diff_Temp4_old - diff_Temp4) < tolerance
        flag_repeat = 1; % Stop the loop
        disp('Convergence achieved.');
    elseif ii >= max_iterations
        warning('Maximum iterations reached without convergence.');
        flag_repeat = 1; % Stop the loop
    else
        % Log current iteration details
        fprintf('Iteration %d: Difference = %.4f\n', ii, diff_Temp4);

        % Update variables for the next iteration
        %T_4 = T_4 + diff_Temp4 * 0.000005;
        T_4 = T_4 + 0.05;
        diff_Temp4_old = diff_Temp4;
        %h_4_old = h_4_old + 10; % Adjust this logic as needed
        ii = ii + 1; % Increment iteration counter
        [h_4,S_4] = refpropm('HS','T',T_4,'P',P_4,Substance); %in J/kg and J/kgK
    end
h_4_old = h_4;
    % Check for divergence or NaN
    if isnan(diff_Temp4) || isinf(diff_Temp4)
        error('Divergence detected: h_11 or h_4_old became NaN or Inf.');
    end


%    h_10 = h_1 - (Q_in / m_t);



% h_9 h_8 h_3 are unknowns 

%h_8 = k * h_12; initial approach - not physicaly correct
%h_9 = (1 - lamda) * h_8 + lamda * h_12;
%h_3 = (1 - lamda) * (h_8 - h_7) + h_4;
%h_10 = h_9 + h_2 - h_3;
%lamda_test = (h_10 - h_2 - h_7 + h_4) / (h_12 - h_7)


end 
h_4 = h_11;  
[T_4,S_4] = refpropm('TS','H',h_4,'P',P_4,Substance); %in J/kg and J/kgK
clear solution system x0 tolerance;

    T_5 = T_4;
    h_5 = h_4;
    S_5 = S_4;

[T_8,S_8] = refpropm('TS','P',P_8,'H',h_8,Substance); %in K and J/kgK
 
[T_9,S_9] = refpropm('TS','P',P_9,'H',h_9,Substance); %in K and J/kgK

[T_3,S_3] = refpropm('TS','P',P_3,'H',h_3,Substance); %in K and J/kgK

[T_10,S_10] = refpropm('TS','P',P_10,'H',h_10,Substance); %in K and J/kgK


[T_12,S_12] = refpropm('TS','P',P_12,'H',h_12,Substance); %in K and J/kgK

[T_11,S_11] = refpropm('TS','P',P_11,'H',h_11,Substance); %in K and J/kgK

S_12_is = S_4;

[h_12_is,T_12_is] = refpropm('HT','P',P_12,'S',S_4,Substance); %considering isentropic apotonosi s12'= s4


%%%%%%%%%%%%%%%%%%%%%%% MASS FLOW CALCULATION %%%%%%%%%%%%%%%%%%%%%%%

% Q_in = m_t * (h_1 - h_10); % Heat added by reactor
%m_t = Q_in / (h_1 - h_10);

% lamda = m_b / m_t;
m_b = lamda * m_t;

% m_t = m_a + m_b
m_a = m_t - m_b;

%%%%%%%%%%%%%%%%%%%%%%% MASS FLOW CALCULATION PER POINT %%%%%%%%%%%%%%%%%%%%%%%
m_1 = m_t;
H_1 = m_1 * h_1;

m_2_is = m_t; 
m_2 = m_t;
H_2_is = m_2_is * h_2;
H_2 = m_2 * h_2;

m_6 = m_a; % after the cooler
H_6 = m_a * h_6;

m_7_is = m_a;
m_7 = m_a;
H_7_is = m_7_is * h_7;
H_7 = m_7 * h_7;

m_5 = m_a;
H_5 = m_5 * h_5;

m_4 = m_t;
H_4 = m_4 * h_4;

m_11 = m_b;
H_11 = m_11 * h_11;

m_12_is = m_b;
m_12 = m_b; 
H_12_is = m_12_is * h_12;
H_12 = m_12 * h_12;

m_8 = m_a;
H_8 = m_8 * h_8;

m_9 = m_t; 
H_9 = m_9 * h_9;

m_3 = m_t; 
H_3 = m_3 * h_3;

m_10 = m_t;
H_10 = m_10 * h_10;




%%%%%%%%%%%%%%%%%%%%%%% POWER CALCULATION %%%%%%%%%%%%%%%%%%%%%%%
W_t = m_t * (h_1-h_2); % mechanical power of turbine in W
W_c = m_a * (h_7 - h_6); % work of main compressor in W
W_rc = m_b * (h_12 - h_11); % work of re-compressor in W
disp("W_turbine = "+W_t/10^6+" MW" + "          W_compressor = "+W_c/10^6+" MW" + "         W_re-compressor = "+W_rc/10^6+" MW")

Q_c = m_a * (h_5 - h_6); % Heat deducted by cooler in W
Q_htr = m_t * (h_2 - h_3); % heat energy of High Temp Recuperator in W
Q_ltr = m_t * (h_3 - h_4); % heat energy of Low Temp Recuperator in W
disp("Q_cooler = "+Q_c/10^6+" MW" + "          Q_HTR = "+Q_htr/10^6+" MW" + "         Q_LTR = "+Q_ltr/10^6+" MW")


Q_in = m_t * (h_1 - h_10);

%%%%%%%%%%%%%%%%%%%%%%% OVERALL SYSTEM EFFICIENCY %%%%%%%%%%%%%%%%%%%%%%%

n_th = (W_t - W_c - W_rc) / (Q_in);




    %%%%%%%%%%%%%%%%%%%%%%% Heat Exchange SYSTEM REQ AREA CALCULATION %%%%%%%%%%%%%%%%%%%%%%%


        % function Cp_a_b = Cp_a_b(h_a,h_b,T_a,T_b)
        %     Cp_a_b = (h_a - h_b) / (T_a - T_b);
        % end

    Cp_9_10 = Cp_a_b(h_9,h_10,T_9,T_10);
    Cp_2_3 = Cp_a_b(h_2,h_3,T_2,T_3);
    Cp_8_7 = Cp_a_b(h_8,h_7,T_8,T_7);
    Cp_3_4 = Cp_a_b(h_3,h_4,T_3,T_4);

    e_htr_calc = Q_htr / (min(m_t*Cp_9_10,m_t*Cp_2_3) * (T_2 - T_9));
    e_ltr_calc = Q_ltr / (min(m_a*Cp_8_7,m_t*Cp_3_4) * (T_3 - T_7));



        % function Dthetam = D8m(Thot_in,Thot_out,Tcold_in,Tcold_out)
        %     %given that the WHR system works like enalaktis antirrois
        %     Dthetam = ((Thot_out-Tcold_in)-(Thot_in-Tcold_out))/(log((Thot_out-Tcold_in)/(Thot_in-Tcold_out)));
        % end

    Dthetam_LTR = D8m(T_3,T_4,T_7,T_8);
    A_LTR = Q_ltr/(U_ltr*Dthetam_LTR); %in meters squared 

    Dthetam_HTR = D8m(T_2,T_3,T_9,T_10);
    A_HTR = Q_htr/(U_htr*Dthetam_HTR); %in meters squared 

    % Q_ev = m_6 * (h_7-h_6);
    % Thot_out_ev = Thot_out_sh - Q_ev/(m_g*Cp_g);
    % Dthetam_ev = D8m(Thot_out_sh,Thot_out_ev,T_6,T_7);
    % A_ev = Q_ev/(U_ev*Dthetam_ev); %in meters squared
    % 
    % Q_ph = m_5 * (h_6-h_5);
    % Thot_out_ph = Thot_out_ev - Q_ph/(m_g*Cp_g);
    % Dthetam_ph = D8m(Thot_out_ev,Thot_out_ph,T_5,T_6);
    % A_ph = Q_ph/(U_ph*Dthetam_ph); %in meters squared

    A_total = A_LTR + A_HTR; %total heat exchanger area


    % Add the results. prefix to every variable so that they get submitted
    % to the results struct
    varList = whos;
    for ii = 1:length(varList)
        varName = varList(ii).name;
        results.(varName) = eval(varName);
    end


clear tolerance system solution options x0;

% Save each field of the struct as a separate variable in the workspace
fields = fieldnames(results);   % Get all field names in the struct
for ii = 1:numel(fields)
    varName = fields{ii};        % Current field name
    assignin('base', varName, results.(varName)); % Assign to base workspace
end


    % Display results
    fprintf('Initial values:\nlamda = %.4f\nPR = %.4f\n\n', ...
        lamda, PR);
    
   fprintf('Optimal parameters:\nlamda = %.4f\nPR = %.4f\nMaximized n_th = %.4f\n', ...
       lamda, PR, n_th);

   fprintf('OPTIMIZATION RESULTS:\nW_t = %.4f MW\nW_c = %.4f MW\nW_rc = %.4f MW\n\nW_net = %.4f MW\n\nQ_c = %.4f MW\nQ_LTR = %.4f MW\nQ_HTR = %.4f MW\n', ...
       W_t/10^6, W_c/10^6, W_rc/10^6, W_t/10^6 - W_c/10^6 - W_rc/10^6, Q_c/10^6, Q_ltr/10^6, Q_htr/10^6);




EQUATIONS FOR CHECK
h_9 - (1 - lamda) * h_8 - lamda * h_12
h_3 - (1 - lamda) * (h_8 - h_7) - h_4
h_10 - h_9 - h_2 + h_3





function Cp_a_b = Cp_a_b(h_a,h_b,T_a,T_b)
     Cp_a_b = (h_a - h_b) / (T_a - T_b);
end

function Dthetam = D8m(Thot_in,Thot_out,Tcold_in,Tcold_out)
     %given that the WHR system works like enalaktis antirrois
     Dthetam = ((Thot_out-Tcold_in)-(Thot_in-Tcold_out))/(log((Thot_out-Tcold_in)/(Thot_in-Tcold_out)));
end


function F = system_of_equations(vars, knowns)
    % Extract unknowns
    h_3 = vars(1);
    h_8 = vars(2);
    h_9 = vars(3);
    h_10 = vars(4);
    h_12 = vars(5);
    h_11 = vars(6);
    %W_net_target = vars(7);
    m_t = vars(7);

    % Extract known values from the input
    e_ltr = knowns.e_ltr;
    h_7 = knowns.h_7;
    lamda = knowns.lamda;
    %h_4 = knowns.h_4;
    e_htr = knowns.e_htr;
    h_2 = knowns.h_2;
    %n_is_ReComp = knowns.n_is_ReComp;
    %h_12_is = knowns.h_12_is;
    W_net_target = knowns.W_net_target;
    h_1 = knowns.h_1;
    h_6 = knowns.h_6;
    Q_in = knowns.Q_in;
    %m_t = knowns.m_t;
    
    % Define the system of equations
    F(1) =  e_ltr * (h_3 - h_7) - (h_8 - h_7);  
    F(2) = h_3 - h_11 - (1 - lamda) * (h_8 - h_7);   %Replaced h_4 with h_11
    F(3) = e_htr * (h_2 - h_9) - (h_10 - h_9); 
    F(4) = h_2 - h_3 - (h_10 - h_9); 
    F(5) = h_12 - (h_9 - (1 - lamda) * h_8) / lamda;  
    %F(6) = h_12 - ((h_12_is - h_11) / n_is_ReComp) - h_11;
    F(6) = h_12 - h_11 - ((W_net_target/m_t) - (h_1 - h_2) + (1 - lamda)*(h_7 - h_6)) / (-lamda);
    F(7) = Q_in - m_t * (h_1 - h_10);
end






end


function [c, ceq] = nonlinear_constraints(params)
    lamda = params(1);
    PR = params(2);
    %m_t = params(3);
    
    % Calculate n_th and intermediate variables
    [n_th, results] = calculate_n_th(params);
    

    % Inequality constraints: c <= 0
    c(1) = results.e_htr_calc - 1;   % e_htr_calc must be <= 1
    c(2) = results.e_ltr_calc - 1;   % e_ltr_calc must be <= 1
    c(3) = - results.e_htr_calc;     % e_htr_calc must be >= 0
    c(4) = - results.e_ltr_calc;     % e_ltr_calc must be >= 0
    c(5) = -results.A_HTR;           % A_HTR must be > 0
    c(6) = -results.A_LTR;           % A_LTR must be > 0
     
    % Inequality constraints to ensure A_HTR and A_LTR are real (imaginary part is zero)
    c(7) = abs(imag(results.A_HTR));      % Imaginary part of A_HTR should be zero
    c(8) = abs(imag(results.A_LTR));      % Imaginary part of A_LTR should be zero
    c(9) = -n_th;            % n_th must be >= 0
    c(10) = n_th - 1;        % n_th must be <= 1

    % Additional inequality constraints for check
    c(11) = -lamda + 0;         % lambda >= 0
    c(12) = lamda - 1;      % lambda <= 1
    c(13) = -PR + 1.75;         % PR >= 1.75
    c(14) = PR - 3.3;     % PR <= 3.3 = 250/74
    %c(15) = -m_t + 100;       % m_t >= 100
    %c(16) = m_t - 10000;         % m_t <= 10000

    % Inequality constraints for A_total.       |         WARNING!!! 
    % : MAY NEED CHANGE IN SELECTED ALGORITHM WHEN IN FORCE
    c(17) = results.A_total - 2*10^4; % A_total should be smaller than 20.000 m^2
    c(18) = -results.A_total;         % A_total should be greater than 0 m^2    

    % c(19) = -results.h_3 + 10000;
    % c(20) = -results.h_8 + 10000;
    % c(21) = -results.h_9 + 10000;
    % c(22) = results.h_3 - results.h_2;
    % c(23) = results.h_8 - results.h_3;
    % c(24) = results.h_9 - results.h_3;

    c(19) = -(results.h_9 - (1-lamda)*results.h_8 - lamda*results.h_12) - 10;
    c(20) = (results.h_9 - (1-lamda)*results.h_8 - lamda*results.h_12) - 10;

    c(21) = -(results.h_3 - (1-lamda)*(results.h_8-results.h_7) - results.h_4) - 10;
    c(22) = (results.h_3 - (1-lamda)*(results.h_8-results.h_7) - results.h_4) - 10;

    c(23) = -(results.h_10 + results.h_3 - results.h_9 - results.h_2) - 10;
    c(24) = (results.h_10 + results.h_3 - results.h_9 - results.h_2) - 10; 
    % Equality constraints
    ceq = [];
    %ceq(1) = results.h_9 - (1-lamda)*results.h_8 - lamda*results.h_12;
    %ceq(2) = results.h_3 - (1-lamda)*(results.h_8-results.h_7) - results.h_4;
    %ceq(3) = results.h_10 + results.h_3 - results.h_9 - results.h_2;
    


end






Run this code and correct any bugs