\ctxlua{metapost.showlog=true}
\startMPinclusions[+]
	picture rotation;
		theta:=60;
		phi:=10;
		def P primary c = ((redpart c)*cosd(theta)+(greenpart c)*sind(theta),((greenpart c)*cosd(theta)-(redpart c)*sind(theta))*sind(phi)+(bluepart c)*cosd(phi)) enddef;

		R = 2; % average distance from the center of Earth to center of Sun = 1 AU
		rS = 0.2; % drawn (innacurately) with the Sun's radius being 0.1 AU
		rE = 0.1; % drawn (innacurately) with the Earth's radius being 0.05 AU

		u:=216;

		picture earth;
		path equator,celestial;
		s:=15;

		% Celestial plane is ax+by+cz=0
		a:=-sind(23.5)*cosd(s);
		b:=-sind(23.5)*sind(s);
		c:=cosd(23.5);
		T:=81-366/360*(angle(b,-a)+180);

		% Construct an orthonormal basis {v,w} for the celestial plane
		color v,w;
		v:=(1-a*a,-a*b,-a*c);
		w:=(-a*b,1-b*b,-b*c);
		v:=v/sqrt(((redpart v)**2)+((greenpart v)**2)+((bluepart v)**2));
		w:=w-((redpart v)*(redpart w)+(greenpart v)*(greenpart w)+(bluepart v)*(bluepart w))*v;
		w:=w/sqrt(((redpart w)**2)+((greenpart w)**2)+((bluepart w)**2));

		% Celestial plane
		Nm:=1440; % Minutes in a day
		equator:=P(rE*u*v) for n=1 upto (Nm-1): .. P(rE*u*(v*cosd(n/Nm*360)+w*sind(n/Nm*360))) endfor .. cycle;
		celestial:=equator scaled (R/rE);
		path axis;
		axis:=(-1.5*rE*u*P(a,b,c))--(1.5*rE*u*P(a,b,c));
		picture equatorFront,equatorBack;
\stopMPinclusions
\starttext
	\startMPpage
		draw equator;
		clip currentpicture to ((axis rotated 90)--(-infinity*P(a,b,c))--cycle);
		equatorFront:=currentpicture;
		currentpicture:=nullpicture;
		draw equator;
		clip currentpicture to ((axis rotated 90)--(infinity*P(a,b,c))--cycle);
		equatorBack:=currentpicture;
		currentpicture:=nullpicture;
		draw axis;
		fill fullcircle scaled (2*rE*u) withcolor (0.6,0.8,1);
%		draw fullcircle scaled (2*rE*u);
		draw equatorFront;
		draw equatorBack dashed evenly scaled 0.5;
		earth:=currentpicture;
		currentpicture:=nullpicture;

		path ecliptic,eclipticCircle;
		Ne:=366; % Days in a (leap) year
		ecliptic:=P(R*u*cosd(-T/Ne*360),R*u*sind(-T/Ne*360),0) for n=1 upto (Ne-1): .. P(R*u*cosd((n-T)/Ne*360),R*u*sind((n-T)/Ne*360),0) endfor .. cycle;
		eclipticCircle:=(R*cosd(-T/Ne*360),R*sind(-T/Ne*360)) for n=1 upto (Ne-1): .. (R*cosd((n-T)/Ne*360),R*sind((n-T)/Ne*360)) endfor .. cycle;

		draw ecliptic withpen pencircle scaled 1.5bp;

		% Draw the Sun
		fill fullcircle scaled (2*rS*u) withcolor (0.95,0.8,0.15);

		% Draw the Earth
%		draw earth shifted point 80.70 of ecliptic; % March 21
%		draw earth shifted point 172 of ecliptic; % June 21
%		draw earth shifted point 265.66 of ecliptic; % September 22
%		draw earth shifted point 356.96 of ecliptic; % December 22

		rotation:=currentpicture;
		t:=1;
		draw earth shifted point t of ecliptic;
		d:=(a*xpart point t of eclipticCircle + b*ypart point t of eclipticCircle);
		draw (point t of ecliptic)--(d*u*P(a,b,c)) withpen pencircle scaled 1bp withcolor green;
		draw origin--(point t of ecliptic) withpen pencircle scaled 1bp withcolor red;
		draw origin--(d*u*P(a,b,c)) withpen pencircle scaled 1bp;
		fill fullcircle scaled 6bp;
		fill fullcircle scaled 6bp shifted (point t of ecliptic);
		fill fullcircle scaled 6bp shifted (d*u*P(a,b,c));
		setbounds currentpicture to ((unitsquare shifted (-0.5,-0.5)) xscaled 910 yscaled 350);

	\stopMPpage

	\dorecurse{365}{

		% March Equinox:     t = 81
		% June Solstice:     t = 173
		% September Equinox: t = 264
		% December Solstice: t = 356

		\startMPpage
			t:=\recurselevel + 1;
			draw  rotation;
			draw earth shifted point t of ecliptic;
			d:=(a*xpart point t of eclipticCircle + b*ypart point t of eclipticCircle);
%			show (t,d);
			draw (point t of ecliptic)--(d*u*P(a,b,c)) withcolor green;
			draw origin--(point t of ecliptic) withcolor red;
			draw origin--(d*u*P(a,b,c));
%			draw equator scaled (1/rE/(R++d)) shifted (point t of ecliptic) dashed evenly scaled 0.5;
			fill fullcircle scaled 6bp;
			fill fullcircle scaled 6bp shifted (point t of ecliptic);
			fill fullcircle scaled 6bp shifted (d*u*P(a,b,c));
			setbounds currentpicture to ((unitsquare shifted (-0.5,-0.5)) xscaled 910 yscaled 350);
		\stopMPpage
	}
\stoptext