#!/usr/bin/perl
#
# Compute the radial velocity of gas which sits at a radius "R"
#   from the center of the galaxy; the Sun is located at radius "R0".
#   The Sun orbits with speed V0, while the gas orbits at a speed
#   V = some function of R.
#
# They are located as follows:
#
#                      Sun
#                       | \
#                       |  \ 
#              theta ---|>  \     here, angle theta is approx 30 degrees
#                       |    A         from Sun toward point A
#                       |   /   
#                       |  /
#                       | /
#                       |/
#                    MW center
#
# MWR 9/26/2025

use POSIX;

$debug = 1;
$MYPI = 3.14159265;


# distance of Sun from center of MW (kpc)
$R0 = 8.0;
# distance of the gas from center of MW (kpc)
$min_R = 5.0;
$max_R = 11.0;
$d_R = 1.0;
# circular velocity at Sun's distance from center of MW (km/sec)
$V0 = 220.0;



# iterate over a range of distances from the center of the MW
for ($R = $min_R; $R <= $max_R; $R += $d_R) {

		  for ($theta_deg = -180; $theta_deg <= 180; $theta_deg += 5.0) {

			 $theta_rad = $theta_deg*($MYPI/180.0);

			 $sun_vel = $V0;
			 $other_vel = velocity_at_r($R);

			 $xx = $R0 * sin($theta_rad);
			 $yy = ( ($other_vel/$R) - ($sun_vel/$R0) );
			 $vrad = $xx*$yy;
			 
printf "#   theta_deg %6.1f  sine %5.2f \n", $theta_deg, sin($theta_rad);
			 printf " R %6.2f V %6.1f  theta_deg %6.1f  theta_rad %6.2f  vrad  %8.3f \n",
					$R, $other_vel, $theta_deg, $theta_rad, $vrad;

		  }

   printf "\n\n";
}

exit 0;


################################################################
# PROCEDURE: velocity_at_r
#
# DESCRIPTION: Compute the circular orbital velocity of gas
#              at a distance "r" from the center of the MW.
#
# USAGE:     velocity_at_r     rad
#
#            where
#                     rad          is r = distance of gas from center of MW
#
# RETURNS:
#            velocity of gas at distance "r", in km/sec
#
sub velocity_at_r {
  my($rad);
  my($theta_0, $theta, $vel, $fixed_vel);

  $rad = $_[0];

  my($theta_0) = $V0/$R0;
 
  $theta = $theta_0*sqrt($rad/$R0);

  $vel = $theta*$rad;

printf "# theta_0 %8.4f  theta %8.4f  vel %6.1f \n", 
$theta_0, $theta, $vel;


  return($vel);

}