[pycalcal] push by enrico.spinielli - coord conversion on 2010-01-31 13:16 GMT
pyca...@googlecode.com
Author: Enrico Spinielli <enrico.s...@gmail.com>
Date: Sun Jan 31 05:16:34 2010
Log: coord conversion
http://code.google.com/p/pycalcal/source/detail?r=511a10d50d
Modified:
/pycalcal.nw
=======================================
--- /pycalcal.nw Wed Jan 27 12:10:30 2010
+++ /pycalcal.nw Sun Jan 31 05:16:34 2010
@@ -3296,111 +3296,87 @@
# Time and Astronomy #
######################
def `ecliptical_from_equatorial(ra, declination, obliquity):
- """Convert equatorial coordinates to ecliptical ones.
+ """Convert equatorial coordinates (in degrees) to ecliptical ones.
'declination' is the declination,
'ra' is the right ascension and
'obliquity' is the obliquity of the ecliptic.
NOTE: if 'apparent' right ascension and declination are used,
then 'true'
obliquity should be input.
"""
+ co = cos_degrees(obliquity)
+ so = sin_degrees(obliquity)
+ sa = sin_degrees(ra)
lon = degrees_from_radians(
- atan2(sin_degrees(ra)*cos_degrees(obliquity) +
- tan_degrees(declination)*sin_degrees(obliquity),
- cos_degrees(ra)))
+ atan2(sa*co + tan_degrees(declination)*so, cos_degrees(ra)))
lat = arcsin_degrees(
- sin_degrees(declination)*cos_degrees(obliquity) -
-
cos_degrees(declination)*sin_degrees(obliquity)*sin_degrees(ra))
- return lon, lat
-
-
-
-# def ecl_to_equ(longitude, latitude, obliquity):
-# """Convert ecliptic to equitorial coordinates.
-
-# [Meeus-1998: equations 13.3, 13.4]
-
-# Parameters:
-# longitude : ecliptic longitude in radians
-# latitude : ecliptic latitude in radians
-# obliquity : obliquity of the ecliptic in radians
-
-# Returns:
-# Right accension in radians
-# Declination in radians
-
-# """
-# cose = cos(obliquity)
-# sine = sin(obliquity)
-# sinl = sin(longitude)
-# ra = modpi2(atan2(sinl * cose - tan(latitude) * sine,
cos(longitude)))
-# dec = asin(sin(latitude) * cose + cos(latitude) * sine * sinl)
-# return ra, dec
+ sin_degrees(declination)*co -
+ cos_degrees(declination)*so*sa)
+ return [lon, lat]
def `equatorial_from_ecliptical(longitude, latitude, obliquity):
- """Convert ecliptical coordinates to equatorial ones.
+ """Convert ecliptical coordinates (in degrees) to equatorial ones.
'longitude' is the ecliptical longitude,
'latitude' is thestarting ecliptical latitude and
'obliquity' is the obliquity of the ecliptic.
NOTE: resuting ra and declination will be referred to the same equinox
as the one of input ecliptical longitude and latitude.
"""
+ co = cos_degrees(obliquity)
+ so = sin_degrees(obliquity)
+ sl = sin_degrees(longitude)
ra = degrees_from_radians(
- atan2(sin_degrees(longitude)*cos_degrees(obliquity) -
- tan_degrees(latitude)*sin_degrees(obliquity),
- cos_degrees(longitude)))
+ atan2(sl*co - tan_degrees(latitude)*so,
+ cos_degrees(longitude)))
dec = arcsin_degrees(
- sin_degrees(latitude)*cos_degrees(obliquity) +
-
cos_degrees(latitude)*sin_degrees(obliquity)*sin_degrees(longitude))
- return ra, dec
+ sin_degrees(latitude)*co +
+ cos_degrees(latitude)*so*sl)
+ return [ra, dec]
-
-# def horizontal_from_equatorial(H, declination):
-# """Convert equitorial to horizontal coordinates.
-# NOTE: azimuth is measured westward from the South.
-# NOTE:
-# This is not a good formula for using near the poles.
-
-# Parameters:
-# H : hour angle in radians
-# decl : declination in radians
-
-# Returns:
-# azimuth in radians
-# altitude in radians
-
-# """
-# cosH = cos(H)
-# sinLat = sin(astrolabe.globals.latitude)
-# cosLat = cos(astrolabe.globals.latitude)
-# A = atan2(sin(H), cosH * sinLat - tan(decl) * cosLat)
-# h = asin(sinLat * sin(decl) + cosLat * cos(decl) * cosH)
-# return A, h
-
-
-
-def `horizontal_from_ecliptical(ra, dec, ha, location):
- """Convert equatorial coordinates to horizontal ones.
- ra is the right ascension, decl is the declination,
- location is the observer's position on Earth,
- hs is the local hour angle (ha=theta - ra; theta is local sidereal
time).
- NOTE: if ra is affected by nutation, then sidereal time too
- must be affected by it.
+def `horizontal_from_equatorial(H, declination, latitude):
+ """Convert equatorial coordinates (in degrees) to horizontal ones.
+ Return azimuth and altitude.
+ 'H' is the hour angle,
+ 'declination' is the declination,
+ 'latitude' is the observer's geographic latitude.
+ NOTE: azimuth is measured westward from the South.
+ NOTE: This is not a good formula for using near the poles.
+ """
+ ch = cos_degrees(H)
+ sl = sin_degrees(latitude)
+ cl = cos_degrees(latitude)
+ A = degrees_from_radians(
+ atan2(sin_degrees(H),
+ ch * sl - tan_degrees(declination) * cl))
+ h = arcsin_degrees(sl * sin_degrees(declination) +
+ cl * cos_degrees(declination) * ch)
+ return [A, h]
+
+def `equatorial_from_horizontal(azimuth, altitude, latitude):
+ """Convert equatorial coordinates (in degrees) to horizontal ones.
+ Return hour angle and declination.
+ 'azimuth' is the azimuth,
+ 'altitude' is the altitude,
+ 'latitude' is the observer's geographic latitude.
+ NOTE: azimuth is measured westward from the South.
"""
- phi = latitude(location)
- az = degrees_from_radians(atan2(sin_degrees(H),
- cos_degrees(H)*sin_degrees(phi) -
- tan_degrees(delta)*cos_degrees(phi)))
- el = arcsin_degrees(sin_degrees(phi)*sin_degrees(delta) +
- cos_degrees(phi)*cos_degrees(delta)*cos_degrees(H))
- return az, el
+ ca = cos_degrees(azimuth)
+ sl = sin_degrees(latitude)
+ cl = cos_degrees(latitude)
+ H = degrees_from_radians(
+ atan2(sin_degrees(azimuth),
+ ca * sl + tan_degrees(altitude) * cl))
+ d = arcsin_degrees(sl * sin_degrees(altitude) -
+ cl * cos_degrees(altitude) * ca)
+ return [H, d]
+
<<astronomical algorithms tests>>=
def `testEquatorialFromEcliptical(self):
# from the values in the Ch 13 Astronomical Algorithms
- lo, la = ecliptical_from_equatorial(28.026183, 116.328942, 23.4392911)
+ lo, la = ecliptical_from_equatorial(116.328942, 28.026183, 23.4392911)
self.assertAlmostEqual(lo, mpf(113.215630), 6)
self.assertAlmostEqual(la, mpf(6.684170), 6)
@@ -3410,6 +3386,17 @@
self.assertAlmostEqual(ra, 116.328942, 5)
self.assertAlmostEqual(de, 28.026183, 6)
+def `testHorizontalFromEquatorial(self):
+ # from the values in the Ch 13 Astronomical Algorithms
+ A, h = horizontal_from_equatorial(64.352133, -6.719892, 38.9214)
+ self.assertAlmostEqual(A, 68.0337, 4)
+ self.assertAlmostEqual(h, 15.1249, 4)
+
+def `testEquatorialFromHorizontal(self):
+ # from the values in the Ch 13 Astronomical Algorithms
+ H, d = horizontal_from_equatorial(68.0337, 15.1249, 38.9214)
+ self.assertAlmostEqual(H, 64.352133, 4)
+ self.assertAlmostEqual(d, -6.719892, 4)
<<time and astronomy>>=
@@ -3671,8 +3658,8 @@
def `precise_obliquity(tee):
"""Return precise (mean) obliquity of ecliptic at moment tee."""
u = julian_centuries(tee)/100
- assert(abs(u) < 1,
- 'Error! This formula is valid for +/-10000 years around
J2000.0')
+ #assert(abs(u) < 1,
+ # 'Error! This formula is valid for +/-10000 years around
J2000.0')
return (poly(u, [angle(23, 26, mpf(21.448)),
angle(0, 0, mpf(-4680.93)),
angle(0, 0, mpf(- 1.55)),