%function omega=flow2omega(p1,p2,u)

% Does not solve the system of equations.  Maybe using some numeric
% approach, but not actually necessary to solve for this case study. 

%% Main function - for testing
%Parameters for dimentionless exponential compressor flow model (Table B1)
%Decided to keep indexing from thesis (starting from 0)
a3= -0.2035;
a2=  0.6938;
a1= -0.7479;
a0= 0.3722;
b2= -11.0505;
b1= 27.5087;
b0= -10.4700;
c2= 5.5616;
c1= -15.0297;
c0= 5.0098;

%Parameters for sizing of compressor. Table 3.1
rho_a=      1.23;   %kg/m^3  --- %should be calculated
kappa_a=    1.4;    %ratio of specific heats
cp_a=       1005;   %J/(kg K)
R_a=        287;    %J/(kg K)
T_amb=      293.2;  %K
d_c=        0.5;    %m (compressor diameter)

%MVs
%N_c = 13370;           %rpm;  range:9500-15300
%U_c=350;                %tip speed m/s
P_ratio= 2.1;
V_flow=500;             %Volumetric flow %400-1200 m^3/s

% Parameters for simulation
p1=[a0,a1,a2,a3,b0,b1,b2,c0,c1,c2];  %dimentionless compressor map
p2=[rho_a, kappa_a, cp_a, R_a, T_amb, d_c]; %compressor sizing
u=[V_flow, P_ratio];   %manipulated variables
%%

syms Uc M Psi Phi a b c  N_c 
% Extraction of parameters
a0=p1(1);
a1=p1(2);
a2=p1(3);
a3=p1(4);
b0=p1(5);
b1=p1(6);
b2=p1(7);
c0=p1(8);
c1=p1(9);
c2=p1(10);

rho_a=      p2(1);  %kg/m^3  
kappa_a=    p2(2);  %ratio of specific heats
cp_a=       p2(3);  %J/(kg K)
R_a=        p2(4); %J/(kg K)
T_amb=      p2(5);  %K
d_c=        p2(6);  %m (compressor diameter)

%Manipulated variables
P_ratio=    u(2);   %Pressure ratio
Qc=         u(1);   %rpm;  range:9500-15300
%Uc =      u(1);   %tip speed

% Equations

eq1=Uc/sqrt(kappa_a*R_a*T_amb)-M==0;  %Eq. 3.5 (Mach number)
eq2=cp_a*T_amb*(P_ratio^((kappa_a-1)/kappa_a)-1)/(0.5*Uc^2)-Psi==0; %Eq. 3.3 (dimensionless pressure head, psi)
eq3=a3*M^3+a2*M^2+a1*M+a0-a==0;        %Eq. 3.13 
eq4=b2*M^2+b1*M+b0-b==0;               %Eq. 3.16
eq5=c2*M^2+c1*M+c0-c==0;               %Eq. 3.18
eq6=a*(1-exp(Psi^c-b))-Phi==0;         %Equation 3.20 (dimentionless/normalized flow rate) 
eq7=Phi*(pi/4)*d_c^2*Uc-Qc==0;         %Eq. 3.4 (Actual volumetric flow)
eq8=pi*d_c*N_c/60-U_c==0;              %Tip speeed (m/s)


%Solve system of equations within constraints (obtained from compressor
%map)
S=solve(eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, 'ReturnConditions', true); %200<Uc, Uc<500, 0<M, 0.5<Psi, Psi<1.4, 0.05<Phi, Phi<0.13, 


%end