If the planet is assumed isothermal, the infrared input is constant along the orbit. The view factor from a small sphere to a much larger sphere is found in View factor tabulations.

>	Qp:=alpha*A*Fbp*epsilon*sigma*Tp^4;A:=Pi*D^2;Fbp:=(1-sqrt(1-1/h^2))/2;eqh:=h=(H+R)/R;LEOh:=subs(H=H_LEO,dat,eqh);GEOh:=subs(H=H_GEO,dat,eqh);Fbp_LEO:=subs(LEOh,Fbp);Fbp_GEO:=subs(GEOh,Fbp);Qp_LEO:=subs(dat,evalf(subs(alpha=1,LEOh,dat,Qp)));Qp_GEO:=subs(dat,evalf(subs(alpha=1,GEOh,dat,Qp)));

i.e. IR inputs are 259 W at LEO and 4 W at GEO.

e) Albedo input along the orbit, assuming an albedo of 0.3 and a simple albedo model.

Maximum albedo input occurs at subsolar position and is Qa_max:=alpha*A*Fbp*rho*E.

For the variation of albedo input with orbit angle, a simple model may be that it is proportional to the lit area seen (i.e. the crescent) squared (to account for the fact that the smaller the crescent, the smaller the irradiance got per normal surface). However, as the fraction of lit area seen depends on orbit altitude, a simpler model may be to compute the crescent area from afar (i.e. its orthogonal proyection) and multiply it (its square) by a parabolic function that ends at the eclipse; namely:

>

Qa_max:=alpha*A*Fbp*rho*E;Qa_max_LEO:=subs(dat,evalf(subs(alpha=1,LEOh,dat,Qa_max)));Qa_max_GEO:=subs(dat,evalf(subs(alpha=1,GEOh,dat,Qa_max)));Fcrescent:=(1+cos(phi))/2;Fa:=Fcrescent^2*piecewise(p(phi)<phi_es,(1-(p(phi)/phi_es)^2),p(phi)>(2*Pi-phi_es),(1-((2*Pi-p(phi))/phi_es)^2),0);Qa:='Qa_max*Fa';plot(subs(alpha=1,LEOh,LEO_es,dat,SI0,Qa),phi=0..2*Pi,Qa_W_m2=0..500);plot(subs(alpha=1,GEOh,GEO_es,dat,SI0,Qa),phi=0..2*Pi,Qa_W_m2=0..10);

Notice that planet effects (IR and albedo) can be neglected at very high orbits like GEO.

f) Periodic temperature evolution, assuming the body is isothermal and has a mass of 100 kg with an average thermal capacity of 1000 J/(kg·K)..

The energy balance is dE/dt=W+Q=0+Qs+Qa+Qp-Qout, since there is no mechanical or electrical work. For instance, in the LEO case:

>	eqEB:=dE/dt=Q[net,inp]-Q[net,out];eqEB:=m*c*diff(T(t),t)=Qs+Qa+Qp-A*sigma*T(t)^4:eqEB_:=subs(alpha=1,LEOh,LEO_es,Afrontal=A/4,dat,SI0,eqEB);

a first-order ordinary differential equation which must be solved with some initial conditions (e.g. T(0)=300 K) until transients decay and a periodic solution remains.

This can be done by Euler's method or better by some Runge-Kutta method.

Let us run the numerical simulation for N=5 orbits, with Dt=100 s, starting with T0=300 K. The number of discrete t-points are N*Tperiod/Dt.

>

To:=subs(SI0,rhs(T_LEO));N:=5;Dt:=100;TE[1]:=300;Np:=round(N*To/Dt);for i from 1 to Np-1 do TE[i+1]:=TE[i]+evalf(subs(phi=2*Pi*Dt*i/To,alpha=1,LEOh,LEO_es,Afrontal=A/4,dat,SI0,Dt*(Qs+Qa+Qp-A*sigma*TE[i]^4)/(m*c)));od:plEuler:=plot([seq([Dt*i,TE[i]],i=1..Np)]):sol_:=dsolve([evalf(subs(phi=2*Pi*t/To,SI0,eqEB_)),T(0)=300],type=numeric,range=0..N*To);plRunge:=odeplot(sol_,color=black):display([plEuler,plRunge]);

We see the 5th orbit is already periodic and plot it aside (notice that the initial and final T-values must coincide to be periodic, but this is not the mean temperature, Tmean.

>	i:='i':Tmean:=sum(TE[i],i=round(Np*(N-1)/N)..Np)/(Np+1-round(Np*(N-1)/N));Tnax:=max(seq(TE[i],i=round(Np*(N-1)/N)..Np));Tnin:=min(seq(TE[i],i=round(Np*(N-1)/N)..Np));plmean:=plot([[0,Tmean],[N*To,Tmean]]):display([plEuler,plRunge,plmean],view=[(N-1)*To..N*To,250..300]);

i.e, once in periodic state, the body temperature oscillates between 268 K and 284 K with a period average of 276 K (i.e. between -5 ºC and 11 ºC, with a mean of 3 ºC).

In the case of GEO

>

To_:=subs(SI0,rhs(T_GEO));N_:=3;Dt_:=1000;TE_[1]:=300;Np_:=round(N_*To_/Dt_);for i from 1 to Np_-1 do TE_[i+1]:=TE_[i]+evalf(subs(phi=2*Pi*Dt_*i/To_,alpha=1,GEOh,GEO_es,Afrontal=A/4,dat,SI0,Dt_*(Qs+Qa+Qp-A*sigma*TE_[i]^4)/(m*c)));od:plEuler_:=plot([seq([Dt_*i,TE_[i]],i=1..Np_)]):sol_:=dsolve([evalf(subs(phi=2*Pi*t/To_,alpha=1,GEOh,GEO_es,Afrontal=A/4,dat,SI0,eqEB)),T(0)=300],type=numeric,range=0..N_*To_);plRunge:=odeplot(sol_,color=black):display([plEuler_,plRunge]);i:='i':Ni_:=round(Np_*(N_-1)/N_);Tmean:=sum(TE_[i],i=Ni_..Np_)/(Np_+1-Ni_);Tnax:=max(seq(TE_[i],i=Ni_..Np_));Tnin:=min(seq(TE_[i],i=Ni_..Np_));plmean:=plot([[0,Tmean],[N_*To_,Tmean]]):display([plEuler_,plRunge,plmean],view=[(N_-1)*To_..N_*To_,240..300]);

Notice how pronunced the effect of eclipse may be (and how much Euler's solution departs from the more accurate R-K) on geostationary orbits (fortunately, there are only GEO eclipses near the equinocces, as wxplained above, contrary to LEO, where eclipse duration varies little with season).

>

Comparison with the linear sinusoidal approximation for LEO: Tlin=275+10.8*cos(phi-1.32) in [K], corresponding to a sinusoidal input (Qs_max_+Qa_max_LEO)*(1+cos(phi))/2+Qp.

>

TlinLEO:=275+10.8*cos(phi-1.32);plEuler_:=plot([seq([2*Pi*(Dt*i-(Np-54)*Dt)/5429,TE[i]],i=Np-54..Np)],color=green):plin:=plot(TlinLEO,phi=0..2*Pi,'T'=250..300):display([plEuler_,plin],view=[0..2*Pi,250..300]);plot(subs(alpha=1,LEOh,LEO_es,dat,SI0,[(Qs_max_+Qa_max_LEO)*(1+cos(phi))/2+Qp,(Qs+Qa+Qp)]),phi=0..2*Pi,Q=0..2000);

The great advantage of this approximation is that it can be analytically solved.

>