#!/usr/bin/perl
#
# We want to generate random velocities of gas particles so that they
#   have a Maxwellian distribution.  Our method:
#
#     1. calculate the distribution function N(v)
#     2. normalize N(v) 
#     3. pick a total number of atoms -- Ntot
#     4. for all speeds between 0 and vmax
#            calc number of atoms with this speed
#            calc doppler shift in wavelength for these atoms
#            print number of atoms and doppler shift
#

# MWR 2/19/2004

$debug = 1;

# some physical constants (MKS)
$pi = 3.14159;
$c = 2.9979e8;
$kb = 1.38E-23;
$mh = 1.66E-27;
$gasr = 8.31;
# for hydrogen atoms, molar mass (kg/mole)
$molar_mass = 0.001;

# central wavelength  (Angstroms)
$lambda = 6563;

# we'll use this as delta-v in calculating distribution function N(v)
#    units are m/s
$delta_v = 100.0;

# all velocities lie within this range (m/s)
$min_v = 0.0;
$max_v = 100000.0;

# total number of atoms
$ntot = 1e5;

# put values from sliders here ...
$temp = 6000.0;

# Given temperature $temp and mass $mass, calculate the distribution function.
#    as we go, we'll determine the cumulative distribution function.
$sum = 0.0;
$index = 0;
for ($v = $min_v; $v <= $max_v; $v += $delta_v) {
  $xx = 4.0*$pi*$v*$v;
  $yy = ($molar_mass/(2.0*$pi*$gasr*$temp))**1.5;
  $zz = ($molar_mass*$v*$v)/(2.0*$gasr*$temp);
  $zz = exp(-$zz);

  $pv = $xx*$yy*$zz;

  $sum += $pv;
  $probv[$index] = $pv;
  $velocity[$index] = $v;
  $index++;

  if ($debug > 1) {
    printf " %7.1f  %12.7e  %12.7e  %5d \n", $v, $pv, $sum, 
              $index-1;
  }
}
$array_size = $index;


# We now normalize the cumulative array, so that its values run
#    from 0.0 to 1.0.
for ($index = 0; $index < $array_size; $index++) {
  $probv[$index] /= $sum;
  if ($debug > 1) {
    printf " %5d %9.4e %12.7e\n", $index, $velocity[$index], $probv[$index];
  }
}



# Now, let's generate some shifts in wavelength due to these velocities
for ($index = 0; $index < $array_size; $index++) {

  $num_atoms = $ntot*$probv[$index];
  $v = $velocity[$index];

  # we need to come up with a random distribution of directions
  #   relative to the line of sight.  We can do this by
  #   picking two random angles (uniformly distributed
  #   between 0 and pi radians) and calculating the fraction
  #   of velocity in the direction of the observer: it will
  #   be  cos(i1)*cos(i2)
  #
if (0 == 1) {
  $v = $velocity[$index];
  for ($j = 0; $j < $num_atoms; $j++) {
    $i1 = rand(1)*$pi;
    $i2 = rand(1)*2.0*$pi;
    $frac = cos($i1)*cos($i2);
    $this_v = $frac*$v;
    $delta_lambda = $lambda*($this_v/$c);
    $new_lambda = $lambda + $delta_lambda;

    printf "%6d %6d  %13.6e %13.7e %13.7e \n", 
         $index, $j, $this_v, $delta_lambda, $new_lambda;
  }
 }
  
  if ($debug > -1) {
    printf " %6d  %12.5e %9.6e  %14.7f  %14.7e  %14.7e %14.7e \n", 
         $index, $probv[$index], $num_atoms, $v, $delta_lambda, 
         $new_lambda;
  }
  
}


exit 0;

