function [KPr,phi1,OPr,phi2,OTr,phi3]=turbomachinery(w)

%%%%%%%%%%%%%%%%%%%%%%%%%material properties
% Gall-Tough stainless steel, tested at 871 C
SSyield = 131;  %MPa
SSdensity = 7850;      %kg/m^3
SSnu = .29;           %poisson's ratio for steel

% fluid properties
LOXcp = 920;    %J/(kg*K)
LOXg = 1.4;     %ratio of specific heats
LOXrho = 1140;  %kg/m^3
KEROcp = 2100;  %J/(kg*K)
KEROrho = 810;  %kg/m^3
KEROg = 1.11;     %ratio of specific heats

%%%%%%%%%%%%%%%%%%%%%%%parameters
% temperature and pressure at stages (0 if unknown)
% [in tank, pump outlet, turbine inlet, turbine outlet]
LOXP = [1e5 6e6 6e6 3.5e6];
LOXT = [140 0 300 0];
KEROP = [1e5 3.5e6 3.5e6 3.5e6];
KEROT = [300 0 0 0];

% all values are [KEROpump, LOXpump, LOXturbine]
n = [.35 .35 .35];      %efficiency
Tmin = [0 155 155];     %LOX condensation temperature
mdot = [.17 .44 .44];   %mass flow rates kg/s

%for the moment, lets assume impulse turbines (no reaction)
sigma = [0 0 0];
%Blade loading (psi) and flow coefficient (phi) get values in the middle of
%Prof. Spakovszky's range
psi = [2.5 2.5 2.5];
phi = [.75 .75 .75];

%%%%%%%%%%%%%%%%%%%%%%%%%Power Requirements
%Pump power necessary = (mdot*dP)/(rho*n)
KEROpumpPower = (mdot(1) * (KEROP(2)-KEROP(1))) / (KEROrho * n(1));
KEROT(2) = KEROT(1) + KEROpumpPower/(mdot(1)*KEROcp);

LOXpumpPower = (mdot(2) * (LOXP(2)-LOXP(1))) / (LOXrho * n(2));
LOXT(2) = LOXT(1) + LOXpumpPower/(mdot(2)*LOXcp);

%Turbine Power necessary = pump power / n
LOXturbinePower = (LOXpumpPower + KEROpumpPower) / n(3);

%%%%%%%%%%%%%%%%%%%%%%%%%Design Turbomachinery
%convert rpm to rad/s
w = w * 2 * pi / 60

disp('Kerosene Pump')
%definition of blade loading
r = sqrt(KEROpumpPower/(psi(1)*w^2));
%Euler Turbine Equation (find a = radius at inlet/radius at outlet)
a = sqrt(1 + KEROpumpPower/(w^2*r^2));
rl = a*r
%tip speed
u = w*rl
%minimum hub thickness, assuming a thin, uniform, spinning disk
%maximum stress (both radial and tangential): at the center
stress = ((3+SSnu)/8) * SSdensity * w^2 * rl^2
%just as a placeholder, assume a 1/8" (.3175cm) thick disk
mass = SSdensity * pi * rl^2 * .003175
%calculate the absolute velocity of the fluid at the stator exit
Vx = phi(1)*w*r;
[a2, a3, b2, b3] = angles(phi(1), psi(1), sigma(1))
c2 = Vx/sin(pi/2 - a2)

disp('LOX Pump')
%definition of blade loading
r = sqrt(LOXpumpPower/(psi(2)*w^2));
%Euler Turbine Equation (find a = radius at inlet/radius at outlet)
a = sqrt(1 + LOXpumpPower/(w^2*r^2));
rl = a*r
%tip speed
u = w*rl
%minimum hub thickness, assuming a thin, uniform, spinning disk
%maximum stress (both radial and tangential): at the center
stress = ((3+SSnu)/8) * SSdensity * w^2 * rl^2
%just as a placeholder, assume a 1/8" (.3175cm) thick disk
mass = SSdensity * pi * rl^2 * .003175
%calculate the absolute velocity of the fluid at the stator exit
Vx = phi(2)*w*r;
[a2, a3, b2, b3] = angles(phi(2), psi(2), sigma(2))
c2 = Vx/sin(pi/2 - a2)

disp('LOX Turbine')
%definition of blade loading
r = sqrt(LOXturbinePower/(psi(3)*w^2))
%Euler Turbine Equation (find a = radius at inlet/radius at outlet)
a = sqrt(1 + LOXturbinePower/(w^2*r^2))
rl = a*r
%tip speed
u = w*rl'
%minimum hub thickness, assuming a thin, uniform, spinning disk
%maximum stress (both radial and tangential): at the center
stress = ((3+SSnu)/8) * SSdensity * w^2 * rl^2
%just as a placeholder, assume a 1/8" (.3175cm) thick disk
mass = SSdensity * pi * rl^2 * .003175
%calculate the absolute velocity of the fluid at the stator exit
Vx = phi(3)*w*r
[a2, a3, b2, b3] = angles(phi(3), psi(3), sigma(3))
c2 = Vx/sin(pi/2 - a2)
O2density = LOXP(3)/(8.314*32*LOXT(3))
LOXturbineHeight = mdot(3)/(2*pi*rl*O2density*Vx)

%rotational speed only affects overall radius, indirectly
%overall radius dependent on blade loading, indirectly
%a is dependent on blade loading, directly
    %smaller a reduces rotor blade length needed
%stress is dependent on blade loading, indirectly

function [a2, a3, b2, b3] = angles(phi, psi, sigma)
b2 = atan((psi - 2*sigma)/(2*phi));
b3 = atan((psi + 2*sigma)/(2*phi));
a2 = atan(tan(b2) + 1/phi);
a3 = atan(tan(b3) - 1/phi);