- #1
FiberOptix
- 12
- 0
Hi all, I've written code that takes J2000 equatorial coordinates for 2 objects and calculates (1) the projected separation of between object 1 and object 2 at the distance of object 1, and the Euclidean (3D) distance between objects 1 and 2. The projected separations should be <= the Euclidean distance, but most of the time they aren't. I've checked both formulae several times and am left scratching my head. As a last ditch effort, perhaps someone here can spot the error. Projected separations are calculated via the Haversine formula. See the C/C++ functions below, where I've included the struct:
Code:
typedef struct _galaxy{
char Name[12];
double RA, DE;
double D;
} galaxy;
double separation (galaxy gal1, galaxy gal2)
{
// convert from eq coords to Cartesian
double x1 = gal1.D * cos(gal1.DE * DEG_TO_RAD) * cos(gal1.RA * DEG_TO_RAD);
double x2 = gal2.D * cos(gal2.DE * DEG_TO_RAD) * cos(gal2.RA * DEG_TO_RAD);
double y1 = gal1.D * cos(gal1.DE * DEG_TO_RAD) * sin(gal1.RA * DEG_TO_RAD);
double y2 = gal2.D * cos(gal2.DE * DEG_TO_RAD) * sin(gal2.RA * DEG_TO_RAD);
double z1 = gal1.D * sin(gal1.DE * DEG_TO_RAD);
double z2 = gal2.D * sin(gal2.DE * DEG_TO_RAD);
// compute 3d separation
double out = sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1) );
return out;
}
double proj_separation(galaxy gal1, galaxy gal2)
{
double dDE = ( gal1.DE - gal2.DE ) * DEG_TO_RAD;
double dRA = ( gal1.RA - gal2.RA ) * DEG_TO_RAD;
// calculate projected separation of galaxy 2 at the distance of galaxy 1 (Haversine)
double out = gal1.D * ( 2.0 * asin(
sqrt( ( sin( dDE/2.0 ) * sin( dDE/2.0 ) )
+ ( cos( gal1.DE * DEG_TO_RAD )
* cos( gal2.DE * DEG_TO_RAD )
* sin( dRA / 2.0 )
* sin( dRA / 2.0 )
) ) ) ) ;
return out;
}