[pycalcal] push by enrico.spinielli - added mpost macros for figures (from Roegel) in astro.mp. Added sectio... on 2009-12-27 17:33 GMT

[pyca...@googlecode.com]

Revision: 4dfa055450
Author: Enrico Spinielli <enrico.s...@gmail.com>
Date: Sun Dec 27 09:32:53 2009
Log: added mpost macros for figures (from Roegel) in astro.mp. Added
sections for astronomical coordinates.
http://code.google.com/p/pycalcal/source/detail?r=4dfa055450

Added:
/astro.mp
Modified:
/figure.mp
/pycalcal.nw

=======================================
--- /dev/null
+++ /astro.mp Sun Dec 27 09:32:53 2009
@@ -0,0 +1,343 @@
+% astro.mp
+% practically wholly inspired by
+% D. Roegel "Spheres, great circles and parallels", TUGBoat,
+% Vol 30 (2009), Num. 1
+%
+% I organised Roegel's macros in a package waiting for him to
+% publish his own (probably on CTAN?)
+
+% Don't load this package twice:
+if known astro_version: expandafter endinput; fi;
+
+numeric astro_version;string astro_date;
+astro_version=0.1;
+astro_date="2009/12/25";
+% The banner:
+message "** astro " & decimal (astro_version) &
+ " (c) E. Spinielli (wholly inspired by D. Roegel) (" &
+ astro_date & ") **";message "";
+
+tracingstats:=1;
+
+
+
+let vector=color;
+let Xp=redpart; let Yp=greenpart; let Zp=bluepart;
+
+r=5cm;
+theta=70;
+phi=-15;
+
+%==== dot product ====
+def dotproduct(expr Vi,Vj)=
+ (Xp(Vi)*Xp(Vj)+Yp(Vi)*Yp(Vj)+Zp(Vi)*Zp(Vj))
+enddef;
+
+%==== vector product ====
+def vecproduct(expr Vi,Vj)=
+ (Yp(Vi)*Zp(Vj)-Zp(Vi)*Yp(Vj),
+ Zp(Vi)*Xp(Vj)-Xp(Vi)*Zp(Vj),
+ Xp(Vi)*Yp(Vj)-Yp(Vi)*Xp(Vj))
+enddef;
+
+
+vector V[];
+V1=(cosd theta,sind theta,0);
+V2=(sind(phi)*sind(theta),-sind(phi)*cosd(theta),cosd(phi));
+V3=vecproduct(V1,V2);
+
+%==== unit vector ====
+def norm(expr V)= sqrt(dotproduct(V,V)) enddef;
+def normed(expr V)= (V/norm(V)) enddef;
+
+%==== projection ====
+def project(expr V,Va,Vb)=
+ (dotproduct(V,Va),
+ dotproduct(V,Vb))
+enddef;
+
+%==== equator ====
+def f_equator(expr r,t)=
+ (r*cosd(t),r*sind(t),0)
+enddef;
+path equator;
+equator=
+ project(f_equator(r,0),V1,V2)
+ for t=10 step 10 until 350:
+ ..project(f_equator(r,t),V1,V2)
+endfor ..cycle;
+
+%==== ellipse =====
+def ellipse(expr ra,rb,an)=
+ (fullcircle xscaled 2ra yscaled 2rb rotated an)
+enddef;
+
+%==== major angle ellipse ====
+vardef ellipse_major_angle(expr p,a)=
+ save pa,pc,pi,ra,rb,rc,an;path pc[];pair pa,pi[];ra=.5a;rb=a;
+ forever: %================ split ===============
+ rc:=.5[ra,rb];pc0:=subpath(0,4) of fullcircle scaled 2rc;
+ pa:=pc0 intersectiontimes p;exitif pa<>(-1,-1);ra:=rc;
+ endfor;
+ %======= find two intersections ======
+ pi1=p intersectiontimes pc0;
+ pc1=subpath(0,ypart(pi1)-0.01) of pc0;
+ pc2=subpath(ypart(pi1)+0.01,length(pc0)) of pc0;
+ pi1:=p intersectionpoint pc0;pi2:=p intersectiontimes pc1;
+ if pi2=(-1,-1):pi2:=p intersectionpoint pc2;
+ else:pi2:=p intersectionpoint pc1;fi;
+ pi3=pi1 rotated 180;pi4=pi2 rotated 180; % other intersections
+ %======= orientation ======
+ pi5=p intersectionpoint (origin--(unitvector(pi2-pi1)*2a));
+ pi6=p intersectionpoint (origin--(unitvector(pi1-pi4)*2a));
+ if arclength(origin--pi5)>arclength(origin--pi6):an=angle(pi1-pi2);
+ else:an=angle(pi1-pi4);fi;
+ an % result
+enddef;
+
+%==== minor angle ellipse ====
+vardef ellipse_minor_axis(expr p,a,an)=
+ save pa;pair pa;
+ pa=p intersectionpoint (origin--(dir(an+90)*2a));
+ arclength(origin--pa) % result
+enddef;
+
+%==== rotate around ====
+% rotates Va around Vb by the angle a
+vardef rotatearound(expr Va,Vb,a)=
+ save v;vector v[];
+ v0=normed(Vb);v1=dotproduct(Va,v0)*v0;
+ v2=Va-v1;v3=vecproduct(v0,v2);
+ v4=v2*cosd(a)+v3*sind(a)+v1;
+ v4 % result
+enddef;
+
+%==== draw_equator ====
+% both is bool to draw bothsides
+% side is a bool to decide which side to dash
+vardef draw_equator(expr both,side)=
+ save pa,pb,pc,alpha,equatorMinorAxis;
+ path pa,pb,pc;
+ alpha=ellipse_major_angle(equator, r);
+ equatorMinorAxis=ellipse_minor_axis(equator, r, alpha);
+ pa=ellipse(r,equatorMinorAxis,alpha);
+ pb=subpath(0,4) of pa;
+ pc=subpath(4,8) of pa;
+ if both=1:
+ if side=1:
+ draw pb dashed evenly; % hidden
+ draw pc; % visible
+ else:
+ draw pc dashed evenly; % hidden
+ draw pb; % visible
+ fi;
+ else:
+ if side=1:
+ draw pc; % visible
+ else:
+ draw pb; % visible
+ fi;
+ fi;
+enddef;
+
+%==== draw_parallel ====
+% phi=latitude, col=color, side=1 or -1 depending on the dashes
+vardef draw_parallel(expr phi,col,side)=
+ save p;path p[];p0=project(f_parallel(a,0,phi),V1,V2)
+ for t=0 step 10 until 360 :..project(f_parallel(a,t,phi),V1,V2) endfor;
+ % we now search for the intersections of this parallel
+ % with the projection plane:
+ % plane: V3x*x+V3y*y+V3z*z=0
+ % parallel: x=r*cos(phi)*cos(theta), y=r*cos(phi)*sin(theta),
z=r*sin(phi)
+ % we search theta:
+ save A,B,C,X,Y,ca,cb,cc,delta,nx,tha,thb;
+ numeric X[],Y[];ca=Xp(V3);cb=Yp(V3);cc=Zp(V3);
+ if cb=0:X1=-(cc/ca)*sind(phi)/cos(phi);nx=1;
+ else:
+ A=1+(ca/cb)**2;B=2*ca*cc*sind(phi)/(cb*cb);
+ C=((cc/cb)*sind(phi))**2-cosd(phi)*cosd(phi);delta=B*B-4A*C;
+ if delta<0:nx=0;% no intersection
+ else:
+ X1=((-B-sqrt(delta))/(2A))/cosd(phi); % = cos(theta)
+ X2=((-B+sqrt(delta))/(2A))/cosd(phi); % = cos(theta)
+ Y1=-((ca*X1+cc*sind(phi)/cosd(phi))/cb); % = sin(theta)
+ Y2=-((ca*X2+cc*sind(phi)/cosd(phi))/cb); % = sin(theta)
+ tha=angle(X1,Y1);thb=angle(X2,Y2);nx=2;
+ fi;
+ fi;
+ if nx=0: % totally (in)visible parallel
+ if side=1:draw p0 withcolor col;
+ else:draw p0 withcolor col dashed evenly;fi;
+ message "NO INTERSECTION";
+ elseif nx=1:X10=angle(X1,1+-+X1);X11=360-X10;
+ else: % general case
+ if tha<thb:X10=tha;X11=thb;else:X10=thb;X11=tha;fi;
+ fi;
+ if nx>0: % determination of the two paths
+ p1=project(f_parallel(a,X10,phi),V1,V2)
+ for t=X10+1 step 10 until X11:..project(f_parallel(a,t,phi),V1,V2)
+ endfor;
+ p2=project(f_parallel(a,X11,phi),V1,V2)
+ for t=X11+1 step 10 until
X10+360:..project(f_parallel(a,t,phi),V1,V2)
+ endfor;
+ % drawing the two paths
+ if side=1:draw p1 withcolor col;
+ else:draw p1 withcolor col dashed evenly;fi;
+ if side=1:draw p2 withcolor col dashed evenly;
+ else:draw p2 withcolor col;fi;
+ fi;
+enddef;
+
+%==== ecliptic ====
+ec_angle=23.5;
+def f_ecliptic(expr r,t)=
+ (r*(cosd(t),sind(t)*cosd(ec_angle),
+ sind(t)*sind(ec_angle)))
+enddef;
+path ecliptic;
+ecliptic=
+ project(f_ecliptic(r,0),V1,V2)
+ for t=10 step 10 until 350:
+ ..project(f_ecliptic(r,t),V1,V2)
+endfor ..cycle;
+
+%==== draw_ecliptic ====
+% both 1 to draw both sides
+% side is a bool to decide which side to dash
+vardef draw_ecliptic(expr both,side)=
+ save pd,pe,pf,beta,eclipticMinorAxis;
+ path pd,pe,pf;
+ beta=ellipse_major_angle(ecliptic, r);
+ eclipticMinorAxis=ellipse_minor_axis(ecliptic, r, beta);
+ pd=ellipse(r,eclipticMinorAxis,beta);
+ pe=subpath(0,4) of pd;
+ pf=subpath(4,8) of pd;
+ if both=1:
+ if side=1:
+ draw pe dashed evenly; % hidden
+ draw pf; % visible
+ else:
+ draw pf dashed evenly; % hidden
+ draw pe; % visible
+ fi;
+ else:
+ if side=1:
+ draw pf; % visible
+ else:
+ draw pe; % visible
+ fi;
+ fi;
+enddef;
+
+%==== lunar orbit ====
+% moon_angle=5.145; % true value
+moon_angle=15; % exagerated value
+Ln=18; % lunar node (angle)
+vector B; B=f_ecliptic(r,Ln);
+def f_lunar(expr r,t)=
+ rotatearound(r*(cosd(t+Ln),
+ sind(t+Ln)*cosd(ec_angle),
+ sind(t+Ln)*sind(ec_angle)),B,moon_angle)
+enddef;
+path lunar;
+lunar=
+ project(f_lunar(r,0),V1,V2)
+ for t=10 step 10 until 350:
+ ..project(f_lunar(r,t),V1,V2)
+endfor ..cycle;
+
+%==== draw_lunar_orbit ====
+% both 1 to draw both sides
+% side is a bool to decide which side to dash
+vardef draw_lunar_orbit(expr both, side)=
+ save pn,po,pp,rho,lunarMinorAxis;
+ path pn,po,pp;
+ rho=ellipse_major_angle(lunar, r);
+ lunarMinorAxis=ellipse_minor_axis(lunar, r, rho);
+ pn=ellipse(r,lunarMinorAxis,rho);
+ po=subpath(0,4) of pn;
+ pp=subpath(4,8) of pn;
+ if both=1:
+ if side=1:
+ draw po dashed evenly; % hidden
+ draw pp; % visible
+ else:
+ draw pp dashed evenly; % hidden
+ draw po; % visible
+ fi;
+ else:
+ if side=1:
+ draw pp; % visible
+ else:
+ draw po; % visible
+ fi;
+ fi;
+enddef;
+
+%==== North, North Ecliptic ====
+vector North,North_Ec, North_Moon;
+North=r*(0,0,1);
+North_Ec=rotatearound(North,(1,0,0),ec_angle);
+North_Moon=rotatearound(North,(1,0,0),moon_angle);
+
+%==== ecliptic meridian ====
+% A is a point in space on the ecliptic
+% t is an angle
+def f_ec_meridian(expr t,A)=
+ (A*cosd(t)+North_Ec*sind(t))
+enddef;
+path ec_meridian;
+vector A; A=f_ecliptic(r,48); % point (in space) on ecliptic
+ec_meridian=
+ project(f_ec_meridian(0,A),V1,V2)
+ for t=10 step 10 until 350:
+ ..project(f_ec_meridian(t,A),V1,V2)
+endfor ..cycle;
+
+%==== meridian ====
+% r radius
+% lon is an angle of longitude
+% t is an angle (latitude)
+def f_meridian(expr r,t,lon)=
+ (r*(cosd(lon)*sind(t),sind(lon)*sind(t),cosd(t)))
+enddef;
+
+%==== draw_meridian ====
+% phi is the longitude of the meridian
+% both is bool to draw bothsides
+% side is a bool to decide which side to dash
+vardef draw_meridian(expr phi,both,side)=
+ save pk,pl,pm,meridian,lon,v,lambda,meridianMinorAxis;
+ vector v;
+ path pk,pl,pm,meridian;
+ v=(cosd(phi),sind(phi),0);
+ lon=angle(Xp(v),Yp(v));
+ meridian=
+ project(f_meridian(r,0,lon),V1,V2)
+ for t=10 step 10 until 350:
+ ..project(f_meridian(r,t,lon),V1,V2)
+ endfor ..cycle;
+ lambda=ellipse_major_angle(meridian, r);
+ meridianMinorAxis=ellipse_minor_axis(meridian, r, lambda);
+ pk=ellipse(r,meridianMinorAxis,lambda);
+ pl=subpath(0,4) of pk;
+ pm=subpath(4,8) of pk;
+ if both=1:
+ if side=1:
+ draw pl dashed evenly; % hidden
+ draw pm; % visible
+ else:
+ draw pm dashed evenly; % hidden
+ draw pl; % visible
+ fi;
+ else:
+ if side=1:
+ draw pm;
+ else:
+ draw pl;
+ fi;
+ fi;
+enddef;
+
+
+endinput
=======================================
--- /figure.mp Sun Dec 13 04:52:02 2009
+++ /figure.mp Sun Dec 27 09:32:53 2009
@@ -1,318 +1,161 @@
verbatimtex
%&latex
\documentclass[12pt]{article}
-\usepackage{wasysym}
+ \usepackage{wasysym}
\begin{document}
etex

