ast_coord.c

/******************************************************************************
AST_COORD: Convert spherical coordinates to new system.

This procedure converts the longitude-latitude coordinates (a1, b1)
of a point on a sphere into corresponding coordinates (a2, b2) in a
different coordinate system that is specified by the coordinates of its
origin (ao, bo).  The range of a2 will be from -pi to pi.

double	ao, bo		# Origin of new coordinates (radians)
double	ap, bp		# Pole of new coordinates (radians)
double	a1, b1		# Coordinates to be converted (radians)
double	a2, b2		# Converted coordinates (radians)

******************************************************************************/

#include <math.h>

int ast_coord(double ao, double bo, double ap, double bp, double a1, double b1,
	      double *a2, double *b2) {

  double sao=0.0,cao=0.0,sbo=0.0,cbo=0.0,sbp=0.0,cbp=0.0;
  double x=0.0,y=0.0,z=0.0,xp=0.0,yp=0.0,zp=0.0,temp=0.0;

  x=cos(a1)*cos(b1); y=sin(a1)*cos(b1); z=sin(b1);
  xp=cos(ap)*cos(bp); yp=sin(ap)*cos(bp); zp=sin(bp);

  /* Rotate the origin about z */
  sao=sin(ao); cao=cos(ao); sbo=sin(bo); cbo=cos(bo);
  temp=-xp*sao+yp*cao; xp=xp*cao+yp*sao; yp=temp; temp=-x*sao+y*cao;
  x=x*cao+y*sao; y=temp;
  
  /* Rotate the origin about y */
  temp=-xp*sbo+zp*cbo; xp=xp*cbo+zp*sbo; zp=temp; temp=-x*sbo+z*cbo;
  x=x*cbo+z*sbo; z=temp;

  /* Rotate pole around x */
  sbp=zp; cbp=yp; temp=y*cbp+z*sbp; y=y*sbp-z*cbp; z=temp;
  
  /* Final angular coordinates */
  *a2=atan2(y,x); *b2=asin(z);

  return 1;

}