当前位置:首页 >> 能源/化工 >>

The Astronomical Almanac's algorithm for approximate solar position (1950–2050)


Solar Energy Vol, 40, No. 3. pp. 227-235, 1988
Printed in the U.S.A.

0038-092X/X8 $3.0~) + .00 Copyright ~ 1988 Pergamon Journals Ltd.

THE ASTRONOMICAL ALMANAC'

S ALGORITHM FOR APPROXIMATE SOLAR POSITION (1950-2050)
JOSEPH J. MICHALSKY Atmospheric Sciences Research Center, State University of New York, Albany, NY 12222, U.S.A. Abstract--This article implements an algorithm for the calculation of solar position that has a stated. and partially demonstrated, accuracy of 0.01 deg until the year 2050. The algorithm is taken from The Astronomical Almanac, which has published it as an addendum to their very accurate tabulations since 1984. It uses the same approach as an earlier paper by Walraven, but has a more simplified form. This subject of the calculation of the solar position is revisited because of the somewhat confusing sequence of notes and letters that followed Walraven's publication in this journal and because of a report challenging the accuracy claimed by him. The Almanac approximation is compared to exact values and to some commonly used approximations. An appendix contains an annotated subroutine in FORTRAN, which should be easily converted to other programming languages. I. INTRODUCTION The solar literature does not have a shortage of papets on the subject of solar position calculation. The motivation for yet another one is proposed in the next paragraphs, Many formulas have appeared in papers with only vague comments regarding their accuracy. Typically, a statement to the effect that the calculations are "adequate for most engineering applications" accompahies the presentation of same. Indeed, these formulas may be well within the accuracy requirements for some applications. As an example, flat-plate tracking photovoltaic panels do not require extremely accurate tracking. Even if the panels are tilted 5 ° with respect to the normal to the sun their output is within 0.4% of the maximum response to the direct component, The actual response may not improve even that much for perfect alignment since a significant portion of the radiation received by the panel can be skylight and reflected radiation, which may not be affected significantly by realignment. On the other hand, there are stringent tracking requirements for some applications. Solar thermal power systems with high concentration ratios require accurate tracking. Vant-Hull and Hildebrandt[l] used 1 mrad (0.06 °) in their calculations, however, the more accurate the pointing, the better the efficiency up to the mechanical limitations of the tracking mechanism, Walraven[2] published an algorithm almost a decade ago with a stated accuracy well within this limit, His algorithm differs from most in that it does not require the sun to have the same position at the same time each year. Rather, it is based on the sun's astronomical position with respect to a given point in time. Walraven derived simplified forms for the high precision calculations used in The American Ephemeris and Nautical Almanac[3]. The latter formulas are accurate to 0.1 arcseconds. His formulas have a stated accuracy of 0.01 °. Why then should there be interest in writing or even reading another solar position paper? There are 227 two reasons. The primary one has to do with the response evoked by Walraven's paper. At least 10 notes and letters to Solar Energy followed its publication[4-13]. Some provided what can be considered useful changes, and others did not. Each time a user implements the solar position calculation, it is necessary to read 11 contributions to decide which changes are meaningful. My initial plan in addressing this subject was to send a brief note to Solar Energy with Walraven's subroutine slightly revised to incorporate the changes that I considered useful. The second reason to address this subject was a report by Zimmerman[14]. In that report, Walraven's accuracy was disputed. The claim was that the accuracy was 0.06 ° in azimuth and 0.0130 in elevation between 1979 and 1986. The implication was that it was less accurate outside this time frame. Since 1986 has passed, I considered submitting a note that would use Walraven's algorithm, but would base the approximation on a more current reference year. In attempting to reproduce Walraven's results, I consulted a recent copy of The Astronomical Almanac[15] (hereafter, The Almanac). Beginning with the 1984 issue of that reference (this title replaced The American Ephemeris and Nautical Almanac[3] in 1981), I found that a set of formulas are published as an addendum to the very accurate tabulated data. These formulas have an accuracy of 0.01 ° for the century that began in 1950. The Almanac's algorithm for the calculation of approximate solar position is presented in the next section. It closely resembles Walraven's[2], but has a simpler form. The coordinate systems' angular ranges differ slightly from those adopted by Walraven, but are more representative of the conventions used by astronomers, when referring to astronomical quantities, and by meteorologists, when referring to the local coordinates. A few ancilliary calculations are added for those who may want to use them. The third section contains comparisons with a few of the more commonly used expressions for solar position. These should allow one to decide whether certain approx-