+input astro
+input textpath
beginfig(0);
- r=5cm;
- theta=70;phi=-15;
-
- let vector=color;
- let Xp=redpart; let Yp=greenpart; let Zp=bluepart;
-
-
- %==== dot product ====
- def dotproduct(expr Vi,Vj)=
- (Xp(Vi)*Xp(Vj)+Yp(Vi)*Yp(Vj)+Zp(Vi)*Zp(Vj))
- enddef;
-
- %==== vector product ====
- def vecproduct(expr Vi,Vj)=
- (Yp(Vi)*Zp(Vj)-Zp(Vi)*Yp(Vj),
- Zp(Vi)*Xp(Vj)-Xp(Vi)*Zp(Vj),
- Xp(Vi)*Yp(Vj)-Yp(Vi)*Xp(Vj))
- enddef;
-
-
- vector V[]; % tableau de vecteurs
- V1=(cosd theta,sind theta,0);
- V2=(sind(phi)*sind(theta),-sind(phi)*cosd(theta),cosd(phi));
- V3=vecproduct(V1,V2);
-
- %==== unit vector ====
- def norm(expr V)= sqrt(dotproduct(V,V)) enddef;
- def normed(expr V)= (V/norm(V)) enddef;
-
- %==== projection ====
- def project(expr V,Va,Vb)=
- (dotproduct(V,Va),
- dotproduct(V,Vb))
- enddef;
-
-
-
- %==== equator ====
- def f_equ(expr r,t)=
- (r*cosd(t),r*sind(t),0)
- enddef;
- path equator;
- equator=
- project(f_equ(r,0),V1,V2)
- for t=10 step 10 until 350:
- ..project(f_equ(r,t),V1,V2)
- endfor ..cycle;
-
-
- %==== ellipse =====
- def ellipse(expr ra,rb,an)=
- (fullcircle xscaled 2ra yscaled 2rb rotated an)
- enddef;
-
-
- %==== major angle ellipse ====
- vardef ellipse_major_angle(expr p,a)=
- save pa,pc,pi,ra,rb,rc,an;path pc[];pair pa,pi[];ra=.5a;rb=a;
- forever: %================ dichotomie ===============
- rc:=.5[ra,rb];pc0:=subpath(0,4) of fullcircle scaled 2rc;
- pa:=pc0 intersectiontimes p;exitif pa<>(-1,-1);ra:=rc;
- endfor;
- %======= calcul de deux intersections ======
- pi1=p intersectiontimes pc0;
- pc1=subpath(0,ypart(pi1)-0.01) of pc0;
- pc2=subpath(ypart(pi1)+0.01,length(pc0)) of pc0;
- pi1:=p intersectionpoint pc0;pi2:=p intersectiontimes pc1;
- if pi2=(-1,-1):pi2:=p intersectionpoint pc2;
- else:pi2:=p intersectionpoint pc1;fi;
- pi3=pi1 rotated 180;pi4=pi2 rotated 180; % autres intersections
- %======= orientation ======
- pi5=p intersectionpoint (origin--(unitvector(pi2-pi1)*2a));
- pi6=p intersectionpoint (origin--(unitvector(pi1-pi4)*2a));
- if arclength(origin--pi5)>arclength(origin--pi6):an=angle(pi1-pi2);
- else:an=angle(pi1-pi4);fi;
- an % r?ɬ©sultat de la macro
- enddef;
-
- %==== minor angle ellipse ====
- vardef ellipse_minor_axis(expr p,a,an)=
- save pa;pair pa;
- pa=p intersectionpoint (origin--(dir(an+90)*2a));
- arclength(origin--pa) % r?ɬ©sultat
- enddef;
-
- %==== rotate around ====
- % rotates Va around Vb by the angle ?Ǭ°a?Ǭ¢
- vardef rotatearound(expr Va,Vb,a)=
- save v;vector v[];
- v0=normed(Vb);v1=dotproduct(Va,v0)*v0;
- v2=Va-v1;v3=vecproduct(v0,v2);
- v4=v2*cosd(a)+v3*sind(a)+v1;
- v4 % result
- enddef;
-
- %==== draw_parallel ====
- % phi=latitude, col=color, side=1 or -1 depending on the dashes
- vardef draw_parallel(expr phi,col,side)=
- save p;path p[];p0=project(f_parallel(a,0,phi),V1,V2)
- for t=0 step 10 until 360 :..project(f_parallel(a,t,phi),V1,V2)
endfor;
- % we now search for the intersections of this parallel
- % with the projection plane:
- % plane: V3x*x+V3y*y+V3z*z=0
- % parallel: x=r*cos(phi)*cos(theta), y=r*cos(phi)*sin(theta),
z=r*sin(phi)
- % we search theta:
- save A,B,C,X,Y,ca,cb,cc,delta,nx,tha,thb;
- numeric X[],Y[];ca=Xp(V3);cb=Yp(V3);cc=Zp(V3);
- if cb=0:X1=-(cc/ca)*sind(phi)/cos(phi);nx=1;
- else:
- A=1+(ca/cb)**2;B=2*ca*cc*sind(phi)/(cb*cb);
- C=((cc/cb)*sind(phi))**2-cosd(phi)*cosd(phi);delta=B*B-4A*C;
- if delta<0:nx=0;% no intersection
- else:
- X1=((-B-sqrt(delta))/(2A))/cosd(phi); % = cos(theta)
- X2=((-B+sqrt(delta))/(2A))/cosd(phi); % = cos(theta)
- Y1=-((ca*X1+cc*sind(phi)/cosd(phi))/cb); % = sin(theta)
- Y2=-((ca*X2+cc*sind(phi)/cosd(phi))/cb); % = sin(theta)
- tha=angle(X1,Y1);thb=angle(X2,Y2);nx=2;
- fi;
- fi;
- if nx=0: % totally (in)visible parallel
- if side=1:draw p0 withcolor col;
- else:draw p0 withcolor col dashed evenly;fi;
- message "NO INTERSECTION";
- elseif nx=1:X10=angle(X1,1+-+X1);X11=360-X10;
- else: % general case
- if tha<thb:X10=tha;X11=thb;else:X10=thb;X11=tha;fi;
- fi;
- if nx>0: % determination of the two paths
- p1=project(f_parallel(a,X10,phi),V1,V2)
- for t=X10+1 step 10 until X11:..project(f_parallel(a,t,phi),V1,V2)
- endfor;
- p2=project(f_parallel(a,X11,phi),V1,V2)
- for t=X11+1 step 10 until
X10+360:..project(f_parallel(a,t,phi),V1,V2)
- endfor;
- % drawing the two paths
- if side=1:draw p1 withcolor col;
- else:draw p1 withcolor col dashed evenly;fi;
- if side=1:draw p2 withcolor col dashed evenly;
- else:draw p2 withcolor col;fi;
- fi;
- enddef;
-
- %==== ecliptic ====
- ec_angle=23.5;
- def f_ecliptic(expr r,t)=
- (r*(cosd(t),sind(t)*cosd(ec_angle),
- sind(t)*sind(ec_angle)))
- enddef;
- path ecliptic;
- ecliptic=
- project(f_ecliptic(r,0),V1,V2)
- for t=10 step 10 until 350:
- ..project(f_ecliptic(r,t),V1,V2)
- endfor ..cycle;
-
- %==== lunar orbit ====
-% moon_angle=5.145; % true value
- moon_angle=15; % exagerated value
- Ln=18; % lunar node (angle)
- vector B; B=f_ecliptic(r,Ln);
- def f_lunar(expr r,t)=
- rotatearound(r*(cosd(t+Ln),
- sind(t+Ln)*cosd(ec_angle),
- sind(t+Ln)*sind(ec_angle)),B,moon_angle)
- enddef;
- path lunar;
- lunar=
- project(f_lunar(r,0),V1,V2)
- for t=10 step 10 until 350:
- ..project(f_lunar(r,t),V1,V2)
- endfor ..cycle;
-
- %==== North, North Ecliptic ====
- vector North,North_Ec, North_Moon;
- North=r*(0,0,1);
- North_Ec=rotatearound(North,(1,0,0),ec_angle);
- North_Moon=rotatearound(North,(1,0,0),moon_angle);
-
- %==== ecliptic meridian ====
- % A is a point in space on the ecliptic
- % t is an angle
- %
- def f_ec_meridian(expr t,A)=
- (A*cosd(t)+North_Ec*sind(t))
- enddef;
- path ec_meridian;
- vector A; A=f_ecliptic(r,48); % point (in space) on ecliptic
- ec_meridian=
- project(f_ec_meridian(0,A),V1,V2)
- for t=10 step 10 until 350:
- ..project(f_ec_meridian(t,A),V1,V2)
- endfor ..cycle;
-
- %==== meridian ====
- % lon is an angle of longitude
- % t is an angle (latitude)
- def f_meridian(expr t,lon)=
- (r*(cosd(lon)*sind(t),sind(lon)*sind(t),cosd(t)))
- enddef;
- path meridian;
- Alon=angle(Xp(A),Yp(A));
- meridian=
- project(f_meridian(0,Alon),V1,V2)
- for t=10 step 10 until 350:
- ..project(f_meridian(t,Alon),V1,V2)
- endfor ..cycle;
-
-
-
-
-%%%=============== drawings ============
-
-
+draw fullcircle scaled 2r;
+draw_equator(1,1);
+draw_ecliptic(1,1);
+
+
+z0=(0,0); % origin
+z1=project((r,0,0),V1,V2); % x-axis
+z2=project((0,r,0),V1,V2); % y-axis
+z3=project((0,0,r),V1,V2); % z-axis
+z4=project(North_Ec,V1,V2);
+z6=project(f_ecliptic(r,75),V1,V2);
+z7=project(f_equator(r,70),V1,V2);
+
+path p[],q[];
+vector ra[], dec[];
+ra0=f_equator(r,0);
+ra1=f_equator(r,15);
+ra2=f_equator(r,30);
+ra3=f_equator(r,45);
+ra4=f_equator(r,60);
+ra5=f_equator(r,285);
+ra6=f_equator(r,300);
+ra7=f_equator(r,315);
+ra8=f_equator(r,330);
+ra9=f_equator(r,345);
+p0:=project(ra0,V1,V2)..project(ra1,V1,V2)..project(ra2,V1,V2)..project(0.5[ra2,ra3],V1,V2);
+draw_meridian(0,1,1);
+draw project(ra1-(0,0,2),V1,V2)..project(ra1+(0,0,2),V1,V2);
+label.bot(btex $1^h$ etex, project(ra1-(0,0,2),V1,V2));
+draw project(ra2-(0,0,2),V1,V2)..project(ra2+(0,0,2),V1,V2);
+label.bot(btex $2^h$ etex, project(ra2-(0,0,2),V1,V2));
+draw project(ra3-(0,0,2),V1,V2)..project(ra3+(0,0,2),V1,V2);
+label.bot(btex $3^h$ etex, project(ra3-(0,0,2),V1,V2));
+draw project(ra3-(0,0,2),V1,V2)..project(ra3+(0,0,2),V1,V2);
+
+label.bot(btex $19^h$ etex, project(ra5-(0,0,2),V1,V2));
+draw project(ra5-(0,0,2),V1,V2)..project(ra5+(0,0,2),V1,V2);
+label.bot(btex $20^h$ etex, project(ra6-(0,0,2),V1,V2));
+draw project(ra6-(0,0,2),V1,V2)..project(ra6+(0,0,2),V1,V2);
+label.bot(btex $21^h$ etex, project(ra7-(0,0,2),V1,V2));
+draw project(ra7-(0,0,2),V1,V2)..project(ra7+(0,0,2),V1,V2);
+label.bot(btex $22^h$ etex, project(ra8-(0,0,1),V1,V2));
+draw project(ra8-(0,0,2),V1,V2)..project(ra8+(0,0,2),V1,V2);
+label.bot(btex $23^h$ etex, project(ra9-(0,0,1),V1,V2));
+draw project(ra9-(0,0,2),V1,V2)..project(ra9+(0,0,2),V1,V2);
+
+picture pa;
+%pa = thelabel.bot(btex $24^h\!\!=\!0^h$ etex, project(ra0-(0,0,2),V1,V2));
+pa = thelabel.bot(btex $0^h$ etex, project(ra0-(0,-2.3,2.3),V1,V2));
+unfill bbox pa;
+draw pa;
+
+dec0=f_meridian(r,90,0);
+dec1=f_meridian(r,80,0);
+dec2=f_meridian(r,70,0);
+dec3=f_meridian(r,60,0);
+dec4=f_meridian(r,50,0);
+dec5=f_meridian(r,40,0);
+dec6=f_meridian(r,30,0);
+dec7=f_meridian(r,100,0);
+dec8=f_meridian(r,110,0);
+dec9=f_meridian(r,120,0);
+q0:=project(dec0,V1,V2)..project(dec1,V1,V2)..project(dec2,V1,V2)..project(dec3,V1,V2)..project(0.5[dec3,
dec4],V1,V2);
+draw project(dec1-(0,2,0),V1,V2)..project(dec1+(0,2,0),V1,V2);
+label.rt(btex $10^{\circ}$ etex, project(dec1+(0,0,2),V1,V2));
+draw project(dec2-(0,2,0),V1,V2)..project(dec2+(0,2,0),V1,V2);
+label.rt(btex $20^{\circ}$ etex, project(dec2+(0,0,2),V1,V2));
+draw project(dec3-(0,2,0),V1,V2)..project(dec3+(0,2,0),V1,V2);
+label.rt(btex $30^{\circ}$ etex, project(dec3+(0,0,2),V1,V2));
+draw project(dec4-(0,2,0),V1,V2)..project(dec4+(0,2,0),V1,V2);
+label.rt(btex $40^{\circ}$ etex, project(dec4+(0,0,2),V1,V2));
+draw project(dec5-(0,2,0),V1,V2)..project(dec5+(0,2,0),V1,V2);
+label.rt(btex $50^{\circ}$ etex, project(dec5+(0,0,2),V1,V2));
+draw project(dec6-(0,2,0),V1,V2)..project(dec6+(0,2,0),V1,V2);
+label.rt(btex $60^{\circ}$ etex, project(dec6+(0,0,2),V1,V2));
+draw project(dec7-(0,2,0),V1,V2)..project(dec7+(0,2,0),V1,V2);
+label.lft(btex $-10^{\circ}$ etex, project(dec7+(0,1,0),V1,V2));
+draw project(dec8-(0,2,0),V1,V2)..project(dec8+(0,2,0),V1,V2);
+label.lft(btex $-20^{\circ}$ etex, project(dec8+(0,1,0),V1,V2));
+draw project(dec9-(0,2,0),V1,V2)..project(dec9+(0,2,0),V1,V2);
+label.lft(btex $-30^{\circ}$ etex, project(dec9+(0,1,0),V1,V2));
+
+dec0:=f_meridian(r,130,0);
+draw project(dec0-(0,2,0),V1,V2)..project(dec0+(0,2,0),V1,V2);
+label.lft(btex $-40^{\circ}$ etex, project(dec0+(0,1,0),V1,V2));
+
+drawarrow z0--0.2z1; label.bot(btex $x$ etex, 0.2z1);
+drawarrow z0--0.2z2; label.rt(btex $y$ etex, 0.2z2);
+drawarrow z0--0.2z3; label.lft(btex $z$ etex, 0.2z3);
+pickup pencircle scaled 1.2pt;
+drawarrow p0;
+label.top(btex $RA$ etex, (point infinity of p0)+(-1,1));
+drawarrow q0;
+label.lft(btex $dec$ etex, point infinity of q0);
+
+dotlabel.ulft(btex $\vernal$ etex, z1);
+dotlabel.lrt(btex $\libra$ etex, -z1);
+dotlabel.llft(btex $N$ etex, z3);
+dotlabel.urt(btex $S$ etex, -z3);
+dotlabel.llft(btex $E$ etex, z2);
+dotlabel.urt(btex $W$ etex, -z2);
+dotlabel.llft(btex $N^*$ etex, z4);
+dotlabel.urt(btex $S^*$ etex, -z4);
+dotlabel.urt(btex $Sol_w$ etex, project(f_ecliptic(r,270),V1,V2));
+dotlabel.llft(btex $Sol_s$ etex, project(f_ecliptic(r,90),V1,V2));
+
+pickup defaultpen;
+label.rt(btex Ecliptic etex, 1.01z6);
+label.rt(btex Celestial equator etex, 1.01z7);
+endfig;
+
+beginfig(1);
draw fullcircle scaled 2r;
-
- z0=(0,0);
- z1=project((r,0,0),V1,V2);
- z2=project((0,r,0),V1,V2);
- z3=project((0,0,r),V1,V2);
+ draw_equator(1,1);
+ draw_ecliptic(1,1);
+ draw_lunar_orbit(1,1);
+ draw_meridian(0,1,1);
+
+ z0=(0,0); % origin
+ z1=project((r,0,0),V1,V2); % x-axis
+ z2=project((0,r,0),V1,V2); % y-axis
+ z3=project((0,0,r),V1,V2); % z-axis
+ z4=project(North_Ec,V1,V2);
+ z5=project(f_ecliptic(r,Ln),V1,V2);
+ z6=project(f_ecliptic(r,75),V1,V2);
+ z7=project(f_equator(r,70),V1,V2);
+ z8=project(f_ecliptic(r,310),V1,V2);
+ z9=project(f_ecliptic(r,315),V1,V2);
+
+ drawarrow z0--0.2z1; label.bot(btex $x$ etex, 0.2z1);
+ drawarrow z0--0.2z2; label.rt(btex $y$ etex, 0.2z2);
+ drawarrow z0--0.2z3; label.lft(btex $z$ etex, 0.2z3);
draw -z1--z1 dashed evenly;
draw -z3--z3 dashed evenly;
- drawdot z0;
- dotlabel.top(btex $\vernal$ etex, z1);
- label.rt(btex equator etex, 1.05z2);
+
pickup pencircle scaled 2pt;
- drawdot z1;
- drawdot z3;
- drawdot -z3;
+ dotlabel.urt(btex $\vernal$ etex, z1);
+ dotlabel.bot(btex $\libra$ etex, -z1);
dotlabel.llft(btex $N$ etex, z3);
dotlabel.urt(btex $S$ etex, -z3);
- z4=project(North_Ec,V1,V2);
- drawdot z4;
- drawdot -z4;
dotlabel.llft(btex $N^*$ etex, z4);
dotlabel.urt(btex $S^*$ etex, -z4);
- z5=project(B,V1,V2);
- drawdot z5;
dotlabel.lrt(btex $\ascnode$ etex, z5);
- pickup defaultpen;
-
-
-
- %-- economic equator --
- path pa,pb,pc;
- alpha=ellipse_major_angle(equator, r);
- equatorMinorAxis=ellipse_minor_axis(equator, r, alpha);
- pa=ellipse(r,equatorMinorAxis,alpha);
- pb=subpath(0,4) of pa;
- pc=subpath(4,8) of pa;
- draw pb dashed evenly; % hidden
- draw pc; % visible
-
- %-- economic ecliptic --
- path pd,pe,pf;
- beta=ellipse_major_angle(ecliptic, r);
- eclipticMinorAxis=ellipse_minor_axis(ecliptic, r, beta);
- pd=ellipse(r,eclipticMinorAxis,beta);
- pe=subpath(0,4) of pd;
- pf=subpath(4,8) of pd;
- draw pe dashed evenly; % hidden
- draw pf; % visible
-
- %-- economic ecliptic meridian thru A --
- path pg,ph,pj;
- gamma=ellipse_major_angle(ec_meridian, r);
- eclipticMeridianMinorAxis=ellipse_minor_axis(ec_meridian, r, gamma);
- pg=ellipse(r,eclipticMeridianMinorAxis,gamma);
- ph=subpath(0,4) of pg;
- pj=subpath(4,8) of pg;
-% draw pj dashed evenly; % hidden
-% draw ph; % visible
-
- %-- economic meridian thru A --
- path pk,pl,pm;
- lambda=ellipse_major_angle(meridian, r);
- meridianMinorAxis=ellipse_minor_axis(meridian, r, lambda);
- pk=ellipse(r,meridianMinorAxis,lambda);
- pl=subpath(0,4) of pk;
- pm=subpath(4,8) of pk;
-% draw pl dashed evenly; % hidden
-% draw pm; % visible
-
- %-- economic lunar orbit --
- path pn,po,pp;
- rho=ellipse_major_angle(lunar, r);
- lunarMinorAxis=ellipse_minor_axis(lunar, r, rho);
- pn=ellipse(r,lunarMinorAxis,rho);
- po=subpath(0,4) of pn;
- pp=subpath(4,8) of pn;
- draw po dashed evenly; % hidden
- draw pp; % visible
-
-
-
- pickup pencircle scaled 2pt;
- pair Lone;
-% Lone=project(A,V1,V2);
-% drawdot Lone;
-% dotlabel.lrt(btex $L_1$ etex, Lone);
+
pickup defaultpen;
-
-
-
- % ==========
-% Lp2=moon intersectionpoint ec_meridian;
-% L2=m_x*moon_x+m_y*moon_y;
-% Lp2=m_x*project(moon_x,V1,V2)+m_y*project(moon_y,V1,V2);
+ label.rt(btex Celestial equator etex, 1.01z7);
+ label.rt(btex Ecliptic etex, 1.01z6);
+ label.rt(btex Lunar orbit etex, 1.02*project(f_lunar(r,70),V1,V2));
+ drawarrow z8..z9;
endfig;
end
=======================================
--- /pycalcal.nw Sun Dec 13 04:52:02 2009
+++ /pycalcal.nw Sun Dec 27 09:32:53 2009
@@ -1,7 +1,8 @@
%-*- mode: Noweb; noweb-code-mode: python-mode; tab-width: 4 -*-
% draft is convinient to show over/under boxes,
% BUT it does NOT inlude graphics nor create links
-\documentclass[a4paper,draft,pdftex]{report}
+%\documentclass[a4paper,draft,pdftex]{report}
+\documentclass[a4paper,pdftex]{report}

