#!/bin/bash

# pm2gal [-h] ra dec [pmra pmdec]
# transform equatorial proper motion vector to Galactic frame
# input: ra, dec ... RA and declination in degrees (J2000)
#        pmra  .... proper motion in RA 
#        pmdec .... proper motion in declination 
#output: l, b ... galactic longitude & latitude in degrees
#        if pm_ra,pm_dec are given the output includes 
#        mu_l, mu_b ... proper motion in l and b coordinates
#        The units of mu_l,mu_b follow that of pmra,pmdec
#
#reference: 
# https://arxiv.org/abs/1306.2945 (for proper motion)
# https://scienceworld.wolfram.com/astronomy/GalacticCoordinates.html    

HELP=0 INVALID=0 

while getopts h optval
do
   case $optval in
	h) HELP=1;;
	*) INVALID=1;;
   esac
done

if [ $# -eq 0 ]; then HELP=1; fi

if [ $HELP -eq 1 ]; then
   echo "help: needs to be written up"
   shift $((OPTIND-1))
   exit
fi

if [ $INVALID -eq 1 ]; then exit; fi


if [ $# -ne 4 ]; then
    echo "need four parameters: ra,dec,pmra,pmdec"; exit -2
fi

echo $1,$2
echo $* | \
awk 'BEGIN{pi=3.141592653589793; 
        torad=pi/180;
	alpha0=192.85948;  #Right Ascension of NGP
	delta0=27.1285;	   #declination of NGP
	l0 = 122.93192;	   #Galactic longitude of NGP
	s0=sin(delta0*torad); c0= cos(delta0*torad)
        }
{ ra=$1; dec=$2; pmra=$3; pmdec=$4
  S0=sin((ra-alpha0)*torad); C0=cos((ra-alpha0)*torad);
  C1 = s0*cos(dec*torad)-c0*sin(dec*torad)*C0
  C2 = c0*S0;
  cosb=sqrt(C1^2+C2^2);
  mu_l = (+C1*pmra+C2*pmdec)/cosb
  mu_b = (-C2*pmra+C1*pmdec)/cosb

  
  A = cos(dec*torad)*c0*C0+sin(dec*torad)*s0;
  B = cos(dec*torad)*S0;
  C = sin(dec*torad)*c0-cos(dec*torad)*s0*C0;
  b = atan2(A,sqrt(1-A^2))/torad;
  lp =atan2(B,C)/torad;
  l = l0-lp;
  printf "%f %f %0.5f\t%0.5f\n", l, b, mu_l,mu_b}'