228

J.J. MICHALSKY Note that the specification (0 - L < 360 °) implies that the quantity L is to have multiples of 360 ° added or subtracted until is has a value in the stated range. 2.3 Celestial coordinates To calculate the right ascension and declination we use the following formulas: tan(ra) = cos(ep) × sin (l)/cos(l) sin(dec) = sin(ep) × sin(I). The formula for the right ascension is written in this form (note that sin(l)/cos(l) = tan(l)) because we wish to assign the right ascension to an angle between 0 and 360 °. Having the signs of the numerator and denominator allows us to accomplish this. A twoargument inverse tangent function eases the task somewhat if one is available in the particular language of the subroutine. Although FORTRAN has this function, I have chosen to demonstrate how one derives the correct angle with a single argument inverse tangent function for use in other programming languages. 2.4 Local coordinates In order to convert to azimuth and elevation, we need to know the declination and hour angle. In order to calculate the hour angle, we need to know the local mean sidereal time (Imst) and the right ascension. The Imst can be defined as the right ascension of the point on the meridian at the celestial equator. With this definition, the hour angle is found by subtracting the right ascension from the lmst. The Almanac tabulates Greenwich mean sidereal time (gmst) at midnight and the formula to approximate gmst at an arbitrary time is given by gmst = 6.697375 + 0.0657098242 x n + hour(UT) (0 -< gmst < 24 h).

imations are adequate for their applications. In particular, Walraven's[2] algorithm is compared to The Almanac algorithm to test whether there is a need to implement the new algorithm for those currently using the one from Walraven. Comments summarizing the accuracy of the various algorithms are offered in the last section of the paper. An appendix completes the article. It contains a FORTRAN subroutine with ample commentary that should allow straightforward implementation in other languages, 2. THE ALMANAC ALGORITHI~I One of the most difficult aspects in describing or understanding the movement of the sun is in understanding the various coordinate systems within which its position is specified. These are described verbally and graphically in Walraven's paper[2] and need not be repeated here. 2.1 Time In The Almanac algorithm, the time argument is the difference in days between the current Julian date, jd, and JD 2451545.0, which is noon 1 January 2000 Universal Time (UT). The current Julian date may be specified directly or calculated from more familiar expressions of time using jd = 2432916.5 + delta x 365 + leap + day + hour/24, where delta = year - 1949, and leap = integer portion of (delta/4). Day is the day of the year (e.g., Feb 1 = 32), and hour is the UT in hours including fractions thereof. The leading number in the formula forjd is the Julian date for midnight 0 January 1949 UT. The choice of 1949 is arbitrary, but it does allow one to work with positive numbers in the calculation of jd for the century in question. The details of specifying local time in terms of UT and day of the year are left to the reader to include in the code that calls the subroutine,
?

The variable "hour" is not multiplied by a factor to convert it to sidereal hours because n already includes the fractional day and this correction is made in that product. The east longitude of the site in question is added to the gmst in order to specify lmst. Since longitude is usually expressed in degrees, the longitude should be divided by 15 to convert to hours. East longitude is considered positive. Therefore, lmst = gmst + (east long.)/15 (0 --< lmst < 24 h),

2.2 Ecliptic coordinates The ecliptic coordinates are calculated from the time argument according to the following steps: n = jd - 2451545.0

(2)

L (mean long.) = 280.460 + 0.9856474 x n (0 -< L < 360*) g (mean anomaly) = 357.528 + 0.9856003 × n (0 -< g < 360 °) l (ecliptic long.) = L + 1.915 × sin (g) + 0.020 x sin (2 × g) ep (obliquity of the ecliptic) = 23.439 - 0.0000004 x n (degs) (0 ---< e < 360 °)

