view release on metacpan or  search on metacpan

libnova-0.15.0/src/transform.c  view on Meta::CPAN

	(struct ln_helio_posn *object,  
	struct ln_rect_posn * position)
{
	double sin_e, cos_e;
	double cos_B, sin_B, sin_L, cos_L;
	
	/* ecliptic J2000 */
	sin_e = 0.397777156;
	cos_e = 0.917482062;

	/* calc common values */
	cos_B = cos(ln_deg_to_rad(object->B));
	cos_L = cos(ln_deg_to_rad(object->L));
	sin_B = sin(ln_deg_to_rad(object->B));
	sin_L = sin(ln_deg_to_rad(object->L));
	
	/* equ 37.1 */
	position->X = object->R * cos_L * cos_B;
	position->Y = object->R * (sin_L * cos_B * cos_e - sin_B * sin_e);
	position->Z = object->R * (sin_L * cos_B * sin_e + sin_B * cos_e);
}

/*! \fn void ln_get_hrz_from_equ (struct ln_equ_posn * object, struct ln_lnlat_posn * observer, double JD, struct ln_hrz_posn * position)
* \param object Object coordinates.
* \param observer Observer cordinates.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects equatorial coordinates into horizontal coordinates
* for the given julian day and observers position.
* 
* 0 deg azimuth = south, 90 deg = west.
*/
/* Equ 12.1,12.2 pg 88 
*
* TODO:
* Transform horizontal coordinates to galactic coordinates.
*/

void ln_get_hrz_from_equ (struct ln_equ_posn * object, struct ln_lnlat_posn * observer, double JD, struct ln_hrz_posn * position)
{
	double sidereal;
	
	/* get mean sidereal time in hours*/
	sidereal = ln_get_mean_sidereal_time (JD);
	ln_get_hrz_from_equ_sidereal_time (object, observer, sidereal, position);
}


void ln_get_hrz_from_equ_sidereal_time (struct ln_equ_posn * object, struct ln_lnlat_posn * observer, double sidereal, struct ln_hrz_posn * position)
{
	long double H, ra, latitude, declination, A, Ac, As, h, Z, Zs;

	/* change sidereal_time from hours to radians*/
	sidereal *= 2.0 * M_PI / 24.0;

	/* calculate hour angle of object at observers position */
	ra = ln_deg_to_rad (object->ra);
	H = sidereal + ln_deg_to_rad (observer->lng) - ra;

	/* hence formula 12.5 and 12.6 give */
	/* convert to radians - hour angle, observers latitude, object declination */
	latitude = ln_deg_to_rad (observer->lat);
	declination = ln_deg_to_rad (object->dec);

	/* formula 12.6 *; missuse of A (you have been warned) */
	A = sin (latitude) * sin (declination) + cos (latitude) * cos (declination) * cos (H);
	h = asin (A);

	/* convert back to degrees */
	position->alt = ln_rad_to_deg (h);   

	/* zenith distance, Telescope Control 6.8a */
	Z = acos (A);

	/* is'n there better way to compute that? */
	Zs = sin (Z);

	/* sane check for zenith distance; don't try to divide by 0 */
	if (fabs(Zs) < 1e-5) {
		if (object->dec > 0)
			position->az = 180;
		else
			position->az = 0;
		if ((object->dec > 0 && observer->lat > 0)
		   || (object->dec < 0 && observer->lat < 0))
		  	position->alt = 90;
		else
		  	position->alt = -90;
		return;
	}

	/* formulas TC 6.8d Taff 1991, pp. 2 and 13 - vector transformations */
	As = (cos (declination) * sin (H)) / Zs;
	Ac = (sin (latitude) * cos (declination) * cos (H) - cos (latitude) * sin (declination)) / Zs;

	// don't blom at atan2
	if (Ac == 0 && As == 0) {
	        if (object->dec > 0)
			position->az = 180.0;
		else
			position->az = 0.0;
		return;
	}
	A = atan2 (As, Ac);

	/* convert back to degrees */
	position->az = ln_range_degrees(ln_rad_to_deg (A));
}

/*! \fn void ln_get_equ_from_hrz (struct ln_hrz_posn * object, struct ln_lnlat_posn * observer, double JD, struct ln_equ_posn * position)
* \param object Object coordinates.
* \param observer Observer cordinates.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects horizontal coordinates into equatorial coordinates
* for the given julian day and observers position.
*/
void ln_get_equ_from_hrz (struct ln_hrz_posn * object, struct ln_lnlat_posn * observer, double JD, struct ln_equ_posn * position)
{
	long double H, longitude, declination, latitude, A, h, sidereal;

	/* change observer/object position into radians */

	/* object alt/az */
	A = ln_deg_to_rad (object->az);
	h = ln_deg_to_rad (object->alt);

	/* observer long / lat */
	longitude = ln_deg_to_rad (observer->lng);
	latitude = ln_deg_to_rad (observer->lat);

	/* equ on pg89 */
	H = atan2 (sin (A), ( cos(A) * sin (latitude) + tan(h) * cos (latitude)));
	declination = sin(latitude) * sin(h) - cos(latitude) * cos(h) * cos(A);
	declination = asin (declination);

	/* get ra = sidereal - longitude + H and change sidereal to radians*/
	sidereal = ln_get_apparent_sidereal_time(JD);
	sidereal *= 2.0 * M_PI / 24.0;

	position->ra = ln_range_degrees(ln_rad_to_deg (sidereal - H + longitude));
	position->dec = ln_rad_to_deg (declination);
}

/*! \fn void ln_get_equ_from_ecl (struct ln_lnlat_posn * object, double JD, struct ln_equ_posn * position)
* \param object Object coordinates.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects ecliptical coordinates into equatorial coordinates
* for the given julian day.