function [Z,mp,ms,lambda]=nma_rocket_mutliStageSolutionLagrange(lagrange,Ve,epsilon,mL)
%design for a multi-stage rocket.
 
%function [Z,mp,ms,lambda]=nma_rocket_mutliStageSolutionLagrange(lagrange,Ve,mL)
%
% design for a multi-stage rocket.
%
% Finds the mp, ms, lambda, and the Z for each rocket stage based on
% optimization based on using the lagrange multiplier method.
%
%INPUT:
%  lagrange: the optimal lagrange multiplier. use the function
%            nma_rocket_getLagrangeMultiplier to determine.
%
%  Ve: A vector whose elements are the effective exhaust velosity in km/sec.
%      Ve = Isp*g
%      Watch out. Ve must be in km/sec.
%      so for an Isp=300 seconds, Ve= 300*9.8066/1000
%      
%  epsilon: A vector whose elements are the epsilon for each rocket stage.
%           This is the structural mass ratio defined as (ms/(mp+ms)) where
%           ms=structural mass, mp=fuel mass
%
%  mL : The payload. in kg
%
%OUTPUT:
%  Z: A vector of the Z's. One element per stage. Z is defined as 
%     m0/mf  where m0 is the initial totoal mass of the rocket at the
%     start of the stage, and mf is the final rocket mass at the end of
%     the stage. (which will be the same as m0 but with mp removed).
%
%  mp: A vector of the minumum mp needed for each stage. mp is the mass of
%      the fuel. in kg.
%
%  ms: A vector of the ms needed for each stage. ms is the mass of
%      the structure of the stage. in kg.
%
% lambda: A vector of the lambda value for each rocket stage.
%      Lambda is defined as the payload ratio = mL/(mp+ms)
%
%EXAMPLE USAGE:
%  This example below solves the 3 stage rocket shown in the above text
%  on page 96.
%
%  First find the lagrange  multiplier:
%
%  deltaV=10.420; Ve=[3.7 3.7 4.2]; epsilon=[0.08 0.09 0.1]; 
%  lagrange = nma_rocket_getLagrangeMultiplier(deltaV,Ve,epsilon)
%
%  Next find the mimumum required fuel for each stage
%  mL=1000;
%
% [Z,mp,ms,lambda]=nma_rocket_mutliStageSolutionLagrange(lagrange,Ve,epsilon,mL)
%
% Z =
% 
%           2.38346597494207
%           2.11863642217073
%           2.87025221091155
% 
% 
% mp =
% 
%           13573.7257938552
%           4557.18875204095
%           2360.85067749024
% 
% 
% ms =
% 
%           1180.32398207437
%           450.710975476578
%            262.31674194336
%            
% lambda =
% 
%         0.0677780009683472
%          0.199684509357321
%          0.381218519485852
         
nStages = length(epsilon);
Z  = zeros(nStages,1);
mp = zeros(nStages,1);
ms = zeros(nStages,1);
lambda = zeros(nStages,1);

for(i = 1:nStages)
    Z(i)  = (lagrange*Ve(i) - 1 )/(lagrange*Ve(i)*epsilon(i));    
end


massOfUpperStage=0;
for(i = nStages:-1:1)
    m     = (1-Z(i))*(mL+massOfUpperStage) / (Z(i)*epsilon(i) - 1);
    if(m<=0)
        error('Failed to find solution. Mass is 0 or negative');
    end
    mp(i) = m-m*epsilon(i);
    ms(i) = m-mp(i);    
    lambda(i) = mL/(mp(i)+ms(i));    
    massOfUpperStage = massOfUpperStage + m;
end