Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#include <cmath>
#include <string>
#include <iostream>
#include <sash/HESSArray.hh>
#include <crash/RotatedSystem.hh>
#include <stash/Coordinate.hh>
/**
Convert position angles between Equatorial (icrs) and Galactic
coordinates ( all parameters in degrees ):
Input parameters are the original position, position angle and system:
l = longitude coordinate
b = latitude coordinate
pa = position angle
sys = system (can be "icrs" or "galactic")
eps = step size (in deg) used to construct point 2 in the algorithm
( default eps should be fine if you only need a precision
of eps and are not near a pole. )
Output parameter is the position angle in the other system,
restricted to be in the range 0 to 360 deg:
new_pa
Note: In both systems position angles are defined "east of north",
i.e. the positive b direction has pa = 0 deg and the positive
l direction has pa = 90 deg.
I have tested this function for a few random values back and forth,
it gave identical results to the NED Coordinate Transformation Calculator
at http://nedwww.ipac.caltech.edu/forms/calculator.html
*/
double ConvertPositionAngle( double l, double b,
double pa, std::string sys,
double eps = 1.e-3 ) {
/**
The following algorithm is used:
Create a second position in the original system by
stepping from the first position a small angle eps into the direction
specified by the position angle.
Then transform both positions and measure the position
angle of the second position relative to the first position.
*/
// Convert all angles from deg to rad
const double RAD = M_PI / 180; // rad per deg
// l *= RAD; b *= RAD; pa *= RAD; eps *=RAD;
// Construct second position
// Note: The extra factor of cos(new_b) is necessary such that
double l2 = l + eps * sin( RAD*pa ) / cos( RAD*b );
double b2 = b + eps * cos( RAD*pa );
double dist = sqrt( pow(l2-l,2) + pow(b2-b,2) );
std::cout << "l = " << l << std::endl;
std::cout << "b = " << b << std::endl;
std::cout << "l2 = " << l2 << std::endl;
std::cout << "b2 = " << b2 << std::endl;
std::cout << "dist = " << dist << std::endl;
std::cout << "pa = " << pa << std::endl;
// Construct original and new coordinate systems
Sash::HESSArray *hess = &Sash::HESSArray::GetHESSArray();
Stash::SystemRef* sys_ref = (Stash::SystemRef*)
&Crash::RotatedSystem::GetGalacticSystem();
Stash::SystemRef* new_sys_ref = (Stash::SystemRef*)
&hess->GetRADecJ2000System();
if ( sys == "galactic") {
} else if ( sys == "icrs" ) {
// switch them
Stash::SystemRef* temp = sys_ref;
sys_ref = new_sys_ref;
new_sys_ref = temp;
} else {
std::cerr << "Coordinate system not supported: "
<< sys << std::endl;
return -1;
}
std::cout << "sys_ref = " << **sys_ref << endl;
std::cout << "new_sys_ref = " << **new_sys_ref << endl;
// Construct coordinates
Stash::Lambda lc ( l, Stash::Angle::Degrees );
Stash::Beta bc ( b, Stash::Angle::Degrees );
Stash::Lambda lc2( l2, Stash::Angle::Degrees );
Stash::Beta bc2( b2, Stash::Angle::Degrees );
Stash::Coordinate p( lc, bc, *sys_ref );
Stash::Coordinate p2( lc2, bc2, *sys_ref );
// Convert coordinates to new system
Stash::Coordinate new_p = p.GetCoordinate ( **new_sys_ref );
Stash::Coordinate new_p2 = p2.GetCoordinate( **new_sys_ref );
// Get new position angle vector as the difference
// between the new coordinates
double new_l = new_p.GetLambda().GetDegrees();
double new_l2 = new_p2.GetLambda().GetDegrees();
double new_b = new_p.GetBeta().GetDegrees();
double new_b2 = new_p2.GetBeta().GetDegrees();
// Note: The extra factor of cos(new_b) is necessary to
double delta_l = ( new_l2 - new_l ) * cos( RAD*new_b );
double delta_b = ( new_b2 - new_b );
double new_dist = sqrt( pow(delta_l,2) + pow(delta_b,2) );
// Compute new position angle
double new_pa = 0;
const double almost_zero = eps * 1.e-6; // to prevent dividing by almost zero
if ( abs(delta_b) > almost_zero ) {
new_pa = atan( delta_l / delta_b ) / RAD; // atan returns range -pi to pi
if ( new_pa < 0 ) new_pa += 360; // force range 0 to 360 deg
} else {
if ( delta_l > 0 )
new_pa = 90;
else
new_pa = 270;
}
std::cout << "p = " << p << std::endl;
std::cout << "p2 = " << p2 << std::endl;
std::cout << "new_p = " << new_p << std::endl;
std::cout << "new_p2 = " << new_p2 << std::endl;
std::cout << "new_l = " << new_l << std::endl;
std::cout << "new_b = " << new_b << std::endl;
std::cout << "new_l2 = " << new_l2 << std::endl;
std::cout << "new_b2 = " << new_b2 << std::endl;
std::cout << "delta_l = " << delta_l << std::endl;
std::cout << "delta_b = " << delta_b << std::endl;
std::cout << "new_dist= " << new_dist << std::endl;
std::cout << "new_pa = " << new_pa << std::endl;
return new_pa;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.