[pycalcal] push by enrico.spinielli - Added coordinate conversion and some tests. Added errata pdf for Dec 2... on 2010-01-19 16:26 GMT
pyca...@googlecode.com
Author: Enrico Spinielli <enrico.s...@gmail.com>
Date: Tue Jan 19 08:26:19 2010
Log: Added coordinate conversion and some tests. Added errata pdf for Dec
2009
http://code.google.com/p/pycalcal/source/detail?r=d4e84ff2ca
Added:
/errata20091230.pdf
Modified:
/pycalcal.nw
=======================================
--- /dev/null
+++ /errata20091230.pdf Tue Jan 19 08:26:19 2010
Binary file, no diff available.
=======================================
--- /pycalcal.nw Sat Jan 16 07:53:28 2010
+++ /pycalcal.nw Tue Jan 19 08:26:19 2010
@@ -92,6 +92,7 @@
\def\nw{\textsc{Noweb}}
\def\py{\textsc{Python}}
\def\cl{\textsc{Common LISP}}
+\def\pcc{\textsc{PyCalCal}}
\def\thankstext{I want to thank my wife, Gilda, for the patience
of having looked after all kids and things while I was
@@ -161,11 +162,9 @@
implementation of the calendrical algorithms described in
the book \textit{Calendrical Calculations}~\cite{calendar:calcal}.
According to the authors, the book is textit{the companion} text of
-the algorithms
-implemented in the \cl\ package \texttt{CC3}.
+the algorithms implemented in the \cl\ package \texttt{CC3}.
The source code for \texttt{CC3}, \texttt{calendrica-3.0.cl}, is made
-available electronically
-by the publisher, see \cite{calendar:calcal}.
+available electronically by the publisher, see \cite{calendar:calcal}.
I provide full reference to the original \cl\ and credit (and deep
respect) to its authors, Prof.s~Dershowitz and Reingold.
@@ -3205,89 +3204,166 @@
\rule{\textwidth}{0.005in}
\begin{center}
\includegraphics[width=\textwidth]{fig_alt-az.png}
- \caption{Alt-az coordinate system. $N$ is the North (Celestial) Pole;
- $x$ points at the intersection between Greenwich meridian and the
equator;
- $A$ is the azimuth of the celestial body $B$ as seen from observer
- positioned in $P$ (at geographical longitude $L$ and geographical
latitude
- $\varphi$; $\varphi'$ is the geocentric latitude [mesured from
Earth's
+ \caption{Alt-az coordinate system. $N$ is the North (Celestial)
+ Pole;
+ $x$ points at the intersection between Greenwich meridian and
+ the equator;
+ $A$ is the azimuth of the celestial body $B$ as seen from
+ observer positioned in $P$ (at geographical longitude $L$ and
+ geographical latitude $\varphi$;
+ $\varphi'$ is the geocentric latitude [mesured from Earth's
center])).The plane is tangent in $P$ to the Earth's
- surface.}\label{fig:ALT-AZ}
+ surface.}\label{fig:ALT-AZ}
\end{center}
\rule{\textwidth}{0.005in}
\end{figure}
-Any coordinates given in the alt-az system depend on the place of
observation
-(because the sky appears different from different points on Earth) and on
the
-time of observation (because the Earth rotates, and each star appears to
trace
-out a circle centred on North Celestial Pole).
+Any coordinates given in the alt-az system depend on the place of
+observation (because the sky appears different from different points
+on Earth) and on the time of observation (because the Earth rotates,
+and each star appears to trace out a circle centred on North Celestial
+Pole).
%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{$HA-dec$ Coordinate systems}
\label{sec:HA-dec}
A system of celestial coordinates which is fixed on the sky and
-independent of the observer's time and place can be defined by selecting
-the celestial equator as its fundamental circle.
+independent of the observer's time and place can be defined by
+selecting the celestial equator as its fundamental circle.
The North Celestial Pole (NCP) and the South Celestial Pole (SCP)
lie directly above Earth's North and South Poles.
-The NCP and SCP form the poles of a great circle on the celestial sphere,
-analogous to the equator on Earth: it is called the celestial equator and
it
-lies directly above the Earth's equator.
-
-Any great circle between the NCP and the SCP is a meridian.
-The one which also passes through the zenith and the nadir
-is "the" celestial meridian, or the observer's meridian.
-(It is identical to the principal vertical.)
-This provides our new zero-point;
-in this case, we use the point where it crosses the southern half of the
equator.
+The NCP and SCP form the poles of a great circle on the celestial
+sphere, analogous to the equator on Earth: it is called the celestial
+equator and it lies directly above the Earth's equator.
+
+Any great circle between the NCP and the SCP is a meridian.
+The one which also passes through the zenith and the nadir
+is "the" celestial meridian, or the observer's meridian.
+(It is identical to the principal vertical.)
+This provides our new zero-point;
+in this case, we use the point where it crosses the southern half of
+the equator.
\subsection{$RA-dec$ Coordinate systems}
\label{sec:RA-dec}
-Figure~\ref{fig:RA-DEC}\ displays the celestial sphere and the RA-dec
coordinate
-system (also named second equtorial).
+Figure~\ref{fig:RA-DEC}\ displays the celestial sphere and the RA-dec
+coordinate system (also named second equtorial).
\begin{figure}[h]
\rule{\textwidth}{0.005in}
\begin{center}
\includegraphics{fig_ra-dec-0.pdf}
\caption{RA-dec coordinate system. $N^*$ is the North pole of the
- ecliptic; $N$ is the North Celestial Pole; $\vernal$ is the first
point of
- Aries and the zero-point for this coordinate system. Right Ascension
($RA$
- or $\alpha$) is measured in hours towards East ($E$); declination
($dec$
- or $\delta$) is measured in degrees, positive if north ($N$) of the
- Celestial equator.
- $Sol_w$ is the winter solstice point of $18^h RA$ and $-23.5^{\circ}$
- $dec$; $Sol_s$ is the summer solstice point of $6^h RA$ and
$23.5^{\circ}$
- $dec$.}\label{fig:RA-DEC}
+ ecliptic; $N$ is the North Celestial Pole; $\vernal$ is the
+ first point of Aries and the zero-point for this coordinate
+ system.
+ Right Ascension ($RA$ or $\alpha$) is measured in hours towards
+ East ($E$); declination ($dec$ or $\delta$) is measured in
+ degrees, positive if north ($N$) of the Celestial equator.
+ $Sol_w$ is the winter solstice point of $18^h RA$ and
+ $-23.5^{\circ}$ $dec$;
+ $Sol_s$ is the summer solstice point of $6^h RA$ and
+ $23.5^{\circ} dec$.}\label{fig:RA-DEC}
\end{center}
\rule{\textwidth}{0.005in}
\end{figure}
-Figure~\ref{fig:ecliptic}\ displays ecliptic, lunar orbit (angle is
exagerated),
-the lunar (ascending) node $\ascnode$\ and the origin \vernal\ (at vernal
equinox).
-It is intersting to note graphically what is later described
algorithmically.
+Figure~\ref{fig:ecliptic}\ displays ecliptic, lunar orbit (angle is
+exagerated), the lunar (ascending) node $\ascnode$\ and the origin
+\vernal\ (at vernal equinox).
+It is intersting to note graphically what is later described
+algorithmically.
\begin{figure}[h]
\rule{\textwidth}{0.005in}
\begin{center}
\includegraphics{fig_ecliptic-0.pdf}
\caption{Celestial coordinate system. $N^*$ is the North pole of the
- ecliptic; $N$ is the North Celestial Pole; $\vernal$ is the first
point of
- Aries or the ascending node of the mean ecliptic. $\ascnode$ is the
- ascending node of the lunar orbit.}\label{fig:ecliptic}
+ ecliptic; $N$ is the North Celestial Pole; $\vernal$ is the first
+ point of Aries or the ascending node of the mean
+ ecliptic. $\ascnode$ is the ascending node of the lunar
+ orbit.}\label{fig:ecliptic}
\end{center}
\rule{\textwidth}{0.005in}
-\end{figure}
-@
-
+\end{figure}
<<time and astronomy>>=
######################
# Time and Astronomy #
######################
+def `ecliptical_from_equatorial(decl, ra, obli):
+ """Convert equatorial coordinates to ecliptical ones.
+ Decl is the declination, ra is the right ascension and
+ obl is the obliquity of the ecliptic.
+ NOTE: if 'apparent' right ascension and declination are
+ used, then 'true' obliquity should be input.
+ """
+ lon = degrees_from_radians(
+ atan2(sin_degrees(ra)*cos_degrees(obli) +
+ tan_degrees(decl)*sin_degrees(obli),
+ cos_degrees(ra)))
+ lat = arcsin_degrees(sin_degrees(decl)*cos_degrees(obli) -
+
cos_degrees(decl)*sin_degrees(obli)*sin_degrees(ra))
+ return lon, lat
+
+def `equatorial_from_ecliptical(lam, beta, obli):
+ """Convert ecliptical coordinates to equatorial ones.
+ Lam is the ecliptical longitude, beta is ecliptical latitude and
+ epsi is the obliquity of the ecliptic.
+ NOTE: if 'apparent' right ascension and declination are to be
+ used, then 'true' obliquity should be input.
+ NOTE: beta is positive if north of ecliptic.
+ """
+ ra = degrees_from_radians(
+ atan2(sin_degrees(lam)*cos_degrees(obli) -
+ tan_degrees(beta)*sin_degrees(obli),
+ cos_degrees(lam)))
+ dec = arcsin_degrees(sin_degrees(beta)*cos_degrees(obli) +
+
cos_degrees(beta)*sin_degrees(obli)*sin_degrees(lam))
+ return ra, dec
+
+def `horizontal_from_ecliptical(ra, dec, tee, location, nut=False):
+ """Convert equatorial coordinates to horizontal ones.
+ ra is the right ascension, decl is the declination,
+ location is the geographical position of the observer,
+ tee is a moment in universal time.
+ NOTE: if ra is affected by nutation, then sideral time too
+ must be affected by it.
+ """
+ theta_0 = sidereal_from_moment(tee)
+ if nut:
+ theta_0 += nutation(tee)*cos_degrees(obliquity(tee))/15.0
+ L = longitude(location)
+ H = theta_0 - L - ra
+ 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
+
+
+<<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)
+ 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)
+
+
+
+<<time and astronomy>>=
# see lines 2667-2670 in calendrica-3.0.cl
def `hr(x):
"""Return the number of days given x hours."""
@@ -3347,11 +3423,18 @@
from math import cos
return cos(radians_from_degrees(theta))
+# from errata20091230.pdf entry 112
+cos_degrees=cosine_degrees
+
+
# see lines 2723-2726 in calendrica-3.0.cl
def `tangent_degrees(theta):
"""Return tangent of theta (given in degrees)."""
return tan(radians_from_degrees(theta))
+# from errata20091230.pdf entry 112
+tan_degrees=tangent_degrees
+
def `signum(a):
if a > 0:
@@ -3528,7 +3611,7 @@
# see lines 2872-2880 in calendrica-3.0.cl
def `obliquity(tee):
- """Return obliquity of ecliptic at moment tee."""
+ """Return (mean) obliquity of ecliptic at moment tee."""
c = julian_centuries(tee)
return (angle(23, 26, mpf(21.448)) +
poly(c, [mpf(0),
@@ -3536,6 +3619,29 @@
angle(0, 0, mpf(-0.00059)),
angle(0, 0, mpf(0.001813))]))
+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')
+ return (poly(u, [angle(23, 26, mpf(21.448)),
+ angle(0, 0, mpf(-4680.93)),
+ angle(0, 0, mpf(- 1.55)),
+ angle(0, 0, mpf(+1999.25)),
+ angle(0, 0, mpf(- 51.38)),
+ angle(0, 0, mpf(- 249.67)),
+ angle(0, 0, mpf(- 39.05)),
+ angle(0, 0, mpf(+ 7.12)),
+ angle(0, 0, mpf(+ 27.87)),
+ angle(0, 0, mpf(+ 5.79)),
+ angle(0, 0, mpf(+ 2.45))]))
+
+def `true_obliquity(tee):
+ """Return 'true' obliquity of ecliptic at moment tee.
+ That is, where nutation is taken into accout."""
+ pass
+
+
# see lines 2882-2891 in calendrica-3.0.cl
def `declination(tee, beta, lam):
"""Return declination at moment UT tee of object at
@@ -5853,7 +5959,6 @@
temp = mod(k + 15, 30) - 15
else:
temp = mod(k - 15, 30) + 15
-
est = s + day - temp
tau = (est -
mod(hindu_lunar_day_from_moment(est + hr(6)) - day + 15, 30) +
@@ -7590,6 +7695,35 @@
-chmod ug+x extractcalcalsignatures
@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\chapter{Tutorial}
+\label{sec:tutorial}
+This chapter goes into some practical examples of using \pcc.
+
+%%%%%%%%
+\subsection*{Using dates and times}
+Convert \texttt{30^{th} Jan 1967 11:20:00 UT} to RD
+\begin{verbatim}
+>>> d=gregorian_date(1967, JANUARY, 30)
+>>> t=time_from_clock([11, 20, 00])
+>>> r=fixed_from_gregorian(d) + t
+718096.47222222225
+\end{verbatim}
+
+Print a RD in human readable format
+\begin{verbatim}
+>>> clock_from_moment(r)
+[11, 20, 2.2351741790771484e-06]
+>>> gregorian_from_fixed(fixed_from_moment(r))
+[1967, 1, 30]
+\end{verbatim}
+
+%%%%%%%%
+\subsection*{Converting coordinates}
+
+
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Future evolutions}
@@ -7603,7 +7737,7 @@
service for astronomical applications as made available by the US Naval
Observatory web services \cite{USNO} or the IMCCE~\cite{imcce}.
-I created a google app for pycalcal.
+I created a google app for \pcc.
In a diirectory named 'calendrica' I did put \texttt{calendrica.py}
and \texttt{app.yaml}. I also create a 'lib' directory where I copied
the whole \texttt{mpmath} distribution.