#!/usr/bin/perl
#
# Given values for velocity for satellite galaxies from data given in
#   Battaglia et al., A&A 657, A54 (2022).
#   compute the mass of MW interior to its position.
#
# MWR 12/1/2025


$debug = 0;
$G = 6.67E-11;
$solar_mass = 1.99E30;
$kpc_to_m = 3.08E19;

while (1) {

		  # get the input values for position
					 printf " gimme dist (kpc), v_los (km/s): ";
					 $line = <STDIN>;
					 @words = split(/\s+/, " " . $line);
					 $dist_kpc = $words[1];
					 $v_los = $words[2];
                $dist_m = $dist_kpc * $kpc_to_m;
					 if ($debug > 0) {
						printf "  dist_m  %10.4e  v_los km/s %10.4e \n" , $dist_m, $v_los;
					 }

		  # get the input values for proper motion
					 printf " gimme mu_alpha, mu_delta: ";
					 $line = <STDIN>;
					 @words = split(/\s+/, " " . $line);
					 $mu_alpha = $words[1];
					 $mu_delta = $words[2];
					 if ($debug > 0) {
						printf " mu_alpha  %10.5f  mu_delta %10.5f \n", 
                      $mu_alpha, $mu_delta;
					 }


		  # compute the alpha, delta components of velocity
		  $v_alpha_ms = $dist_m * (($mu_alpha*0.001) / 205265.0) / (3.15E7);
		  $v_delta_ms = $dist_m * (($mu_delta*0.001) / 205265.0) / (3.15E7);
        $v_los_ms = $v_los * 1000.0;
        if ($debug > 0) {
          printf "  v_alpha km/s  %10.2f  v_delta %10.2f  v_los %10.2f \n",
                 $v_alpha_ms*0.001, $v_delta_ms*0.001, $v_los;
        }
		  $total_vel_ms = sqrt($v_alpha_ms*$v_alpha_ms 
             + $v_delta_ms*$v_delta_ms + $$v_los_ms*$v_los_ms);

		  # compute the radial distance from center of MW (m)
        #    assumes distance from center of MW = distance from Sun
		  $rad_dist_kpc = $dist_kpc;
		  $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  v_alpha %10.2f v_dec %10.2f v_los %10.2f  \n", 
            $rad_dist_kpc, $int_mass_print,
            $v_alpha_ms/1000.0, $v_delta_ms/1000.0, $v_los;


}




exit 0;