Algorithm for approximate solar position and the hour angle (ha) is then given by ha = lmst - ra ( - 12 < ha -< 12 h). (3)

229

from the earth. In terms of the astronomical unit, which is the mean sun-earth distance and assigned a value 1, the distance varies according to R = 1.00014 - 0.01671 x cos(g) - 0.00014 x cos(2 × g). The radius is then calculated from solar radius = 0.2666/R. These corrections apply to only one particular structure of the atmosphere, which it seldom, if ever, assumes exactly. It is impractical to expect an accuracy of 0.01 ° when observing at low solar elevations because of this uncertainty. There has been no mention of the equation of time in this discussion. The equation for ecliptic longitude (1) accounts for the differences in the length of the solar day that gives rise to this equation. For completeness, the equation of time can be calculated from the previously derived variables using EOT = L(mean longitude) - ra(right ascension). Frequently, solar researchers wish to know the solar position in order to find the sun's angle of incidence on an arbitrarily oriented flat surface. Rewriting lqbal's[18] eqn 1.6.5a, we use the calculated variables to define cos (aoi) = sin (el)
×

Before the sun reaches the meridian, its hour angle is negative, and its hour angle is positive in the postmeridian hemisphere. (Note that many solar researchers use the opposite sign convention for ha with the eastern hemisphere positive and western hemisphere negative.) To calculate the elevation and azimuth, we use the following transformations: sin (el) = sin (dec) × sin (lat) + cos (dec) x cos (lat) x cos (ha) (degs) and sin (az) = - c o s (dec) x sin (ha)/cos (el) (0 -< az < 360 °)

(4)

(5)

where lat is the latitude of the site and hour angle has been converted to a more conventional angular argument. Azimuth is measured from north, where it is 0 °, through east to 360 °. (Again, note that many in solar energy use the convention that due south is 0 ° azimuth and east is positive while west is negative azimuth.) In order to assign the correct quadrant, a little care must be taken. Walraven[2] showed that at an azimuth of 90 ° sin(el) = sin(dec)/sin(lat), We will call it elc for el critical. To correctly assign azimuth we use the algorithm that if el -> elc, and if el -< elc then and ha > 0, then az = 180 degs - az,

cos(slope) + cos(el)

x sin(slope) x cos(ra - saz). where "slope" is the angle that from the horizontal, ~saz" is the fined similarly to local azimuth, gle of incidence with respect to the surface is tilted surface azimuth deand "aoi" is the anthe surface normal.

3. COMPARISONS OF APPROXIMATIONS 3.1 True Almanac values vs Almanac approximations We cannot check our calculations using azimuth and elevation since these quantities are not tabulated in The Almanac. The quantities needed to calculate azimuth and elevation are the hour angle, the declination, and the latitude (eqns (4) and (5)). Hour angle is calculated from fight ascension and local mean sidereal time (eqn (3)). Local sidereal time is obtained from the longitude and the Greenwich mean sidereal time (eqn (2)). The variables in the calculation of azimuth and elevation for a fixed location then are the fight ascension, the declination, and the Greenwich mean sidereal time. These are tabulated in The Almanac for 0 h UT for each day of the year. To check the accuracy of the approximations, we will compare the tabulated values in The Almanac to the nearest 0.1 s of time or nearest arcsecond of angle. These comparisons will be for 1986.

az = 360 + az.

2.5 Additional considerations The calculations above are for the center of the sun, which is about 1/2 ° in diameter and for a nonrefracting atmosphere. When observing near the horizon, refraction is important, and when calculating sunrise and sunset the angular extent of the sun needs to be taken into account, since the upper limb (or edge) of the solar disk is used as the reference point, The Almanac's refraction correction is given by r = 3.51561 x (0.1594 + 0.0196 x el + 0.00002 × el^2)/(l + 0.505 x el + 0.0845 x el^2), where 3.51561 = 1013.2 m b / 2 8 8 . 2 K. These are the atmospheric pressure and temperature for the U.S. Standard Atmosphere[16]. To calculate the diameter of the sun at any time, we need to know its distance

230
O O

J . J . MICHALSKY

_6
m × o.

O O

go
E r.4

I
u

0d
O O

