Figure 0. Showing astronomic and geocentric coordinate systems.
Solar Ephemeris at 1990 Jan 0.0
-------------------------------
DAY0 726482 Day number of Epoch (origin approx. year BC 1)
M0 356.6349 Mean anomaly at epoch
L0 279.4033 Mean longitude along ecliptic
EQC1 1.9151 Equation of centre harmonic coefficient
EQC2 0.0200 " " " "
OBLQ 23.4406 Obliquity of ecliptic (inclination)
Step 1. Calculate number of days since epoch
--------------------------------------------
Enter Year (e.g.1985), Day number (Jan 1 = 1), hours, minutes, seconds:
YRDAY = INT( (YEAR-1)*365.25 ) - DAY0 Days from year difference
FRDAY = ( HOUR + ( MIN + SEC/60 )/60 )/24 Fraction of day
TIM = YRDAY + DAYNO + FRDAY Time in days since epoch
Note: These formulae only valid 1900 Mar 01 - 2100 Feb 28.
Step 2. Compute GHAA and True Longitude
---------------------------------------
MA = M0 + 360 * TIM / 365.2596413 Sun's Mean Anomaly
ML = L0 + 360 * TIM / 365.2421988 Mean Longitude
GHAA = ML-180 + 360 * FRDAY Greenwich Hour Angle Aries
TL = ML + EQC1*SIN(MA) + EQC2*SIN(2*MA) True Longitude
Note: 1. 'Longitude' means angle from ascending node up the ecliptic.
2. The Sun's argument of perigee (AP) is given by AP = ML - MA
Step 3a. Compute Unit Vector, Geocentric Rectangular Coordinates
----------------------------------------------------------------
Xsun = COS(TL) ( = COS(DEC) * COS(RA) )
Ysun = SIN(TL) * COS(OBLQ) ( = COS(DEC) * SIN(RA) )
Zsun = SIN(TL) * SIN(OBLQ) ( = SIN(DEC) )
Coordinate set is right-handed orthogonal located at the Earth's centre, X
and Y axes in the equatorial plane; X points at the ascending node, or the
first point in Aries (Vernal Equinox), and Y axis is at RA = 90 degrees;
Z axis normal to the plane, pointing approximately at the Pole Star.
Step 3b. Unit Vector, Geocentric, Greenwich Coordinates
-------------------------------------------------------
C = COS(GHAA): S= SIN(GHAA)
Xg = C * Xsun + S * Ysun ( = COS(LAT) * COS(-LON) )
Yg = -S * Xsun + C * Ysun ( = COS(LAT) * SIN(-LON) )
Zg = Zsun ( = SIN(LAT) )
Coordinate set is right-handed orthogonal located at the Earth's centre, X
and Y axes in the equatorial plane. X points at intersection of Greenwich
meridian and Equator, and Y axis points at longitude 90 degrees East. Z
axis normal to the plane, pointing up through North Pole. This set is as
the set above but rotated about the Z axis RAAN degrees East.
Step 4a. Astronomic Coordinates
--------------------------------
DEC = ASN (Zsun) Declination
RA = ATN (Ysun/Xsun) Right Ascension
if Xsun < 0 then RA = RA + 180 correct quadrant
if RA < 0 then RA = RA + 360
ATN is assumed to yield a value in the range +- 90 degrees. RA is in the
same quadrant as true longitude, TL.
Step 4b. Geographic Coordinates
-------------------------------
LAT = DEC Sun's Latitude, +North
LON = GHAA - RA Sun's Longitude, +West
Step 4c. Equation of Time
-------------------------
EQOT = ( ML - RA ) * 4 Minutes
Equation of time is the amount to be subtracted from 'sundial' time to
give local time, maxima being -(14m 16s) around Feb 12th, and +(16m 24s)
around Nov 3rd. The Sun is at maximum elevation at local noon minus EQOT.
Test Data
---------
The following example data are provided to help you to check your own
calculations. EXAMPLE: Evaluate all data for, Aug 12th 1985 at 0145:00
GMT. All angles reduced to range 0 - 360 degrees. From diary: DAYNO = 224
YRDAY = -1826: FRDAY = 0.072917: TIM = -1601.927083 days
MA 217.7751: ML 140.4681: TL 139.3144
GHAA 346.7181: RA 141.7354: DEC 15.0302
Xsun -0.758298: Xg -0.875425: LAT 15.0302 North
Ysun 0.598108: Yg 0.407897: LON 204.9828 West
Zsun 0.259328: Zg 0.259328: EQOT -(5m 4s)
Sources and Accuracy
--------------------
The solar data has been derived from:
1. The Astronomical Ephemeris 1984. ISBN:0-11-886919-1, HMSO & USGPO.
2. Explanatory Supplement to the Astronomical Ephemeris. 1961 Rev. 1974.
ISBN: 0-11-880578-9, HMSO (presently out of print, new edn. in prep.)
Angular errors are unlikely to exceed 0.2' of arc. The data may be used
for a decade forward or backwards from 1990 with graceful degradation.
These equations describe a 'mean' Sun's motion. Discrepancies compared
with an almanac, are created partly by planetary purturbations to the
Earth's orbit, partly by equating Universal time with Ephemeris time and
partly to numerical adjustments deliberately embodied in almanac tables.
It's the Earth/Moon centre of gravity and not the Earth's centre which
orbits the Sun. There are other sources too - consult the references.
SUN'S UP! Part 2
by J.R. Miller G3RUH
This article is a continuation of the first part (O.N 50, December 1984)
on the general theme of the relationship between the Sun and an Earth
satellite. I'm using Oscar-10 in particular, but the methods are in
general applicable to any of the others.
First a few notes about vectors and then how to calculate solar elevation
above the orbit plane, and the solar panel illumination coefficient.
Eclipses will be dealt with in part 3.
Vectors
-------
Before proceeding I think it would be useful to introduce a few things
about vectors. A vector is a neat way of describing a point or a line in
space, and vectors can be added, subtracted and multiplied, as can
ordinary numbers which are merely one-dimensional vectors.
In our context a vector consists of three numbers, representing
coordinates along the three axes. For example a point 1 unit along the
X-axis, and zero units along the Y and Z axes can be compactly written
< 1, 0, 0>.
A line at 45 degrees to the X and Y axes, and in the same plane is
< 0.71, 0.71, 0>
A line equidistant from all three axes is
< 0.58, 0.58, 0.58>, and so on. (Sketch them).
The length of a vector < X, Y, Z > is given by L = SQR( X*X + Y*Y + Z*Z ).
So you will see that all the above vectors are of length 1.0, and are not
suprisingly called Unit Vectors.
The angle between two unit vectors is found by what is called their scalar
(dot) product. Denoting unit vectors A and B by < XA, YA, ZA > and
< XB, YB, ZB> then their included angle ANG is found from:
Cos(ANG) = XA*XB + YA*YB + ZA*ZB = < A > dot < B > = A.B
The length of vector A projected onto unit vector B (or vice versa) is
also given by their scalar product, though in this case A need not
necessarily be of unit length. This is called 'resolving A onto B'.
If the scalar product of two unit vectors is unity, then the vectors are
coincident. If the scalar product of two vectors is zero then they are
at right angles (90 deg) to each other. Vectors also very obligingly obey
some familiar rules A+B = B+A, and A.B = B.A PROVIDED they are
expressed in the same coordinate system.
Verify these easy concepts for yourself using the three unit vectors
above, and inventing some of your own. For further reading try
"Advanced Engineering Mathematics", E. Kreyszig, John Wiley and Sons, ISBN
0-471-88941-5 or almost any similar title.
Solar Elevation above Oscar-10's Orbit
--------------------------------------
As a first example of using the sun's position let us work out the
elevation of the Sun above (or below) Oscar-10's orbit plane.
To do this we need a vector to describe the Sun's position, and another
vector to represent the satellite's orbital plane; a perpendicular to the
plane provides just this.
Then the scalar (dot) product of these two gives the cosine of the angle
'twixt them. We want the angle above the plane which is 90 minus this
(sketch it), i.e. the arcsine of their scalar product.
Denoting the orbit normal components by < Xnorm, Ynorm, Znorm >, the Sun's
vector as before by < Xsun, Ysun, Zsun >, and the solar elevation by SEL
we have:
SEL = ASN( Xsun*Xnorm + Ysun*Ynorm + Zsun*Znorm )
where for a plane, given its RAAN and inclination INCL the normal in
celestial coordinates is:
Xnorm = SIN(INCL) * SIN(RAAN)
Ynorm = - SIN(INCL) * COS(RAAN)
Znorm = COS(INCL)
(Verify this for yourself with some sketches, trying various values for
INCL and RAAN. For example if INCL = 0 then the normal vector becomes
< 0,0,1 >, ie straight up the Z-axis, and so on).
Continuing the example (Pt.1) for 1985 Aug 12th 0145:00 GMT, DATIM =
224.072917, Oscar-10 INCL = 25.6 deg and RAAN = 159.8 - 0.1724 * DATIM =
121.2 deg, we find:
Xnorm = 0.3696 Xsun = -0.7583 norm . Sun
Ynorm = 0.2238 Ysun = 0.5981 = 0.0874
Znorm = 0.9018 Zsun = 0.2593
and hence SEL = ASN(0.0874) = +5 deg.
(Fig. 1) When solar elevation approaches zero, as it does twice a year,
the sun lies in the satellite orbital plane, and hence Oscar-10 passes
through the earth's shadow on each orbit. Also the smallest angular
diameter of the Earth, seen from maximum range is +-ATN(6378/41000) = +-9
deg, so whenever the solar elevation is less than this (as in the above
example) an eclipse is inevitable at some time round the orbit.
You can also see that in 1987 the value hovers around zero for most of the
year, so we can expect some frequent changes to the operating schedule in
consequence.
Figure 1. OSCAR-10 SOLAR ELEVATION 1984-7. The maximum possible value
is the sum of Oscar-10 and the Sun's inclinations, i.e. 25.6 + 23.4 = 49 deg.
Illumination Coefficient
------------------------
From the satellite's point of view, perhaps the single most important
parameter is the solar panel illumination coefficient, for without
sunlight a non-nuclear powered spacecraft dies.
You will know from the pictures how Oscar-10's panels are arranged, and as
the satellite spins at 30 r.p.m. or so about the axis of symmetry each
panel is illuminated in turn.
Obviously for maximum effect the Sun needs to be perpendicular to the spin
axis, whilst if it becomes co-axial then there will be no illumination.
So to analyse this problem we need a vector representing the satellite's
attitude.
If we assume that Oscar-10 is nominally oriented, then the spin axis is
exactly parallel to the orbit semi major axis (pointing in the apogee to
perigee direction), so this gives us our attitude vector.
And, as before, the cosine of the angle between Sun and spin axis is given
by the scalar product of their two unit vectors. If we call this the Sun
Angle, then 0 deg = no illumination, 90 deg = 100% illumination.
We already have the Sun's vector, in celestial coordinates, so it remains
to translate a unit length along the spin axis into the same coordinate
set by rotating it through the various angles that define the orbit
attitude; argument of perigee (AP), inclination (INCL) and right ascension
of ascending node (RAAN). All elliptic orbit calculation programs do
this, though probably a bit obscured. The answer fully worked out is:
Xa = Cos(AP) * Cos(RAAN) - Sin(AP) * Cos(INCL) * Sin(RAAN)
Ya = Cos(AP) * Sin(RAAN) + Sin(AP) * Cos(INCL) * Cos(RAAN)
Za = Sin(AP) * Sin(INCL)
The sun angle (SA) and illumination coefficient (ILL) are thus:
SA = ACS( Xa*Xsun + Ya*Ysun + Za*Zsun)
ILL = 100 * SIN( SA ) (percent)
Worked example: Continuing as before for 1985 August 12 @ 0145 GMT. First
work out argument of perigee for DATIM = 224.072917:
AP = 336.6 + 0.2834*DATIM = 400.1 = 40.1 deg (see O.N. 50), and we already
have RAAN = 121.2 and INCL = 25.6 from above. Evaluating we find
Xa = -0.8931 Xsun = -0.7583 a . sun = 0.9608
Ya = 0.3534 Ysun = 0.5981 SA = ACS(0.9608) = 16 deg
Za = 0.2783 Zsun = 0.2593 ILL = 28%
A value of SA in the range 0 to 90 means that the Sun is illuminating
somwhere between the spacecraft's top and midriff, while a value
90-180 degrees, the lower half. Exactly 90 degrees is the ideal
condition, side-on.
Figure 2. OSCAR-10 NOMINAL ILLUMINATION computed as above for 1984 to 1987.
An illumination factor of 28% is unacceptably low, and the managers will
have had to take corrective action long before this happens, by
reorienting the spacecraft.
Now it is clearly desirable to keep the antennae pointing at the Earth as
far as possible, so the favoured option is to 'twist' the satellite in the
plane of the orbit only.
A twist of about +30 degrees has been in use now since Autumn 1984, and
has meant best access considerably before apogee. In terms of the
mathematics shown above, an in-plane twist of TW is modelled quite simply
by replacing argument of perigee AP in the equations by AP-TW !
Sketch it, and verify that an in terms of attitude an in-plane twist is
appears identical to an AP change.
Figure3. OSCAR-10 ILLUMINATION. The effect of a +30 degree twist and the
nominal curve. In particular note the improvement in illumination obtained
around the beginning of 1985. However you can also see that the +30 degree
twist will be no good for the period of our example, Autumn 1985.
Figure 4. OSCAR-10 ILLUMINATION showing the effect of a -30 degree adjustment
again superimposed on the nominal. Now the null in August 1985 is raised
to 50%, still not too good, but better than the original of less than 20%.
A reorientation like this will have to be implemented around June 1st.
This would also mean peak activity now concentrated after apogee, in much
the same way as is presently done before it. In addition, as I remarked
in O.N.49 page 3, eclipses are due in August 1985 around MA 85-130/256
when the transponder will have to be switched off. All in all an
interesting time ahead.
You might like to anticipate events of 1986-7, using the graphs.
Visualising these effects is made immeasurably easier with my simple 3-D
model (O.N.44 et seq) which comprises a hardboard ellipse and small
globe. It can be made in a couple of hours. A table lamp a few metres
away can be used to represent the Sun, and casts a shadow nicely for
eclipse studies. In part 3: ECLIPSES.
SUN'S UP! Part 3
by J.R. Miller G3RUH
This is the final part of a series in which I've looked at the
relationship between the Sun and a satellite. The first (O.N.50, December
1984) described a proper Sun model, whilst in the second (O.N.51, February
1985) vectors were explained and the results used to investigate sun
angles and their effect on a satellite's solar illumination.
In this part I deal with eclipses. Though I'm using Oscar-10 for my
examples the methods apply to any Earth satellite, low orbiters like
UOSAT, RS and NOAA through to OSCAR-10 and the geostationary comsats.
Satellite Eclipses
------------------
A satellite is said to be 'eclipsed' when it is within the Earth's shadow,
and the term 'shadowed' means the same thing. You can visualise what is
happening in figure 5, which shows the shadow zone and a satellite within
it.
Figure 5. Illustrating conditions of a satellite's solar eclipse. It should be obvious that the only condition necessary for an eclipse is that the distance marked UD (Umbral Distance) should be less than the radius of the Earth, since this is the size of the shadow. From the triangle, distance UD = Rsat * Sin(UA), so we need to work out Rsat and the Umbral Angle UA. Remembering the notes about vectors (pt. 2), the cosine of UA is simply the scalar (dot) product of a unit vector AWAY from the Sun, and a unit vector towards the satellite, i.e., denoting Cos(UA) by CUA: CUA = - (Xsun*Xsat + Ysun*Ysat + Zsun*Zsat) Since it's possible to get UD less than one Earth radius but on the SUNNY side of Earth, we must also check that the satellite is on the dark side, i.e. UA less than 90 degrees which means that CUA must be positive. Sun and Satellite's Unit Position Vector ---------------------------------------- The Sun unit vector calculations were dealt with in part 1 of this series. For the satellite we have to do the full orbital calculation, i.e. find mean anomaly (MA), thence eccentric anomaly (EA) and the in-plane coordinates (X and Y). These are then rotated through argument of perigee (AP), inclination (INCL) and then right ascension of the ascending node (RAAN) until we arrive at the satellite's position vector in celestial coordinates, the same as we used for the Sun. First solve Kepler's equation (in radians) MA = EA - EC*Sin(EA) for EA given MA where EC is eccentricity. This equation is 'transcendental' which means it's a sod, and has to be resolved iteratively (see appendix). Then for compactness in what follows denote: CEA = Cos(EA): SEA = Sin(EA) CAP = Cos(AP): SAP = Sin(AP) CNC = Cos(INCL): SNC = Sin(INCL): where INCL is Inclination. CRA = Cos(RAAN): SRA = Sin(RAAN) Rsat = SMA*(1-EC*CEA): where SMA = Semi Major Axis. X = (CEA-EC)/(1-EC*CEA): Unit satellite vector in Y = SQR(1-EC*EC)*SEA/(1-EC*CEA): orbit plane. Xsat = X*(CAP*CRA - SAP*CNC*SRA) + Y*( -SAP*CRA - CAP*CNC*SRA) Ysat = X*(CAP*SRA + SAP*CNC*CRA) + Y*( -SAP*SRA + CAP*CNC*CRA) Zsat = X*(SAP*SNC) + Y*( CAP*SNC) These calculations are embodied in all elliptic orbit programs (they use longitude of ascending node rather than RAAN because they work in geocentric coordinates). Equations above, though correct are inefficient and in the appendix they are presented more nicely. Eclipse - Worked Example ------------------------ Continuing the example started in parts 1&2, for Oscar-10 on 1985 August 12 at 0145 GMT. If you flog through the arithmetic you will find: MA = 129.3, AP = 40.1, RAAN = 121.2 (all in degrees - see pt1. page 2). Since EC = 0.61 we find EA = 147.9 degrees, and using a semi major axis of SMA = 26100 km you get, in addition to a headache: Xsat = 0.7864 Xsun = -0.7583 (as before. part 1 test data) Ysat = -0.5923 Ysun = 0.5981 Zsat = -0.1755 Zsun = 0.2593 Rsat = 39582 km CUA = 0.9961 UA = 5.09 deg UD = 3510 km The Earth's (equatorial) radius is 6378 km and, since the umbral distance is only 3510 km, at 0145 GMT August 12 1985 Oscar-10 will be eclipsed. Observability ------------- An MA of 129 degrees, (91/256 of the orbit) is 100 minutes before apogee, and the event will be observable for an hour or so (on the telemetry) by stations in Europe, Africa, the middle East and South America. The transponder will not, of course, be switched on! This 'eclipse season' lasts from 10th August to 3rd September 1985,and is illustrated on page 3 of O.N.49, Oct 1984 with an eclipse diagram. OSCAR-10 Planning Software -------------------------- These ideas and others are embodied in an updated version of the program which appeared as a pull-out centrefold in O.N. 45 (December 1983). Written in BBC BASIC and called PLAN10 it's available from AMSAT-UK. [ Updated 1990 Oct ]. It is well commented and serves as a useful comprehensive and comprehensible set of reference equations for solving the elliptic orbit problem, WITH DRAG. It can be dissected and refashioned to individual needs. It can be translated into other BASIC dialects. As coded, a sample PLAN10 output looks typically like the following, though as is usual with BASIC programs, the time increment and most else can be changed at will. You will notice the data is the same as for the examples in this series. SATELLITE PREDICTIONS --------------------- Location: G3RUH Lat 52.208 Lon 0.059 Ephemeris Dated: 1985 Day 224 OSCAR-10 YEAR: 1985 ORBIT: 1627 ILL: 28% SOLAR EL: 5 Deg DAY:HHMM MA/256 MODE RANGE EL AZ SQ UMD ECL? -------------------------------------------------------------- 224:0100 75 B 34114 19 203 25 1.07 + 224:0115 80 B 35320 18 204 22 0.85 ECL 224:0130 86 B 36396 17 205 20 0.66 ECL 224:0145 91 B 37348 16 206 17 0.55 ECL 224:0200 97 B 38178 14 207 15 0.58 ECL 224:0215 102 B 38891 13 208 13 0.73 ECL 224:0230 108 B 39488 12 209 12 0.94 ECL 224:0245 113 B 39971 10 210 10 1.18 + etc Wither? ------- If you have found this little series useful - or found it induced terminal boredom, I should like to know. I am quite happy to scribble more on the numerical aspects of satellites - but you may not be interested! Soon I propose to write about drag, and maybe about the non-spherical Earth. Perhaps you could suggest topics. Please don't forget the SAE though. [ I got 2 replies. One said that the article was "mind-blowing", and the other that I couldn't spell whither. JRM ] APPENDIX - Satellite Position Calculation ----------------------------------------- Equations are given in BASIC, as they've been tested in that form. By separating out the coordinate transformations the number of multiplies and additions are much reduced. The test at line 2140 is not set to umpteen decimal places. Newton's method is a 'second order' algorithm, so as coded will converge to within 0.002*0.002/2 radians, (for Oscar-10, this is equivalent to about 50 metres worst case position error, vastly better than the ephemeris accuracy). For definitions see part 1. 2000 REM Calculate Satellite Position 2010 REM ---------------------------- 2020 REM Assumes mean anomaly (MA), Argument of perigee (AP) and Right 2030 REM Ascension of Ascending Node (RAAN) to have been calculated 2040 REM as well as inclination (INCL), eccentricity (EC), semi-major axis 2050 REM (SMA) and semi-minor axis (SMB). All angles in RADIANS. 2060 2070 REM Solve MA = EA - EC*SIN(EA) for EA given MA. (Newton's Method) 2080 EA = MA: REM Initial solution 2090 REPEAT 2100 C = COS(EA): S = SIN(EA) 2110 DNOM = 1-EC*C 2120 DE = (EA - EC*S - MA)/DNOM: REM Change to EA for better 2130 EA = EA - DE: REM solution by this amount 2140 UNTIL ABS(DE) < 0.002: REM until converged. Coding as 2150 REM fast as possible. 2160 Rsat = SMA*DNOM: REM Geocentric distance 2170 2180 REM Now compute satellite coordinates in plane of ellipse 2190 C = COS(EA): S = SIN(EA) 2200 X1 = SMA*(C - EC) 2210 Y1 = SMB*S: REM Note: SMB=SMA*SQR(1-EC*EC) 2220 Z1 = 0 2230 2240 REM Rotate coordinates about Z through Argument of perigee 2250 C= COS(AP): S= SIN(AP) 2260 X2 = X1*C - Y1*S 2270 Y2 = X1*S + Y1*C 2280 Z2 = Z1 2290 2300 REM Rotate coordinates about X through Inclination 2310 C=COS(INCL):S=SIN(INCL) 2320 X3 = X2 2330 Y3 = Y2*C 2340 Z3 = Y2*S 2350 2360 REM Rotate about Z through RA of Ascending Node, and normalise 2370 C = COS(RAAN): S = SIN(RAAN): REM with respect to Rsat 2380 Xsat = (X3*C - Y3*S)/Rsat 2390 Ysat = (X3*S + Y3*C)/Rsat 2400 Zsat = (Z3 )/Rsat 2410 2420 REM Satellite's position unit vector in celestial coordinates 2430 REM now calculated. For geocentric coordinates use RAAN-GHAA 2440 REM above, instead of RAAN.
