function [outdata] = malkuswheel_nd(varargin);

%% [outdata] = malkuswheel(varargin);
%%
%%  Examples  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  [outdata] = malkuswheel_nd('t', 10, 'unimass');
%%  [outdata] = malkuswheel_nd('t', 1, 'arbmass', massvect, 'R', 0.13);
%%  [outdata] = malkuswheel_nd('pstruct', p, 'randmass');
%%  [outdata] = malkuswheel_nd('t', 20, 'damp', 0.08);
%%  
%%  Input and output arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%   [argID,    (string)  See below. The values which are unambiguous can 
%%    value]    (varies)  be given without the preceeding argID.
%%
%% Valid argument IDs and corresponding values. Values which are
%% unambiguous (marked with '*') can be given without the preceeding argID.
%%
%%  ODE Parameters  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      Individual parameters used by the differential equations, or a
%%      way to pass in the parameter struct
%%  'N'         (integer)  number of cups
%%  'mass'      (real   )  maximum mass of water in a cup (kg)
%%  'sinscale'  (integer)  exponent in dmdt term (sin^x)
%%  'leakrad'   (real   )  radius of hole in cup          (m)
%%  'rc'        (real   )  radius of cup                  (m)
%%  'Rw'        (real   )  radius of wheel                (m)
%%  'phi'       (real   )  angle of elavation of wheel    (deg)
%%  'damp'      (real   )  damping coefficient
%%  'flow'      (real   )  water flow rate                (kg/s)
%%                         Note: garden tap ~1000L/hr, so 1000kg/hr=0.28 kg/s
%%  'pstruct'   (struct )  parameters struct
%%  'lyapunov'  (struct )  outdata struct - run lyapunov experiment
%%
%%  Run Parameters  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      User can define type of initial mass distribution, duration of simulation
%%      initial angular displacement and type of ODE solver
%%  'randmass'  (*string)  random initial mass distribution (default)
%%  'unimass'   (real   )  uniform initial mass distribution
%%  'arbmass'   (vector )  arbitrary initial mass distribution, followed by
%%                         initial distribution of water mass (vector)
%%  't'         (integer)  number of minutes to simulate (default: 20)
%%  'omega'     (real   )  initial angular displacement of wheel
%%  'odetype'   (*string)  type of ODE integrator: ode45 (default), ode15s
%%  'recfig'    (*string)  record figures to disk automatically
%%  'delay'     (real   )  delay coordinate for Poincare Sections
%%
%%  Dimensionless Parameters  (not working presently) %%%%%%%%
%%      These dimensionless parameters are used to derive run parameters
%%  'sigma'     (real   )  prandtl number
%%  'rayleigh'  (real   )  rayleigh number
%%  'K'         (real   )  flow (?)
%%
%%  Output Data-Struct  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  datastruct.p            (struct) structure holding run parameters
%%            .elapsed_time (real  ) time spent simulating water wheel
%%            .rayleigh     (real  ) Rayleigh number
%%            .sigma        (real  ) Prandtl number
%%            .mass         (matrix) matrix of mass distributions as function of time
%%            .omega        (matrix) matrix of angular velocity as function of time
%%            .t            (vector) vector of time steps
%%            .mass_o       (vector) vector of initial mass distribution
%%
%%  Authors  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Erin Miller   (emiller@csss.santafe.edu)
%%  Nicky Grigg   (grigg@csss.santafe.edu)
%%  Aaron Clauset (clauset@csss.santafe.edu)
%%
%%  Santa Fe Institute, Complex Systems Summer School 2003
%%  date created : 06.24.03
%%  last modified: 06.26.03
%%  MATLAB file (not necessarily OCTAVE compatible)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% path(path,'/amd/poptart/root/research/aaron/m_files/');

% default ODE parameter values
p = struct('N'           ,  12,      ...       % number of cups
           'm_max'       ,  0.4,    ...       % max water mass           (kg)
           'large_n'     ,  1000,    ...       % sin-scaling exponent
           'leak_radius' ,  0.0014,  ...       % radius of leak-hole      (m)
           'a'           ,  0,       ...       % area of leak-hole        (m^2)
           'rc'          ,  0.03325, ...       % radius of cup            (m)
           'g'           ,  9.8,     ...       % accel. due to gravity    (m/s^2)
           'Rw'          ,  0.25,    ...       % radius of wheel          (m)
           'rho'         ,  998,     ...       % density of water         (kg/m^3)
           'phi_0'       ,  15/180*pi, ...     % angle of wheel elevation (rad)
           'alpha'       ,  0.1,   ...       % damping constant
           'Q'           ,  0.1081);              % water in-flow rate       (kg/s)
