#!/usr/bin/perl
#
# Given values for position and velocity of a globular cluster
#   from Sun et al., Res. A&A, 23 (2023),
#   compute the mass of MW interior to its position.
#
# MWR 12/1/2025


$debug = 1;
$G = 6.67E-11;
$solar_mass = 1.99E30;

while (1) {

		  # get the input values for position
					 printf " gimme X, Y, Z (kpc) : ";
					 $line = <STDIN>;
					 @words = split(/\s+/, " " . $line);
					 $x = $words[1];
					 $y = $words[2];
					 $z = $words[3];
					 if ($debug > 0) {
						printf " X %6.2f   Y %6.2f   Z %6.2f \n", $x, $y, $z;
					 }

		  # get the input values for velocity
					 printf " gimme U, V, W (km/s) : ";
					 $line = <STDIN>;
					 @words = split(/\s+/, " " . $line);
					 $u = $words[1];
					 $v = $words[2];
					 $w = $words[3];
					 if ($debug > 0) {
						printf " U %6.2f   V %6.2f   W %6.2f \n", $u, $v, $w;
					 }


		  # compute the total velocity (m/s)
		  $total_vel_ms = 1000.0*sqrt($u*$u + $v*$v + $w*$w);

		  # compute the radial distance from center of MW (m)
		  $rad_dist_kpc = sqrt($x*$x + $y*$y + $z*$z);
		  $rad_dist_m = 3.08E19*$rad_dist_kpc;

		  if ($debug > 0) {
			 printf "  rad_dist_m  %10.4e  total_vel_ms %10.4e \n",
						  $rad_dist_m, $total_vel_ms;
		  }


		  # compute the interior mass (solar masses)

		  $int_mass_kg = (  (($total_vel_ms*$total_vel_ms)*$rad_dist_m) / $G );
		  $int_mass_solar = $int_mass_kg / $solar_mass;
		  if ($debug > 0) {
			 printf "  int_mass_kg  %10.4e  int_mass_solar %10.4e \n",
						  $int_mass_kg, $int_mass_solar;
		  }


		  # print one line per object, with radial dist and interior mass
        $int_mass_print = $int_mass_solar / 1.0E10;
		  printf "  dist %10.2f   mass_int_e10 %10.2f  \n", $rad_dist_kpc, $int_mass_print;


}




exit 0;