clear all;
close all;
clc;

%% Taka inn gögn
formatSpec = '%f %f';
size = [2 4];
%hola 1
fileID = fopen('hola1.txt','r');
% þessi loopa á basically að sækja töluna 1900kJ/kg, má taka út og setja
% beint inn töluna?
for i=1:5
     line = fgetl(fileID);
     if line(2) == 'V'
        h1 = str2double(line(20:23)); %Enthalpy
     end
end
hola_1 = fscanf(fileID,formatSpec, size);
fclose(fileID);

%hola 2
fileID = fopen('hola2.txt','r');
for i=1:5
     line = fgetl(fileID);
     if line(2) == 'V'
        h2 = str2double(line(20:23));   %Enthalpy
     end
end
hola_2 = fscanf(fileID,formatSpec, size);
fclose(fileID);

%hola 3
fileID = fopen('hola3.txt','r');
for i=1:5
     line = fgetl(fileID);
     if line(2) == 'V'
        h3 = str2double(line(20:23));   %Enthalpy
     end
end
hola_3 = fscanf(fileID,formatSpec, size);
fclose(fileID);


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

%hola1
p_1 = hola_1(1,1:4);
m_1 = hola_1(2,1:4);
poly_1 = polyfit(p_1,m_1,2);

%hola2
p_2 = hola_2(1,1:4);
m_2 = hola_2(2,1:4);
poly_2 = polyfit(p_2,m_2,2);

%hola3
p_3 = hola_3(1,1:4);
m_3 = hola_3(2,1:4);
poly_3 = polyfit(p_3,m_3,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
    %3.3 í Xsteam.pdf
    hg = XSteam('hV_p',p);
    hf = XSteam('hL_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, poly_1)*x1)+(mass_flow(p, poly_2)*x2)+(mass_flow(p, poly_3)*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;
    Isotropic_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)*Isotropic_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 - Isotropic_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)])

%% Kostnaðarútreikningar

%fjárhagslegar forsendur
kostnadur_per_borun = 5000000;
boranir = 4; %3 af 4 heppnast svo til að fa 3 þa væntum við þess að bora 4
borkostnadur_total = kostnadur_per_borun*boranir;

kostn_raforkukerfi = 100;%USD/kW
kostn_gufulagna = 250;%USD/kW
kostn_annar = 1500;%USD/kW (startkostnaður á ári 0)

%Breytilegur kostnaður
kostn_vidhald = 0.02; %USD/kWh

%Heildar startkostnaður reiknaður miðað við stærð virkjunar
total_start = borkostnadur_total + big_w*(kostn_raforkukerfi + kostn_gufulagna + kostn_annar);



%Other information about the project
nytingarstudull = 0.9;
r = 0.12;%ávöxtunarkrafa 12%
liftimi = 30;
klst_per_yr = 24*365;

%Afl eftir nýtingarstuðul
power = big_w * nytingarstudull;

%Rekstrarkostnaður á ársgrundvelli
rekstrarkostn_per_yr = kostn_vidhald * power * klst_per_yr;



%Find present value of variable cost through all 30 years with given
%discount factor
total_kostnadur = 0;
for i = 0:liftimi
    if (i==0)
        total_kostnadur = total_kostnadur + total_start/((1+r).^i); %discount factor á ári 0 er bara 1
    else
        total_kostnadur = total_kostnadur + rekstrarkostn_per_yr/((1+r).^i);
    end
end


%Annuity jafnan PV = PMT x ((1-(1/(1+r)^n))/r); PMT = arstekjur
ars_tekjur = (total_kostnadur) / ((1 - (1/ (1 + r).^liftimi))/r);

%Find the pricing of kWh by dividing the annual revenue requried by total
%kWh produced in one year
verd_per_kWh = ars_tekjur/(klst_per_yr*power);

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