Modified:
/pycalcal.nw
=======================================
--- /pycalcal.nw Tue Feb 2 06:19:37 2010
+++ /pycalcal.nw Sat Feb 6 06:37:45 2010
@@ -3313,11 +3313,19 @@
cos_degrees(declination)*so*sa)
return [lon, lat]
-
+<<astronomical algorithms tests>>=
+def `testEclipticalFromEquatorial(self):
+ # from the values in the Ch 13 Astronomical Algorithms
+ ra, de = equatorial_from_ecliptical(113.215630, 6.684170, 23.4392911)
+ self.assertAlmostEqual(ra, 116.328942, 5)
+ self.assertAlmostEqual(de, 28.026183, 6)
+
+
+<<time and astronomy>>=
def `equatorial_from_ecliptical(longitude, latitude, obliquity):
"""Convert ecliptical coordinates (in degrees) to equatorial ones.
'longitude' is the ecliptical longitude,
- 'latitude' is thestarting ecliptical latitude and
+ 'latitude' is the starting 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.
@@ -3333,7 +3341,15 @@
cos_degrees(latitude)*so*sl)
return [ra, dec]
-
+<<astronomical algorithms tests>>=
+def `testEquatorialFromEcliptical(self):
+ # from the values in the Ch 13 Astronomical Algorithms
+ 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)
+
+
+<<time and astronomy>>=
def `horizontal_from_equatorial(H, declination, latitude):
"""Convert equatorial coordinates (in degrees) to horizontal ones.
Return azimuth and altitude.
@@ -3353,66 +3369,37 @@
cl * cos_degrees(declination) * ch)
return [A, h]
-def `equatorial_from_horizontal(azimuth, altitude, latitude):
+<<astronomical algorithms tests>>=
+def `testHorizontalFromEquatorial(self):
+ # from the values in the Ch 13 Astronomical Algorithms
+ A, h = horizontal_from_equatorial(64.352133, -6.719892, angle(38, 55,
17))
+ self.assertAlmostEqual(A, 68.0337, 4)
+ self.assertAlmostEqual(h, 15.1249, 4)
+
+
+<<time and astronomy>>=
+def `equatorial_from_horizontal(A, h, phi):
"""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.
+ 'A' is the azimuth,
+ 'h' is the altitude,
+ 'phi' is the observer's geographic latitude.
NOTE: azimuth is measured westward from the South.
"""
- # 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)
- a=altitude #altitude
- A=azimuth #azimuth
- phi=latitude #latitude
-
- #cos delta * sin H = cos a * sin A
- # sin delta = sin phi * sin a + cos phi * cos a * cos A
- #cos delta * cos H = cos phi * sin a + sin phi * cos a * cos A
-
- delta = arcsin_degrees(sin_degrees(phi) * sin_degrees(a) +
- cos_degrees(phi) * cos_degrees(a) *
cos_degrees(A))
- H = arcsin_degrees((cos_degrees(a) * sin_degrees(A)) /
cos_degrees(delta))
-
-
- #alpha= #right ascension
- #delta= #declination
- #H= #local hour angle
+ H = degrees_from_radians(
+ atan2(sin_degrees(A),
+ (cos_degrees(A) * sin_degrees(phi) +
+ tan_degrees(h) * cos_degrees(phi))))
+ delta = arcsin_degrees(sin_degrees(phi) * sin_degrees(h) -
+ cos_degrees(phi) * cos_degrees(h) *
cos_degrees(A))
return [H, delta]
-
-
<<astronomical algorithms tests>>=
-def `testEquatorialFromEcliptical(self):
- # from the values in the Ch 13 Astronomical Algorithms
- 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)
-
-def `testEclipticalFromEquatorial(self):
- # from the values in the Ch 13 Astronomical Algorithms
- ra, de = equatorial_from_ecliptical(113.215630, 6.684170, 23.4392911)
- 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)
+ H, d = equatorial_from_horizontal(68.0337, 15.1249, angle(38, 55, 17))
self.assertAlmostEqual(H, 64.352133, 4)
- self.assertAlmostEqual(d, -6.719892, 4)
+ self.assertAlmostEqual(d, degrees(angle(-6,-43,-11.61)), 4)
<<time and astronomy>>=
@@ -3524,7 +3511,7 @@
def `arcsin_degrees(x):
"""Return arcsine of x in degrees."""
from math import asin
- return degrees(degrees_from_radians(asin(x)))
+ return degrees_from_radians(asin(x))
# see lines 2746-2749 in calendrica-3.0.cl
def `arccos_degrees(x):