function varargout = GKLceta2(dqp, zeta, theta_c, F, mu,varargin)
%ceta = GKLceta2(dqp, zeta, theta_c, F, mu) solves the 
%  differential equation found in T.M. Tran and B.G. Danly "Optimization 
%  of gyroklystron efficiency," Phys. Fluids, 4/1986, pp1274-1281., 
%  equation #44a, using the ODE45 command and then solves for the efficiency 
%  as an average over initial phases 'theta_c'.  The ODE45 command can handle
%  the two vector inputs and is much faster at calculating a solution for each
%  initial phase theta_c in parallel than serially by for-loop.
%
%  The differential equation solved is EQN 44, for nonlinear gyroklystron 
%  analysis (accomodates a prebunched beam):
%
%  dP/d{\zeta} = i({\Delta} - 1 + |P|^2)P + i F exp{-2({\zeta}/{\mu})^2 }
%
%[ceta,P] = GKLceta2(dqp, zeta, theta_c, F, mu) is an alternative usage that 
%  also returns the electron particle momentum evolution, where the size of P 
%  is length(zeta)-by-length(theta_c).  NOTE: if only one output argument is 
%  used, then only the final efficiency value is reported (to work with the 
%  fminsearch algorithm used in GKLcontour2.m).
%
%.... = GKLceta2(dqp, zeta, theta_c, F, mu, Pin) is an optional form that
%  allows the initial normalized momentum to be specified.  Pin is set to 
%  unity by default.  It should be a constant or a vector of size 
%  1-by-length(theta_c).  Velocity spread can be included by initializing 
%  Pin = 1 + spr*rand(1,length(theta_c)); where 'spr' is the amount of 
%  initial momentum spread.
%
%
% Inputs:
%   'dqp'    : contains dqp = [delta, q, psi]  = [scalar detuning factor, 
%      bunching parameter, relative phase ] initial values.
%   'zeta'   : vector containing points to evaluate the solution at.
%   'theta_c': initial phase vector uniformly distributed between [0,2*pi).
%   'F'      : scalar normalized field amplitude parameter.
%   'mu'     : scalar normalized cavity length.
%   'Pin'    : normalized initial momentum (optional).
%
% Outputs:
%   'ceta'   : "complementary eta", that is, ceta = 1 - eta, where eta is 
%     the output efficiency as a function of length 'zeta'.  NOTE: if only one 
%     output argument is used, then only the final efficiency value is reported 
%     (to work with the fminsearch algorithm used in GKLcontour2.m).
%
%   'P'      : electron momentum as a function of 'zeta' and particle number.
%
% Last edit by Colin Joye, 7/22/03 
%  (5/24/04: added 'P' output)
%  (10/21/04: added zeta-dependence to 'ceta' output)
%  (11/3/04: added initial momentum input to function)
tsteps = length(theta_c);
zsteps = length(zeta);
% Initial normalized momentum (EQN 43a) - can add velocity spread here somehow.
if nargin>5,
  pin = varargin{1};
else,
  pin = 1;
end

% Input initial electron phases (EQN 43a):
theta_in = theta_c + dqp(2)*sin(theta_c) - dqp(3); 

% Initial input momentum with phase (EQN 44a):
Pin = pin.*exp(-i * theta_in); 
P   = zeros(zsteps,tsteps);   % output matrix of momentum solutions.

% Solve the DE in parallel for each Pin condition at each point zeta.
% P is thus a matrix of solutions.
%          ODE45( EQN ,TSPAN,Y0,OPT,constants ... )
[zeta,P] = ode45(@myeqn,zeta,Pin,[],F,mu,dqp(1));

% compute complementary efficiency (EQN 45)
%ceta = mean( abs( P(end,:) ).^2, 2 ); % final efficiency value (original)
ceta = mean( abs( P(:,:) ).^2, 2 );    % efficiency vs. zeta (10/21/04)

% assign outputs:
if nargout>1,
  varargout{1} = ceta;
  varargout{2} = P;
else
  varargout{1} = ceta(end);
end

% -- Declare ODE differential equation here --
function Pprime = myeqn(zeta, P, F, mu, delta)
% Equation 44a from Tran & Danly.
Pprime = i*(delta - 1 + abs(P).^2) .* P + ...
    i * F .* exp( -(2 * zeta ./ mu).^2 );
return;

% -------- End GKLceta2.m ------------------