I
f_

o
o

~

~

[

0

~.o0

200 Day o f 1.986

300

400

Fig. 1. The tabulated Almanac declination minus The Almanac's approximate declination as a function of the day of 1986. In Fig. 1 the difference of the declination tabulated in The A l m a n a c and the declination calculated from The A l m a n a c ' s lower precision algorithm is plotted as a function of the day of the year in 1986. The true A l m a n a c declination is recorded to the nearest arcsecond, which could, at most. affect the difference by about 0.0001 °, about the size of the high frequency noise in this plot. The approximation is well within 0.01 ° of the true declination. In Fig. 2 the same kind of comparison is made for right ascension. In this case, the data are in hours since that
?Y O 0 0 Q t

is the normal unit for right ascension. In this plot, the true Ahnanac right ascension is rounded to the nearest 0.1 s. To convert to degs, multiply by 15. The maximum error is about 0.009 °. The error is not symmetric about zero as it would be if the approximation were determined using a least squares procedure. In any case, the error is better than the 0.01 ° accuracy claimed for The A l m a n a c approximation. In Fig. 3, true Greenwich mean sidereal time to the nearest 0.1 s is differenced with the approximated value in hours. Although the noise appears greater in

F

~6
(_ o. e~ u

o

g~
,--4

E

O

O

I

!

O

s

0

I

1

I

0

I00

200

300

400

Day of t 9 8 6

Fig. 2. The tabulated Almanac right ascension minus the Almanac's approximate right ascension as a function of the day of 1986.

Algorithm for approximate solar position
kt3 O O

231

k-

'i
I

o

I

I

I

~0

lOO

200 Day of 1986

300

400

Fig. 3. The tabulated Almanac Greenwich mean sidereal time minus The Almanac's approximate Greenwich mean sidereal time as a function of the day of 1986. this plot, the maximum error is about I/lOth that in the previous plot and, therefore, well within the targeted angular accuracy, 3.2 True Almanac values vs some common approximations In the rest of this section, we compare declination only. This is an arbitrary choice, and any of the other variables could be used to economize on the discussion. One of the more c o m m o n ways to represent the solar position is through the use of trigonometric seLrl

ries. Spencer[17] represented the solar declination, sun-earth distance, and equation of time with truncated Fourier series, lqbal[18] suggests these formulas for those interested in a reasonable yet simple approximation. Figure 4 is a plot of the true declination minus the Spencer approximation for each day of 1986. The Almanac approximation differenced with the true declination is also plotted on this scale. The Almanac approximate formula for declination quite clearly is a more accurate representation of the declination.

m o
C 0 .r-+
E x o F_ , 0

~

. . . .. . . . .

Almanac
F'

~S

Lr3 O

I
U Q

0
I

?-

6

I

I

I

= O

t00

200 Days of 1986

300

400

Fig. 4. The Spencer formula for declination and The Almanac formula for declination compared with the tabulated declination for 1986.

232

J.J. MICHALSKY

~ ~

Spencer Cooper

~o

,
S
N .,~

!
,
W-I 0
I I I I

/l Az
V ..... ~t-. i
I I I I

' 0

200

400

600
oays of

800

iooo

I200

I400

1500

19B5-1988

Fig.5. The Cooper formulafor declination and the Spencer formulafordeclination compared to The Al,nanac approximationfor the years 1985-1988. The vertical lines are drawn through the last day of the .',ear.
Perhaps, even more frequently the approximation of declination is made using the formula of Cooper[19]. This formula is a sine wave approximation of the declination. In Fig. 5 The Almanac approximation for declination is taken as the standard against which the Spencer[17] and Cooper[19] formulas are compared? Four consecutive years are compared to show how the uncertainties vary over the leap period. The differences Will repeat in about the same pattern every four years with only a slight shift. Clearly, the Cooper formula is the least accurote of the three we have discussed. On the other hand, it is quite adequate if we only need to position to within about 2 ° . The Spencer formula can be quite accurate if it is used to approximate a particular year's declination, which is the way it was derived originally (note the third year). If it is based on four years, the maximum inaccuracy can be reduced as shown in Fig. 6, but there is little to be gained with this slight change. Therefore, the Fourier series coefficients given by Spencer are not amended here. 3.3 Walraven vs Almanac approximation A question raised in the Zimmerman review[14] was the accuracy of the Walraven[2] formulas. In Fig. 7, the differences of two approximations from the

m
~ m
C 0 .r4 (0 E

"~

...........

o

Spencer" Spencer Coooer

1
I /

, o u
a,1 0

,

,"

/ .. V'L~
~,

" ' " ".:m, ?" '..

?

1

;--

i1
i

eIu E ?
o

i

i

i

i

f

i

0

200

400

600

80o
of

leO0

1200

14o0

1600

Oays

1985-1988

Fig. 6. The same as Fig. 5 except the Spencer-type approximation is rederived based on four years of declination values?

Algorithm

for approximate

solar position

233

7

0

100

200

300

400

Day
Fig. 7. The tabulated Almanac declination
minus

of

1986
formula and the Walraven formula for days in 1986.

The Almanac

II

I

I

I

0

100

200

300 2021 approximation declination in

400

Days
The difference in the

of

Year

Almanac

and Walraven

202 I.

tabulated declination for 1986 are plotted. As we can see, The Almanac approximation performs slightly better than the Walraven approximation for that year. However, the Walraven approximation is well within the accuracy goal of 0.01”. In Fig. 8, we have calculated the difference in the declination between The Almanac and Walraven approximations for the year 202 1. The range of differences is only slightly worse than that indicated in Fig. 7. 4. CONCLUSIONS This article is intended to promote the use of a very accurate algorithm for the calculation of solar

position. The article also compares commonly used approximations. The comparisons made in the last section have not been exhaustive. We have chosen to compare the declination in most instances. However, this should be representative of the other variables that are calculated with The Almanac algorithm. The examples were selected arbitrarily and not with any intent to amplify the errors in competing algorithms. Cooper’s formula for declination is frequently used because it was included in Duffie and Beckman’s[ZO] book Solar Energy Thermal Processes. It is clearly the least accurate of the approximations discussed, but is adequate for many applications. The accuracy

234

J . J . MICHALSKY 2. R. Walraven. Calculating the position of the sun. Solar Energy 20, 393-397 (1978). 3. The American Ephemeris and Nautical Almanac. 1980 Edition. U.S. Gov't. Printing Office, Washington, DC (1979). 4. R. Walraven, Erratum. Solar Energy 22, 195 (1979). 5. C. B. Archer, Comments on "Calculating the position of the sun." Solar Energy 25, 91 (1980L 6. B.J. Wilkinson. An improved FORTRAN program for the rapid calculation of the solar position. Solar Energy 27, 67-68 (1981). 7. L. R. Muir, Comments on "The effect of atmospheric refraction in the solar azimuth. ~ Solar Energy 30, 295 (1983). 8. B. J. Wilkinson. The effect of atmospheric refraction on the solar azimuth. Solar Energy 30, 295 (1983). 9. M. llyas, Solar position programs: Refraction. sidereal time and sunset/sunrise. Solar Energy 31, 437-438 (1983). 10. D. J. B. Pascoe, Comments on "Solar position programs: Refraction. sidereal time and sunrise/sunset. ~ Solar Energy 34, 205-206 (1984). I I. M. Ilyas, Reply to comments by Dr. D. J. B. Pascoe. Solar Energy 34, 206 (1984). 12. B. J. Wilkinson. Re: "Solar position programs: refraction, sidereal time and sunset/sunrise.- Solar Energy 33, 383 (1984). 13. M. llyas, Reply to B. J. Wilkinson's letter on ~Solar position programs: refraction, sidereal time and sunset/sunrise." Solar Energy 33, 383 ~1984). 14. J. C. Zimmerman. Sun-pointing programs and their accuracy. SAND81-0761. Sandia National Laboratories, Albuquerque. NM (1981). 15. TheAstronomicalAlmanac, 1986 Edition. U.S. Gov't. Printing Office. Washington. DC (1985~. 16. U.S. Standard Atmosphere. U.S. Gov't. Printing Office. Washington. DC (1962). 17. J. w . Spencer. Fourier series representation of the position of the sun. Search 2, 172 (19711. 18. M. lqbal, An Introduction to Solar Radiation. Academic Press, Toronto (1983). 19. P. I. Cooper, The absorption of solar radiation in solar stills. Solar Energy 12, 333-346 (19691. 20. J. A. Duffle and W. A. Beckman, Solar Energy Therreal Processes. John Wiley & Sons, New York (1974).

o f S p e n c e r ' s formulas vary with year in the leap cycle. They were d e r i v e d using only a single y e a r ' s declination as input for the least squares fit. T h e y can be i m p r o v e d in the sense that the m a x i m u m deviation is r e d u c e d by fitting to four years o f data as d e m o n strated in Fig. 6. A n o t h e r possibility is to base the coefficients on a single y e a r s ' data, and then add 1/4 day to the angular a r g u m e n t in each o f the next three s u c c e e d i n g years and then repeat the pattern. This is s o m e w h a t a w k w a r d , but k e e p s the f o r m u l a e simple while r e d u c i n g the error, The c o m p a r i s o n s given in Figs. 7 and 8 w o u l d indicate very little reason to c h a n g e to the present algorithm if one has already i m p l e m e n t e d the Walraven code. T h e c o m p a r i s o n in Fig. 8 is the least satisfactory one in the article in that it is not m a d e to true positions. T h e s e were not available for 2021 since we do not have The A l m a n a c ' s high precision routines and The A l m a n a c is p u b l i s h e d for the f o l l o w i n g year only. A s s u m i n g that The A l m a n a c ' s claims for the low precision f o r m u l a s are correct, h o w e v e r , the W a l r a v e n algorithm appears g o o d for several years to c o m e . In c o n c l u s i o n , The A l m a n a c a l g o r i t h m is the best o f the simple a p p r o x i m a t e f o r m u l a s for calculating solar position. It is superior to the calculations based on least squares fits to a g i v e n set o f data since it is based on the l o n g - t e r m p r o g r e s s i o n o f the sun in the ecliptic plane. It is s o m e w h a t s i m p l e r than the alg o r i t h m o f W a l r a v e n , yet is accurate to 0.01 ° tmtil the year 2050. REFERENCES I. L. L. Vant-Hull and A. F. Hildebrandt, Solar thermal power system based on optical transmission. Solar Energy 18, 31-39 (1976).

APPENDIX c c c c c c c c c c c c c c c c c c c c c c c c subroutine sunae (year, day, hour, lat, long, az, el) This subroutine calculates the local azimuth and elevation of the sun at a specific location and time using an approximation to equations used to generate tables in The Astronomical Almanac. Refraction correction is added so sun position is the apparent one. The Astronomical Almanac, U.S. Gov't Printing Ofrice, Washington, DC (1985) input parameters year = year (e.g., 1986) day = day of year (e.g., Feb. 1 = 32) hour = hours plus fraction in UT (e.g., 8:30 A.M. eastern daylight time is equal to 8.5 + 5 (5 hours west of Greenwich) - 1 (for daylight savings time correction) lat = latitude in degrees (north is positive) long = longitude in degrees (east is positive) output parameters a = sum azimuth angle (measured east from north, 0 to 360 °) c c c c e = sun elevation angle (degs) plus others, but note the units indicated before return subroutine sunae (year, day, hour, lat. long, az, el) work with real variables and define some constants, ineluding one to change between degs and radians implicit real (a - z) data twopi, pi, tad/6.2831853. 3.1415927, .017453293/ get the current julian date (actually add 2,400,000 for jd) delta = year - 1949. leap = aint(delta/4.) jd = 32916.5 + delta*365. + leap - day + hour/ 24. 1st no. is mid. 0 jan 1949minus 2.4e6: leap = leap days since 1949 calculate ecliptic coordinates time = jd - 51545.0 51545.0 + 2.4e6 = noon I jan 2000 force mean longitude between 0 and 360 degs

c c

c c c c c c c

A l g o r i t h m for a p p r o x i m a t e s o l a r p o s i t i o n mnlong = 280.460 + .9856474*time mnlong = mod(mnlong, 360.) if(mnlong.lt.0.)mnlong = mnlong + 360. c m e a n a n o m a l y in r a d i a n s b e t w e e n 0 a n d 2 * p i mnanom = 357.528 + .9856003 * time mnanom = mod(mnanom, 360.) if(mnanom.lt.0.)mnanom = mnanom + 360. mnanom = mnanom*rad c o m p u t e ecliptic l o n g i t u d e a n d o b l i q u i t y o f ecliptic in radians eclong = mnlong + 1.915*(mnanom) + 0.20* sin(2.*mnanom) 1 eclong = mod(eclong.360.) i f ( e c l o n g . l t . O . ) e c l o n g = e c l o n g + 360. oblqec = 23.429 - .O000004*time eclong = eclong*rad oblqec = oblqec*rad c a l c u l a t e right a s c e n s i o n a n d d e c l i n a t i o n num= cos(oblqec)*sin(eclong) den = c o s ( e c l o n g ) ra = a t a n ( n u m / d e n ) force ra b e t w e e n 0 a n d 2*pi if(den.lt.0)then ra = ra + pi elseif(num.lt.0)then ra = ra + t w o p i endif dec in r a d i a n s dec = a s i n ( s i n ( o b l q e c ) * s i n ( e c l o n g ) ) c a l c u l a t e G r e e n w i c h m e a n sidereal t i m e in h o u r s gmst = 6.697375 + .0657098242*time + hour h o u r not c h a n g e d to s i d e r e a l s i n c e ' t i m e ' i n c l u d e s the fractional d a y gmst = mod(gmst, 24.) i f ( g m s t . l t . 0 . ) g m s t = g m s t + 24. calculate local m e a n sidereal t i m e in r a d i a n s lmst = g m s t + l o n g / 1 5 , imst = m o d ( l m s t . 2 4 . ) i f ( l m s t . l t . 0 . ) l m s t = lmst + 24. Imst = I m s t * 1 5 . * r a d c c

235

calculate h o u r a n g l e in r a d i a n s b e t w e e n - p i and pi h a = lmst - ra if(ha.It. - pi)ha = h a + t w o p i if(ha.gt.pi)ha = ha - twopi c h a n g e latitude to r a d i a n s lat = lat*rad calculate azimuth and elevation el = a s i n ( s i n ( d e c ) * s i n ( l a t ) + cos(dec)*cos(latl*cos(ha)) az = a s i n ( - c o s ( d e c ) * s i n ( h a ) / c o s ( e l ) ) if az = 9 0 d e g s . elcritieal = a s i n ( s i n ( d e c ) / s i n ( l a t ) ) elc = a s i n ( s i n ( d e c ) / s i n ( l a t ) ) i f ( e l . g e . e l c ) a z = pi - az i f ( e l . l e . e l c . a n d . h a . g t . O . ) a z = t w o p i + az this puts a z i m u t h b e t w e e n 0 a n d 2*pi r a d i a n s c a l c u l a t e r e f r a c t i o n c o r r e c t i o n f o r US stan. a t m o s p h e r e need to h a v e el in d e g s b e f o r e c a l c u l a t i n g c o r r e c t i o n el = e l / r a d if(el.gt. - . 5 6 ) t h e n refrac = 3 . 5 1 5 6 1 " ( . 1 5 9 4 + . 0 1 9 6 " e l + .00002"e1"'2)/ 1 (1. + . 5 0 5 " e l + . 0 8 4 5 " e 1 " ' 2 ) else refrac = .56 endif note that 3 . 5 1 5 6 1 = 1 0 1 3 . 2 m b / 2 8 8 . 2 C el = el + refrac elevation in d e g s c o n v e r t az a n d lat to d e g s before r e t u r n i n g az = a z / r a d lat = l a t / r a d m n l o n g in d e g s , g m s t in hours, j d in d a y s if 2 . 4 e 6 added; m n a n o m , e c l o n g , o b l q e c , ra, d e c , Imst. a n d h a in radians return end

c c c c

c c c

c c

c c

c c c c

c

c c c c c c

c c c c

c c

c c


相关文章:
更多相关标签: