#!/usr/bin/perl
#
# a simple numerical simulation of a quantum simple harmonic oscillator
#
# MWR 2/5/2005

# pi
$pi = 3.14159;

# spring constant, N/m
$k = 100.0;
# mass attached to spring, kg
$m = 10.0;
# initial displacement from rest, m
$A = 1.0;
# how many steps per amplitude
$num_step = 1000.0;
# how far out should we go in position (multiples of amplitude A)
$num_amplitudes = 5;

# angular frequency of oscillation, rad/s
$omega = sqrt($k/$m);
# period of oscillation, s
$period = (2.0*$pi)/$omega;
# timestep
$dt = $period/$num_step;

# coeff of exponential wave function, in units of (1.0/h-bar)
$alpha = sqrt($k*$m)/2.0;
# stepsize in position
$dx = $A/$num_step;

$start_x = 0.0 - ($num_amplitudes*$A);
$end_x = $num_amplitudes*$A;

# run through the loop twice: first time, calc total_prob, 
#   second time can calc relative prob normalized and print out
$tot_prob = 0.0;
$iter = 0;
for ($x = $start_x; $x <= $end_x; $x += $dx) {

  # this is ground state of oscillator
  $arg = 0.0 - ($alpha*$x*$x);
  $psi = $A*exp($arg);
  $psi_sq = $psi*$psi;
  $tot_prob += $dx*$psi_sq;
  $iter++;

}

$iter = 0;
for ($x = $start_x; $x <= $end_x; $x += $dx) {

  # this is ground state of oscillator
  $arg = 0.0 - ($alpha*$x*$x);
  $psi = $A*exp($arg);
  $psi_sq = $psi*$psi;
  $rel_prob = ($dx*$psi_sq)/$tot_prob;

  # print out information at this time
  # printf " %8d  t %9.4e  pos %9.4e  vel %9.4e  KE %9.4e  SPE %9.4e  tot_e %9.4e \n",
  #     $iter, $t, $pos, $vel, $ke, $spe, $tot_e;
   printf " %8d  x %9.4e  psi %9.4e  prob %9.4e   rel_prob %9.4e \n",
       $iter, $x, $psi, $psi_sq, $rel_prob;

  $iter++;

}


exit(0);