\usepackage{wasysym}
%%\usepackage[a4paper,top=3cm,bottom=3cm,left=3.5cm,right=3.5cm,%
@@ -3140,6 +3141,34 @@
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Time and Astronomy}
+\label{sec:astro}
+\subsection{$alt-az$ Coordinate systems}
+\label{sec:alt-az}
+
+\subsection{$HA-dec$ Coordinate systems}
+\label{sec:HA-dec}
+
+\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).
+\begin{figure}[h]
+ \rule{\textwidth}{0.005in}
+ \begin{center}
+ \includegraphics{figure-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}
+ \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.
@@ -3147,8 +3176,11 @@
\begin{figure}[h]
\rule{\textwidth}{0.005in}
\begin{center}
- \includegraphics{figure-0.pdf}
- \caption{}\label{fig:ecliptic}
+ \includegraphics{figure-1.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}
\end{center}
\rule{\textwidth}{0.005in}
\end{figure}
@@ -3195,7 +3227,8 @@
# see lines 2692-2696 in calendrica-3.0.cl
def `angle(d, m, s):
"""Return an angle data structure
- from d degrees, m arcminutes and s arcseconds."""
+ from d degrees, m arcminutes and s arcseconds.
+ This assumes that negative angles specifies negative d, m and s."""
return d + ((m + (s / 60)) / 60)

# see lines 2698-2701 in calendrica-3.0.cl
@@ -3306,6 +3339,8 @@
# see lines 5898-5901 in calendrica-3.0.cl
`JERUSALEM = location(31.8, 35.2, mt(800), hr(2))

