#!/usr/bin/perl
#
# Compute the phase angle for a forced harmonic oscillator
#
# MWR 2/6/2021

use POSIX;
$debug = 1;
$DEGTORAD = (3.14159/180.0);

# mass of weight (kg)
$m = 0.1;
# force constant of spring (N/m)
$k = 7.5;
# resistive force constant (N/m/s)
$b = 0.02;
# natural frequency of system (rad/s)
$omega_0 = sqrt($k / $m);

if ($debug > 0) {
  printf "# mass m = %6.3f kg \n", $m;
  printf "# force const k = %6.3f N/m \n", $k;
  printf "# b = %6.3f N/m/s \n", $b;
  printf "# omega_0 = %6.3f rad/s \n", $omega_0;
}


# range of driving frequency (rad/s)
$min_omega = 0.1;
$max_omega = $omega_0*10.0;
$d_omega = 0.02;

for ($omega = $min_omega; $omega <= $max_omega; $omega += $d_omega) {

  $top = ($b/$m)*$omega;
  $bot = ($omega_0*$omega_0) - ($omega*$omega);

  $arg = $top / $bot;
  $phi = atan2($top, $bot);
  $phi_deg = $phi/$DEGTORAD;

  printf " omega %8.3f  arg %8.3f  phi %8.3f  phi_deg %8.3f \n",
         $omega, $arg, $phi, $phi_deg;


}



exit 0;


