%% Read data from files
clear all;
close all;
clc;

%Read hola1.txt
fileID = fopen('hola1.txt','r');
formatSpec = '%f %f';
Asize = [2 4];
for i=1:5
     line = fgetl(fileID);
     if line(2) == 'V'
        h1 = str2double(line(20:23)); %Enthalpy
     end
end
A = fscanf(fileID,formatSpec, Asize);
fclose(fileID);

%Read hola2.txt
fileID = fopen('hola2.txt','r');
formatSpec = '%f %f';
Bsize = [2 4];
for i=1:5
     line = fgetl(fileID);
     if line(2) == 'V'
        h2 = str2double(line(20:23));   %Enthalpy
     end
end
B = fscanf(fileID,formatSpec, Bsize);
fclose(fileID);

%Read hola3.txt
fileID = fopen('hola3.txt','r');
formatSpec = '%f %f';
Csize = [2 4];
for i=1:5
     line = fgetl(fileID);
     if line(2) == 'V'
        h3 = str2double(line(20:23));   %Enthalpy
     end
end
C = fscanf(fileID,formatSpec, Csize);
fclose(fileID);


%p# = Array of pressures (bar)
%m# = Array of mass flow (kg/s)
%poly# = Constants of 2. degree polonomial

%hola1
p1 = A(1,1:4);              %(bar)
m1 = A(2,1:4);              %(kg/s)
poly1 = polyfit(p1,m1,2);

%hola2
p2 = B(1,1:4);
m2 = B(2,1:4);
poly2 = polyfit(p2,m2,2);

%hola3
p3 = C(1,1:4);
m3 = C(2,1:4);
poly3 = polyfit(p3,m3,2);


%% 1-1 Find the optimal pressure to maximize the power of the plant

%Function handle for pressure and mass flow.
mass_flow = @(p,poly) poly(1).*p.^2+poly(2).*p+poly(3);

%Variables set to zero
big_w = 0;          %Work
big_m = 0;          %Mass Flow
big_p = 0;          %Pressure
big_q = 0;          %Energy
W_out_ideal = 0;    %Work with efficiency of 100%
Q_in = 0;           %Energy into the turbine
H_in = 0;           %Enthaply into the turbine
H_out = 0;          %Enthalpy into the turbine

%The for loop iterates through pressures of 0 to 50 bar
%and finds the optimal pressure to maximize work from the turbine.
for p = linspace(0,50,800);

    %Finds enthalpy (hf, hg and Hfg) for the given pressure
    hf = XSteam('hL_p',p);
    hg = XSteam('hV_p',p);
    Hfg = hg - hf;

    %Finds steam quality of each flow
    x1 = (h1-hf)/Hfg;
    x2 = (h2-hf)/Hfg;
    x3 = (h3-hf)/Hfg;

    %Finds total mass flow from all three wells
    total_mass_flow = mass_flow(p, poly1)*x1 ...
                    + mass_flow(p, poly2)*x2 ...
                    + mass_flow(p, poly3)*x3;

    %Finds entropy of steam before entering the turbine
    S_in = XSteam('sV_p',p);

    %We use that S_out = S_in because we have an isotropic efficiency of
    %90% and we can take the efficiency into account later, when we calulate
    %the work out of the turbine.
    S_out = S_in;
    Iso_efficiency = 0.9;

    %Pressure out of the turbine is given
    p_out = 0.1;

    %Finds steam quality of output from turbine
    S_g_out =XSteam('sV_p',p_out);
    S_f_out =XSteam('sL_p',p_out);
    S_fg_out = S_g_out - S_f_out;
    x_cont = (S_out - S_f_out)/S_fg_out;

    %Finds enthalpy of output from turbine
    H_out_current = XSteam('h_px',p_out,x_cont);

    %Finds work out of the turbine taking the isotropic efficiency into account
    W_out = total_mass_flow*(hg - H_out_current)*Iso_efficiency;

    %If the work for the current pressure is the most up to this we update all variables
    if(W_out > big_w)
           big_w = W_out;
           big_m = total_mass_flow;
           big_p = p;
           H_in = hg;
           H_out = H_out_current;
           big_q = total_mass_flow * H_in;
           Q_in = total_mass_flow*hg;
           W_out_ideal = total_mass_flow*(hg - H_out);

    end
