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

Notice that face 1 is not seen the planet, but all the four lateral faces get some radiation.

>

Qpi:=alpha*'A'*Fip*epsilon*sigma*Tp^4;F1p:=0;F2p:=(1/Pi)*(arctan(1/sqrt(h^2-1))-sqrt(h^2-1)/h^2);F3p:=1/h^2;F4p:='F2p';F5p:='F2p';F6p:='F2p';Qp:='alpha*A*(F1p+F2p+F3p+F4p+F5p+F6p)*epsilon*sigma*Tp^4';F2p_LEO:=evalf(subs(LEOh,F2p));F3p_LEO:=evalf(subs(LEOh,F3p));Qp_LEO:=subs(dat,evalf(subs(LEOh,eqAface,dat,Qp)));F2p_GEO:=evalf(subs(GEOh,F2p));F3p_GEO:=evalf(subs(GEOh,F3p));Qp_GEO:=subs(dat,evalf(subs(GEOh,eqAface,dat,Qp)));

i.e. IR inputs are 127 W at LEO and 1.5 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:

>	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);

The view factors from a face to the planet are already computed (Fip).

>

Qai:='alpha*A*Fip*rho*E*Fa';Qa1:=0;Qa2:=subs(Fip=F2p,Qai):Qa3:=subs(Fip=F3p,Qai):Qa4:=subs(Fip=F4p,Qai):Qa5:=subs(Fip=F5p,Qai):Qa6:=subs(Fip=F6p,Qai):Qa:='Qa1+Qa2+Qa3+Qa4+Qa5+Qa6';plot(subs(LEOh,LEO_es,eqAface,dat,SI0,[Qa1,Qa2,Qa3,Qa4,Qa5,Qa6,Qa]),phi=0..2*Pi,Qa_W_m2=0..250);N:=100:Qa_max_:=max(seq(evalf(subs(LEOh,LEO_es,phi=2*Pi*i/N,eqAface,dat,SI0,Qa)),i=1..N))*W_;plot(subs(GEOh,GEO_es,eqAface,dat,SI0,[Qa1,Qa2,Qa3,Qa4,Qa5,Qa6,Qa]),phi=0..2*Pi,Qa_W_m2=0..3);N:=100:Qa_max_:=max(seq(evalf(subs(GEOh,GEO_es,phi=2*Pi*i/N,eqAface,dat,SI0,Qa)),i=1..N))*W_;

i.e. the satellite gets a maximum of 223 W when at the subsolar point in LEO. 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+Q, 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-6*A*sigma*T(t)^4:eqEB_:=evalf(subs(LEOh,LEO_es,eqAface,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:=9;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,eqAface,dat,SI0,Dt*(Qs+Qa+Qp-6*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 262 K and 275 K with a period average of 269 K (i.e. between -11 ºC and 2 ºC, with a mean of -4 ºC).

In the case of GEO

>

To:=subs(SI0,rhs(T_GEO));N:=5;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,GEOh,GEO_es,eqAface,dat,SI0,Dt*(Qs+Qa+Qp-6*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,GEOh,GEO_es,eqAface,dat,SI0,eqEB)),T(0)=300],type=numeric,range=0..N*To);plRunge:=odeplot(sol_,color=black):display([plEuler,plRunge]);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,200..300]);

Notice how pronunced the effect of eclipse may be on geostationary orbits (fortunately, there are only GEO eclipses near the equinocces, contrary to LEO, where eclipse duration varies little with season).

>