+`BRUXELLES = location(angle(4, 21, 17), angle(50, 50, 47), mt(800), hr(1))
+
`URBANA = location(40.1,
-88.2,
mt(225),
@@ -4063,7 +4098,8 @@
# see lines 3367-3376 in calendrica-3.0.cl
def `mean_lunar_longitude(c):
"""Return mean longitude of moon (in degrees) at moment
- given in Julian centuries c.
+ given in Julian centuries c (including the constant term of the
+ effect of the light-time (-0".70).
Adapted from eq. 47.1 in "Astronomical Algorithms" by Jean Meeus,
Willmann_Bell, Inc., 2nd ed. with corrections, 2005."""
return degrees(poly(c,deg([mpf(218.3164477), mpf(481267.88123421),
@@ -4200,11 +4236,118 @@
return mod(cap_L_prime + correction + venus +
jupiter + flat_earth + nutation(tee), 360)

+# see lines 3663-3732 in calendrica-3.0.cl
+def `lunar_latitude(tee):
+ """Return the latitude of moon (in degrees) at moment, tee.
+ Adapted from "Astronomical Algorithms" by Jean Meeus,
+ Willmann_Bell, Inc., 1998."""
+ c = julian_centuries(tee)
+ cap_L_prime = mean_lunar_longitude(c)
+ cap_D = lunar_elongation(c)
+ cap_M = solar_anomaly(c)
+ cap_M_prime = lunar_anomaly(c)
+ cap_F = moon_node(c)
+ cap_E = poly(c, [1, mpf(-0.002516), mpf(-0.0000074)])
+ args_lunar_elongation = \
+ [0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 0, 4, 0,
0, 0,
+ 1, 0, 0, 0, 1, 0, 4, 4, 0, 4, 2, 2, 2, 2, 0, 2, 2, 2, 2, 4,
2, 2,
+ 0, 2, 1, 1, 0, 2, 1, 2, 0, 4, 4, 1, 4, 1, 4, 2]
+ args_solar_anomaly = \
+ [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, -1, -1, 1, 0,
1,
+ 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1,
1,
+ 0, -1, -2, 0, 1, 1, 1, 1, 1, 0, -1, 1, 0, -1, 0, 0, 0, -1, -2]
+ args_lunar_anomaly = \
+ [0, 1, 1, 0, -1, -1, 0, 2, 1, 2, 0, -2, 1, 0, -1, 0, -1, -1,
-1,
+ 0, 0, -1, 0, 1, 1, 0, 0, 3, 0, -1, 1, -2, 0, 2, 1, -2, 3, 2,
-3,
+ -1, 0, 0, 1, 0, 1, 1, 0, 0, -2, -1, 1, -2, 2, -2, -1, 1, 1,
-2,
+ 0, 0]
+ args_moon_node = \
+ [1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1,
-1,
+ -1, 1, 3, 1, 1, 1, -1, -1, -1, 1, -1, 1, -3, 1, -3, -1, -1, 1,
+ -1, 1, -1, 1, 1, 1, 1, -1, 3, -1, -1, 1, -1, -1, 1, -1, 1, -1,
+ -1, -1, -1, -1, -1, 1]
+ sine_coefficients = \
+ [5128122, 280602, 277693, 173237, 55413, 46271, 32573,
+ 17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211,
+ 2065, -1870, 1828, -1794, -1749, -1565, -1491, -1475,
+ -1410, -1344, -1335, 1107, 1021, 833, 777, 671, 607,
+ 596, 491, -451, 439, 422, 421, -366, -351, 331, 315,
+ 302, -283, -229, 223, 223, -220, -220, -185, 181,
+ -177, 176, 166, -164, 132, -119, 115, 107]
+ beta = (deg(1/1000000) *
+ sigma([sine_coefficients,
+ args_lunar_elongation,
+ args_solar_anomaly,
+ args_lunar_anomaly,
+ args_moon_node],
+ lambda v, w, x, y, z: (v *
+ pow(cap_E, abs(x)) *
+ sin_degrees((w * cap_D) +
+ (x * cap_M) +
+ (y * cap_M_prime) +
+ (z * cap_F)))))
+ venus = (deg(175/1000000) *
+ (sin_degrees(deg(mpf(119.75)) + c * deg(mpf(131.849)) +
cap_F) +
+ sin_degrees(deg(mpf(119.75)) + c * deg(mpf(131.849)) -
cap_F)))
+ flat_earth = (deg(-2235/1000000) * sin_degrees(cap_L_prime) +
+ deg(127/1000000) * sin_degrees(cap_L_prime -
cap_M_prime) +
+ deg(-115/1000000) * sin_degrees(cap_L_prime +
cap_M_prime))
+ extra = (deg(382/1000000) *
+ sin_degrees(deg(mpf(313.45)) + c * deg(mpf(481266.484))))
+ return beta + venus + flat_earth + extra
+
+
# see lines 192-197 in calendrica-3.0.errata.cl
-def `lunar_node(date):
+def `lunar_node(tee):
"""Return Angular distance of the node from the equinoctal point
- at fixed date, date."""
- return mod(moon_node(julian_centuries(date)) + deg(90), 180) - 90
+ at fixed moment, tee.
+ Adapted from eq. 47.7 in "Astronomical Algorithms"
+ by Jean Meeus, Willmann_Bell, Inc., 2nd ed., 1998
+ with corrections June 2005."""
+ return mod(moon_node(julian_centuries(tee)) + deg(90), 180) - 90
+
+def `alt_lunar_node(tee):
+ """Return Angular distance of the node from the equinoctal point
+ at fixed moment, tee.
+ Adapted from eq. 47.7 in "Astronomical Algorithms"
+ by Jean Meeus, Willmann_Bell, Inc., 2nd ed., 1998
+ with corrections June 2005."""
+ return degrees(poly(julian_centuries(tee), deg([mpf(125.0445479),
+ mpf(-1934.1362891),
+ mpf(0.0020754),
+ mpf(1/467441),
+ mpf(-1/60616000)])))
+
+def `lunar_true_node(tee):
+ """Return Angular distance of the true node (the node of the
instantaneus
+ lunar orbit) from the equinoctal point at moment, tee.
+ Adapted from eq. 47.7 and pag. 344 in "Astronomical Algorithms"
+ by Jean Meeus, Willmann_Bell, Inc., 2nd ed., 1998
+ with corrections June 2005."""
+ c = julian_centuries(tee)
+ cap_D = lunar_elongation(c)
+ cap_M = solar_anomaly(c)
+ cap_M_prime = lunar_anomaly(c)
+ cap_F = moon_node(c)
+ periodic_terms = (deg(-1.4979) * sin_degrees(2 * (cap_D - cap_F)) +
+ deg(-0.1500) * sin_degrees(cap_M) +
+ deg(-0.1226) * sin_degrees(2 * cap_D) +
+ deg(0.1176) * sin_degrees(2 * cap_F) +
+ deg(-0.0801) * sin_degrees(2 * (cap_M_prime -
cap_F)))
+ return alt_lunar_node(tee) + periodic_terms
+
+def `lunar_perigee(tee):
+ """Return Angular distance of the perigee from the equinoctal point
+ at moment, tee.
+ Adapted from eq. 47.7 in "Astronomical Algorithms"
+ by Jean Meeus, Willmann_Bell, Inc., 2nd ed., 1998
+ with corrections June 2005."""
+ return degrees(poly(julian_centuries(tee), deg([mpf(83.3532465),
+ mpf(4069.0137287),
+ mpf(-0.0103200),
+ mpf(-1/80053),
+ mpf(1/18999000)])))
+

# see lines 199-206 in calendrica-3.0.errata.cl
def `sidereal_lunar_longitude(tee):
@@ -4361,67 +4504,9 @@
return invert_angular(lunar_phase, phi, a, b)


-# see lines 3663-3732 in calendrica-3.0.cl
-def `lunar_latitude(tee):
- """Return the latitude of moon (in degrees) at moment tee.
- Adapted from "Astronomical Algorithms" by Jean Meeus,
- Willmann_Bell, Inc., 1998."""
- c = julian_centuries(tee)
- cap_L_prime = mean_lunar_longitude(c)
- cap_D = lunar_elongation(c)
- cap_M = solar_anomaly(c)
- cap_M_prime = lunar_anomaly(c)
- cap_F = moon_node(c)
- cap_E = poly(c, [1, mpf(-0.002516), mpf(-0.0000074)])
- args_lunar_elongation = \
- [0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 0, 4, 0,
0, 0,
- 1, 0, 0, 0, 1, 0, 4, 4, 0, 4, 2, 2, 2, 2, 0, 2, 2, 2, 2, 4,
2, 2,
- 0, 2, 1, 1, 0, 2, 1, 2, 0, 4, 4, 1, 4, 1, 4, 2]
- args_solar_anomaly = \
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, -1, -1, 1, 0,
1,
- 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1,
1,
- 0, -1, -2, 0, 1, 1, 1, 1, 1, 0, -1, 1, 0, -1, 0, 0, 0, -1, -2]
- args_lunar_anomaly = \
- [0, 1, 1, 0, -1, -1, 0, 2, 1, 2, 0, -2, 1, 0, -1, 0, -1, -1,
-1,
- 0, 0, -1, 0, 1, 1, 0, 0, 3, 0, -1, 1, -2, 0, 2, 1, -2, 3, 2,
-3,
- -1, 0, 0, 1, 0, 1, 1, 0, 0, -2, -1, 1, -2, 2, -2, -1, 1, 1,
-2,
- 0, 0]
- args_moon_node = \
- [1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1,
-1,
- -1, 1, 3, 1, 1, 1, -1, -1, -1, 1, -1, 1, -3, 1, -3, -1, -1, 1,
- -1, 1, -1, 1, 1, 1, 1, -1, 3, -1, -1, 1, -1, -1, 1, -1, 1, -1,
- -1, -1, -1, -1, -1, 1]
- sine_coefficients = \
- [5128122, 280602, 277693, 173237, 55413, 46271, 32573,
- 17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211,
- 2065, -1870, 1828, -1794, -1749, -1565, -1491, -1475,
- -1410, -1344, -1335, 1107, 1021, 833, 777, 671, 607,
- 596, 491, -451, 439, 422, 421, -366, -351, 331, 315,
- 302, -283, -229, 223, 223, -220, -220, -185, 181,
- -177, 176, 166, -164, 132, -119, 115, 107]
- beta = (deg(1/1000000) *
- sigma([sine_coefficients,
- args_lunar_elongation,
- args_solar_anomaly,
- args_lunar_anomaly,
- args_moon_node],
- lambda v, w, x, y, z: (v *
- pow(cap_E, abs(x)) *
- sin_degrees((w * cap_D) +
- (x * cap_M) +
- (y * cap_M_prime) +
- (z * cap_F)))))
- venus = (deg(175/1000000) *
- (sin_degrees(deg(mpf(119.75)) + c * deg(mpf(131.849)) +
cap_F) +
- sin_degrees(deg(mpf(119.75)) + c * deg(mpf(131.849)) -
cap_F)))
- flat_earth = (deg(-2235/1000000) * sin_degrees(cap_L_prime) +
- deg(127/1000000) * sin_degrees(cap_L_prime -
cap_M_prime) +
- deg(-115/1000000) * sin_degrees(cap_L_prime +
cap_M_prime))
- extra = (deg(382/1000000) *
- sin_degrees(deg(mpf(313.45)) + c * deg(mpf(481266.484))))
- return beta + venus + flat_earth + extra


+<<astronomical lunar calendars>>=
# see lines 3734-3762 in calendrica-3.0.cl
def `lunar_altitude(tee, location):
"""Return the geocentric altitude of moon at moment, tee,
@@ -4493,6 +4578,13 @@
return mt(385000560) + correction


+def `lunar_position(tee):
+ """Return the moon position (geocentric latitude and longitude [in
degrees]
+ and distance [in meters]) at moment, tee.
+ Adapted from "Astronomical Algorithms" by Jean Meeus,
+ Willmann_Bell, Inc., 2nd ed."""
+ return (lunar_latitude(tee), lunar_longitude(tee), lunar_distance(tee))
+
# see lines 3815-3824 in calendrica-3.0.cl
def `lunar_parallax(tee, location):
"""Return the parallax of moon at moment, tee, at location, location.
@@ -8028,6 +8120,8 @@
# TeX (use predefined rule)
pycalcal.tex: premarkup

+pycalcal.pdf: $(NW_MAIN:.nw=.tex) figure
+
# Python files (for pycalcal.py use predefined rule +
# extra dependency on premarkup)
pycalcal.py: premarkup
@@ -8061,7 +8155,8 @@
STATUS \
COPYRIGHT_DERSHOWITZ_REINGOLD \
makemake.sh \
-
+ figure.mp \
+ astro.mp

.PHONY : distro
distro: all