end

disp(['Holtoppsthrystingurinn ' , num2str(big_p) , ...
    ' bar gefur virkjuninni mesta aflid, sem er: ' , ...
    num2str(big_w) , ' kW.'])


%% 2-2 Mass flow of the cooling water (kaelivatn)

%Finds the actual output of enthalpy
H_out_actual = H_in - Iso_efficiency*(H_in-H_out);

%Temperature of condensed water (after going through the condenser)
T_2 = XSteam('Tsat_p',0.1);
%Enthalpy of condesed water (after going through the condenser)
h_2 = XSteam('hL_p',0.1);

%Temperature of the coling water entering the condenser
T_3 = 20;
%Temperatur of the cooling water exiting the condenser
T_4 = T_2 - 5;
%Enthalpy of the cooling water entering the condenser
h_3 = XSteam('hL_T',T_3);
%Enthalpy of the cooling water exiting the condenser
h_4 = XSteam('hL_T',T_4);

%Mass flow needed of the cooling water
mass_flow_cooling_water = (big_m*(H_out_actual - h_2))/(h_4 - h_3);

disp(['Eimsvalinn tharf ', num2str(mass_flow_cooling_water), ' kg/s af kaelivatni.']);

%% 2-3 Thermal efficiency of the plant

%Work produced in proportion to work supplied to the turbine
Thermal_efficiency = big_w/Q_in;

disp(['Varmafraedileg nytni virkjunarinnar er: ', num2str(Thermal_efficiency)]);

%% 2-4 Second law efficiency of the plant

%Temperature of the high-temperature reservoir (Input into turbine)
T_in = XSteam('Tsat_p',big_p);
%Temperature of the low-temperature reservoir (Output of turbine)
T_out = XSteam('Tsat_p',p_out);

%Maximum efficiency of the two reservoirs
Max_efficiency = (1 - ((T_out + 273)/(T_in + 273)));

%Thermal efficiency in proportion to the maximum efficiency
Second_law_efficiency = Thermal_efficiency/Max_efficiency;

disp(['Annars logmals nytni virkjunarinner er: ', num2str(Second_law_efficiency)])

%% Cost of electric power

format long

%Fixed cost
drill_one_hole = 5000000;
drill_holes = drill_one_hole *4;

electric_system = 100;                      %USD/kW
steam_piping = 250;                         %USD/kW
other_cost = 1500;                          %USD/kW

%Total fixed cost of the project
total_fixed = (electric_system + steam_piping + other_cost)*big_w + drill_holes;

%Variable cost
MaintenanceCost = 0.02;                         %USD/kWh

%Other information about the project
capacity_factor = 0.9;
discount_factor = 0.12;
lifetime = 30;

%Hours in a year
hours_per_year = 24*365;                    %h

%Power of the plant when the capacity factor is considered
power = big_w * capacity_factor;

%Cost of running the plant for one year
CostPerYear = MaintenanceCost * power * hours_per_year;

%Set the total variable cost to zero
total_variable = 0;

%Find present value of variable cost through all 30 years with given
%discount factor
for i = 1:lifetime
    total_variable = total_variable + CostPerYear/((1+discount_factor).^i);
end

%Present value of total cost of the project
total_cost = total_variable + total_fixed;

%Find annual revenue required to reach breakeven point by using present
%value equation of annuity(timabundins sjodstreymis)
Annual_Revenue = (total_cost * discount_factor) / (1 - (1/ (1 + discount_factor).^lifetime));

%Find the pricing of kWh by dividing the annual revenue requried by total
%kWh produced in one year
Price_of_kWh = Annual_Revenue/(power*hours_per_year);

disp(['Raforkuverd sem skilar nullpunkti fyrir verkefnid er: ', num2str(Price_of_kWh)]);