%-------------------------------------------------------------------
% MATLAB routine 
% for tracking the motions of free falling mass shells in  
% Jeans-instable ISM clouds with various density profiles
%
% created: April 2010
% author:  Susanne Hoefner, Uppsala University
%-------------------------------------------------------------------

clear
clf

%-------------------------------------------------------------------
% physical and astronomical constants (SI units)
%-------------------------------------------------------------------

grav   = 6.6732E-11;
boltz  = 1.38062E-23;
pmass  = 1.6726E-27;
m_unit = 1.66E-27;

solm   = 1.989E+30;
syear  = 3.1557E+7;
parsc  = 3.0856E+16;

%-------------------------------------------------------------------
% Jeans criterion and cloud parameters
%-------------------------------------------------------------------

% computing Jeans mass M_J and Jeans length L_J 
% for typical conditions in a dense cloud core

T_test   = 10;                          % temperature [K]
mu_test  = 2.3;                         % mean molecular weight
rho_test = mu_test * m_unit * 1.0e+10;  % density [kg/m^3]

M_J = rho_test * sqrt( ( pi*boltz*T_test/(4*mu_test*m_unit*grav*rho_test) )^3 );
L_J = 2*( M_J / rho_test )^(1/3);

% setting the actual cloud mass and radius (Jeans instable)

M_cloud = M_J * 3.0;      % cloud mass (> M_J for instability)
R_cloud = L_J/2;          % cloud radius 

% setting the density profile (total mass = M_cloud)

i_rho = 0;                % homogeneous sphere
% i_rho = 1;                % density profile ~ 1/r
% i_rho = 2;                % density profile ~ 1/r^2
% i_rho = 3;                % density profile ~ 1/r^3

% computing additional cloud properties for given M_cloud & R_cloud

g_cloud    = grav*M_cloud/R_cloud^2;          % grav. acceleration
u_esc      = sqrt(2.0*grav*M_cloud/R_cloud);  % escape velocity
t_ff_cloud = sqrt(pi^2*R_cloud^3/(8*grav*M_cloud)); % free fall time

%-------------------------------------------------------------------
% varying starting point of track (free falling mass shell)
%-------------------------------------------------------------------

ntr=5; % number of tracks

r0_min = 0.20; % minimum value of r0 [r/R_cloud]
r0_max = 1.0;  % maximum value of r0 [r/R_cloud]

hold on 

for j=1:ntr

  dr0       = (r0_max-r0_min)/(ntr-1);
  r0_val(j) =  r0_max - (j-1)*dr0;      % starting point of track

%-------------------------------------------------------------------
% solving the ODEs (equation of motion) for each track
%-------------------------------------------------------------------

  r_0   = r0_val(j);   % starting point [r/R_cloud]
  u_0   = 0.0;         % initial velocity [u/u_esc]

  % grav. acceleration and t_ff for a track starting at r_0
  % for different radial density (mass) profiles 
  %   (only the mass below r_0 affects the movement of a shell
  %    starting at r_0, and this mass stays constant in the  
  %    co-moving frame of the shell during free fall)

  if i_rho == 3
    g_r0 = g_cloud;
    t_ff_r0 = t_ff_cloud*r_0^(3/2);
  elseif i_rho == 2
    g_r0 = g_cloud *  r_0;      % density ~ 1/r^2
    t_ff_r0 = t_ff_cloud*r_0;
  elseif i_rho == 1
    g_r0 = g_cloud * (r_0)^2;   % density ~ 1/r
    t_ff_r0 = t_ff_cloud*r_0^(1/2);
  else
    g_r0 = g_cloud * (r_0)^3;   % homogeneous sphere
    t_ff_r0 = t_ff_cloud;
  end

  % physical scaling factors for variables 
  %   (dimensionless equations for better numerics)

  t_scale = t_ff_r0;  % Note: varies with r_0 !
  r_scale = R_cloud;
  u_scale = u_esc;
  
  y0(1) = r_0;
  y0(2) = u_0;

  tspan = [0, 0.99999];  
  y0 = y0';

  ode = @(t,y) track_ff(t,y,t_scale,r_scale,u_scale,g_r0);
  [t,y] = ode45(ode, tspan, y0);
  
  y_r = y(:,1);
  y_u = y(:,2);

%-------------------------------------------------------------------
% plotting tracks
%-------------------------------------------------------------------

  plot(t*t_scale/(syear*1.0e+6),y_r*r_scale/parsc,'b-'); 

end

% plot parameters and labels

title ('Free Fall');  
ylabel('Radial distance  [pc]'); 
xlabel('Time  [10^6 yrs]'); 

hold off

%-------------------------------------------------------------------