p.a = pi*p.leak_radius^2;

% default run parameter values
rparm = struct('mdist'   , 'randmass', ...
               'mscale'  ,  0.7,       ...
               'duration',  20,        ...     % run duration              (minutes)
               'omega'   ,  0,        ...     % intial angular velocity   (rad/s)
               'figures' ,  1,         ...
               'odetype' , 'ode45',    ...     % type of integration
               'nof3'    ,  1,         ...
               'recfig'  ,  0,         ...
               'lyap'    ,  0,         ...
               'delay'   ,  30);               % Poincare Section delay coord (s)

% parse command-line parameters
i=1; 
while i<=length(varargin), 
  argok = 1; 
  if ischar(varargin{i}), 
    switch varargin{i},
		% ODE parameters
     case 'N',           p.N           = varargin{i+1}; i = i + 1;
     case 'mass',        p.m_max       = varargin{i+1}; i = i + 1;
     case 'sinscale',    p.large_n     = varargin{i+1}; i = i + 1;
     case 'leakrad',     p.leak_radius = varargin{i+1}; i = i + 1;
     case 'rc',          p.rc          = varargin{i+1}; i = i + 1;
     case 'Rw',          p.Rw          = varargin{i+1}; i = i + 1;
     case 'phi',         p.phi_0       = varargin{i+1}; i = i + 1;
     case 'damp',        p.alpha       = varargin{i+1}; i = i + 1;
     case 'flow',        p.Q           = varargin{i+1}; i = i + 1;
     case 'pstruct',     p             = varargin{i+1}; i = i + 1;
     case 'lyapunov',    indata        = varargin{i+1}; i = i + 1;
                         delta         = varargin{i};
                         p             = outdata.p;
                         rparm.lyap    = 1;
		% Run parameters
     case {'randmass','unimass'}, rparm.mdist = varargin{i};
     case 'arbmass',              rparm.mdist = varargin{i};
	                          mdist       = varargin{i+1}; i = i + 1;
     case 't',                         rparm.duration = varargin{i+1}; i = i + 1;
     case 'omega',                     rparm.omega    = varargin{i+1}; i = i + 1;
     case {'ode45','ode15s', 'ode23'}, rparm.odetype  = varargin{i};
     case 'nof3',                      rparm.nof3     = 1;
     case 'nofigs',                    rparm.figures  = 0;
     case 'recfig',                    rparm.recfig   = 1;
     case 'delay',                     rparm.delay    = varargin{i+1}; i = i + 1;

     otherwise argok=0; 
    end
  end
  if ~argok, 
    disp(['(Malkus Waterwheel) Ignoring invalid argument #' num2str(i+1)]); 
  end
  i = i+1; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% store initial clock value (for datastruct.elapsed_time)
start_time = clock;

%% initialize timespace vector
tspan       = [0, rparm.duration*60];   % timespace vector          (seconds)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize mass and cup distribution
if (rparm.lyap==0)
    switch (rparm.mdist)
    case 'arbmass'          %put amount of water in cups according to mdist
        m = mdist;
    case 'unimass'          %put same amount of water in each cup
        m = ones(p.N,1) * rparm.mscale * p.m_max;
    otherwise % 'randmass'  %put random amount of water in cups
        m = rand(p.N,1) * rparm.mscale * p.m_max;
    end
    p.thetadist = transpose  (linspace(0, 2*pi, p.N+1) );
    p.thetadist = p.thetadist(1:end-1);
else              
    m           = indata.mass(end,:);  % grab final mass distribution
    m(1)        = m(1) + delta;        % add delta mass
    p.thetadist = indata.theta(end,:); % grab final theta distribution
    rparm.omega = indata.omega(end);   % grab final omega distribution
end;
outdata.mass_o = m;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% calculate non-dimensional parameters

p.T=sqrt(p.Rw)/sqrt(p.g*sin(p.phi_0));
p.D1=sqrt(2*p.a^2*p.rho*p.Rw)/sqrt(pi*p.rc^2*p.m_max*sin(p.phi_0));
p.D2=sqrt(p.Q^2*p.Rw)/sqrt(p.m_max^2*p.g*sin(p.phi_0));
p.D3=p.alpha/p.m_max/sqrt(p.Rw^3*p.g*sin(p.phi_0));

m_nd=m/p.m_max;
rparm.omega_nd=rparm.omega*p.T;
tspan_nd=tspan/p.T;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% run ODE solver
%ode_options = odeset('reltol',5e-4);
ode_options = [];
eval(['[t_nd,y_nd]=' rparm.odetype '(@waterwheel_nd, tspan_nd, [m_nd; p.thetadist(1); rparm.omega_nd], ode_options, p);']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% extract ODE results from y-matrix
%% and convert to dimensional quantities
m      = y_nd(:,1:p.N)*p.m_max;
theta1 = y_nd(:,p.N+1);
theta  = ones(length(t_nd),1)*p.thetadist'+theta1*ones(1,p.N);
omega  = y_nd(:,end)/p.T;
t      = t_nd*p.T;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find approximately dodgy dimensionless numbers
mf                   = abs(m(end,:)); % Note: mass occassionally slightly out-of-bounds (integration errors)
outdata.sigma        = p.alpha/(p.a*sum(sqrt(2*p.g*p.rho*mf/pi/p.rc^2))*p.Rw^2);
outdata.rayleigh     = p.Q*pi*p.g*sin(p.phi_0)*p.Rw*sum(mf)^2/(p.alpha*sum(p.a*sqrt(2*p.g*p.rho*mf/pi/p.rc^2)));
outdata.mass         = m;
outdata.omega        = omega;
outdata.theta        = theta;
outdata.elapsed_time = etime(clock,start_time)/60; %time to run script (minutes)
outdata.t            = t;
outdata.p            = p;


fprintf('sigma    = %10.6f\n',outdata.sigma);
fprintf('rayleigh = %10.6f\n',outdata.rayleigh);
fprintf('time     =   %4.2f real minutes for %4.0f simulated minutes\n\n',outdata.elapsed_time, rparm.duration);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot figures of simulation results
if (rparm.figures>0)
 figure(1); clf;

 %% plot cup-mass as a function of time (p.N data series)
 subplot(3,1,1);
 plot(t/60,m);
 ylabel('mass (kg)');
 xlabel('time (minutes)');
 axis([0 rparm.duration 0 ceil(p.m_max*10)/10]);
 for i=1:p.N legstr{i} = int2str(i); end
% legend(legstr,-1);

 %% plot angular velocity as a function of time
 subplot(3,1,2);
 plot(t/60,omega);
 ylabel('omega (rad/s)');
 xlabel('time (minutes)');
 axis([0 rparm.duration 1.1*min(omega) 1.1*max(omega)]);

 %% plot last 2 minutes of omega-vs-time data
 if (t(end)/60>3)
  subplot(3,1,3);
  a = find(t>(t(end)-180));
  plot( t(a(1):end)/60,omega(a(1):end));
  ylabel('omega (rad/s)');
  xlabel('last 3 minutes');
  axis([t(a(1))/60 rparm.duration 1.1*min(omega) 1.1*max(omega)]);
 end;
 if (rparm.recfig>0) print -depsc -f1 malkuswheel_timeseries.eps; end;

 %% plot total mass versus angular velocity (phase diagram?)
 figure(2); clf;
 bound = floor(size(m,1).*0.5);
 plot(omega(end-bound:end),sum(m(end-bound:end,:),2))
 xlabel('omega (rad/s)');
 ylabel('total mass (kg)');
 title(strcat('total mass vs. omega for last: ',int2str(size(omega,1)-bound),' timesteps'));
 legend('data',strcat('sigma = ',num2str(outdata.sigma,5)),strcat('rayleigh = ',num2str(outdata.rayleigh,5)));
 if (rparm.recfig>0) print -depsc -f2 malkuswheel_phaseplot.eps; end;

 %% plot mass distribution as a function of time and cup
 if (rparm.nof3==0)
   figure(3); clf;
   pcolor([transpose(diff(m)); zeros(1,size(m,1)-1)]);
   shading flat; colorbar; colormap(1-gray);
   xlabel('time (minutes)');
   ylabel('cup number');
   if (rparm.recfig>0) print -depsc -f3 malkuswheel_massdiff.eps; end;
 end;
 
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  outdata.p            (struct) structure holding run parameteres
%%         .elapsed_time (real  ) time spent simulating water wheel
%%         .rayleigh     (real  ) Rayleigh number
%%         .sigma        (real  ) Prandtl number
%%         .mass         (matrix) matrix of mass distributions as function of time
%%         .omega        (matrix) matrix of angular velocity as function of time
%%         .theta        (matrix) matrix of angle as function of time
%%         .t            (vector) vector of timestep information
%%         .mass_o       (vector) vector of initial mass distribution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
