FEB 2 3 1965

^ecknlcui ^iote ^0.303

EXTENSION OF PROGRAMS FOR CALCULATIONS

OF GREAT CIRCLE PATHS AND

SUNRISE-SUNSET TIMES

J. H. Crary

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
THE NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards is a principal focal point in the Federal Government for assuring

maximum application of the physical and engineering sciences to the advancement of technology in industry and commerce. Its responsibilities include development and maintenance of the national stand- ards of measurement, and the provisions of means for making measurements consistent with those standards; determination of physical constants and properties of materials; development of methods

for testing materials, mechanisms, and structures, and making such tests as may be necessary, particularly for government agencies; cooperation in the establishment of standard practices for incorporation in codes and specifications; advisory service to government agencies on scientific and technical problems; invention and development of devices to serve special needs of the Government; assistance

to industry, business, and consumers in the development and acceptance of commercial standards and

simplified trade practice recommendations; administration of programs in cooperation with United

States business groups and standards organizations for the development of international standards of practice; and maintenance of a clearinghouse for the collection and dissemination of scientific, technical, and engineering information. The scope of the Bureau's activities is suggested in the following

listing of its four Institutes and their organizational units. Institute for Basic Standards. Electricity. Metrology. Heat. Radiation Physics. Mechanics. Applied Mathematics. Atomic Physics. Physical Chemistry. Laboratory Astrophysics.* Radio Standards Laboratory: Radio Standards Physics; Radio Standards Engineering.** Office of Standard Ref- erence Data. Institute for Materials Research. Analytical Chemistry. Polymers. Metallurgy. Inorganic Materials. Reactor Radiations. Cryogenics.** Office of Standard Reference Materials.

Central Radio Propagation Laboratory.** Ionosphere Research and Propagation. Troposphere

and Space Telecommunications. Radio Systems. Upper Atmosphere and Space Physics.

Institute for Applied Technology. Textiles and Apparel Technology Center. Building Research.

Industrial Equipment. Information Technology. Performance Test Development. Instrumentation.

Transport Systems. Office of Technical Services. Office of Weights and Measures. Office of Engineering Standards. Office of Industrial Services.

* NBS Group, Joint Institute for Laboratory Astrophysics at the University of Colorado. ** Located at Boulder, Colorado.
NATIONAL BUREAU OF STANDARDS

^ecltnical c

T[ote> 303

Issued February 18, 1965

EXTENSION OF PROGRAMS FOR CALCOLATIONS

OF GREAT CIRCLE PATHS AND

SUNRISE-SUNSET TIMES

J. H. Crary

Central Radio Propagation Laboratory

National Bureau of Standards

Boulder, Colorado

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U.S. Government Printing Office Washington D.C. 20402

Price: 25 cents

CONTENTS

Page

1. Introduction 1

2. The Computer Programs 2

2.1 Introduction 2

2.2 Sunrise-Sunset Program (DALITP) 3

2.3 Sunrise-Sunset Program with Bearing Calculations

(DALITB) . . h

2.4 Great Circle Path Subroutine (PATH) h

2.5 Solar Angle Subroutine (SOLANG) 9

3. Some Illustrations of the Results for Long and Short Paths 11

h. Acknowledgments 13

5. References 13

Figures lk

in

EXTENSION OF PROGRAMS FOR CALCULATIONS OF GREAT

CIRCLE PATHS AND SUNRISE-SUNSET TIMES

J, H. Crary

A previous note by Brady and Crombie [1964]

describes methods and a program for calculating and

plotting sunrise and sunset times at specified

heights and distances along a great circle path

determined by the locations of the end points . The

program also plots the distribution of darkness or daylight along the path versus time.

This program has been modified to remove some restrictions on the use of the earlier program, and

to provide the option of calculating the short or long great circle path. In addition, transmitter,

receiver, and solar bearings are also calculated and

plotted. Some interesting characteristics of the

short and long paths from Rugby, England, to Byrd

Station, Antarctica, are shown.

1. Introduction

A previous note by Brady and Crombie [1964] described the methods

and gave the formulas for the calculation of sunrise and sunset times

at specified heights and distances along a great circle path. The

program previously described had some limitations which restricted its

use. Modifications were made to the program to remove the preceding

restrictions and to add some additional features.
These features are the calculation of (l) either the short or long

path, (2) the transmitter bearing at each segment of the path, (3) the

solar bearing angle at sunrise and sunset at each segment, (k) the difference angle between the transmitter and solar bearings, and (5) the cosine of this difference angle. A routine was then added to plot these

results, together with the percentage of darkness along the path. The

necessity of giving the easternmost end of the path first was also removed, since the choice depends upon whether the short or long path is

desired.

These modifications were accomplished by adding three subroutines to

the original program. These programs and subroutines are available in

Fortran language. A description of these will be given which includes

the necessary equations and logical statements to enable a person to calculate these results manually or by means of a computer program, without

first determining how to make the necessary decisions; this is a time consuming process and need not be duplicated.

2 . The Computer Programs

2.1. Introduction

The computer programs consist of two main Fortran programs, DALITP

and DALITB, together with the necessary subroutines (PATH, SOLANG, DATBX,

PLOT). The functions of these subroutines are: (l) PATH--to calculate

the necessary information about the great circle path, such as the total

length, the bearing angles of the transmitter and receiver, and the position of points at specified intervals along the path; (2) SOLANG—to
calculate the solar zenith and bearing angles; (3) DATBX—a function

which furnishes the date on which the calculations are performed to the

calling program; and (k) PLOT--a general output plotting routine. The

functions of these subroutines are discussed in 2.2. - 2.5. The following notation will be used in the description:

L = latitude in degrees (North is +, South is -)

= longitude in degrees (East is +, West is -)

LHA. = solar local hour angle in degrees

DECL = solar declination in degrees

Subscript i or j indicates ith or jth point of the path

R indicates receiver

T indicates transmitter

2.2. Sunrise-Sunset Program (DALITP)

This program is the modified version of the original DALIT program

and controls the data card reading and the output printing. It also performs the calculations of the sunrise and sunset times, using the same

equations given by Brady and Crombie [196U].

In the case of | cos LHA| > 1, the following test is made:

if (L. • DECL) < the path is always dark,

if (L

i

• DECL) 2: the path is always light.

This takes care of both northern and southern hemisphere conditions.

The data inputs are the same as before except for the following
changes: (l) a given number of pairs of ionosphere and screen heights,

up to a maximum of ten, may be used, and (2) the option of calculating

the long or short path is provided.

2.3. Sunrise -Sunset Program with Bearing Calculations (DALITB)

This program is the same as DALITP, with the addition of the subroutine SOLANG to calculate the solar bearing angles between the solar

and transmitter bearing, and the cosine of this difference. The PLOT

subroutine then is used to plot the solar bearings, the cosine of the

difference angle, and the distance along the path where the terminator*

crosses, as a function of time.

2.1+. Great Circle Path Subroutine (PATH)

This subroutine is called by the main program to do the great circle

calculations. The program has several options to calculate for the short

or long path: (l) distance between two specific points; (2) latitude,

longitude, and bearing of transmitter at given intervals on a given path;

(3) coordinates of reflection points for an n-reflection ray between given

points; or (k) coordinates of the end points of a path when given one end

and the jth point on an n-reflection path.

The option of being able to calculate either the short or the long

great circle path between two specified points requires several logical

* Terminator--a useful astronomical term for the line separating the dark

and sunlit portions of a planet, applied here at the ionospheric height

specified.
decisions to be made at different points during the calculation in order

to insure that the proper solution is being obtained. The spherical

triangles involved in these solutions are illustrated in figure 1. The

polar angle of the spherical triangle is the difference in longitude between the two points on the path being considered. The algebraic sign of

the difference is not used in the calculations. Four cases can actually

be defined, depending upon whether the short or long path is desired, and

whether the transmitter or starting end of the path is the eastern-most or

western-most point of the path. This decision depends upon whether the

short or long path is desired.

In addition, decisions must be made as to which of the possible

solutions to equations involving cos"1 will yield the desired path. It

is also necessary to account for some of the cases of the degenerate

spherical triangle, such as where both points of the path lie on a

meridian, etc. For this reason, the necessary equations and logic statements are included, permitting a straightforward calculation of the

desired path. The flagl mentioned in the following description is a

logical variable which indicates a decision made on the basis of the

orientation of transmitter and receiver and whether the short or long

path is desired.

(l) Calculate the polar angle or longitude difference, |^ - $~\

(see figure l)

(2.1) If |fzL -

<f>

| < 180°, continue; if > l80° , go to (2.2).

(2.1.1) If short path, <f> = 360 -

\<f> - </> | and set flagl = 2, go to

(3).
(2.1.2) If long path, = |0T

- R | , and set flagl = 1, go to (3).

(2.2) If |0T

- R

| > 180°.

(2.2.1) If short path, = |0T

- R | , and set flagl = 1, go to (3).

(2.2.2) If long path, = 36O - |0T

- R | , and set flagl = 2.

(3) Calculate D, the total path distance.

(3.1) If cos LR

± 0, go to (3.2); if = 0, continue. If cos % = 0,

then DRAD = |LT

- LJ * tt/180 radians.

(3.1.1) If LR

= 0, impossible branch, error message, exit.

(3.1.2) If LR < 0, T = tt rad, where T is the angle of the spherical

path triangle at the transmitter.

(3.1.3) If LR > 0, T = 0.

(3.1.4) If long path DRAD - DRAD - 2tt. D = DRAD * R where R is the

effective earth radius. R is taken as the temperate latitude

average of the radii determined from the polar and equatorial

circumferences, or (3 Rea + R^ol)/^*

r^ae calculation then

proceeds to (k) .

(3.2) If cos L

R

^ 0, DRAD = cos"1 [sin Lm sin L-^ + cos Lrp cos LR

cos 0] from the cosine law for the spherical triangle shown

in figure 1.

(3.2.1) If sin (DRAD) = 0, then end points are antipodal, and there

are an infinite number of solutions. This causes a message

to be printed and the calculation is terminated.

(3.2.2) If sin (DRAD) ^ 0: for long path DRAD = 2tt - DRAD; for both

paths D = DRAD • R. If only distance is desired, return to

main program.
(3.2.3) Determine number of segments. Complete calculation of angle

T.

(4.1) If cos L

T

= 0, continue; if # 0, go to (4.3)

(U.l.l) If L = 0, error message, exit.

(4.1.2) If L

T < 0, T = 0, go to (4.1.4).

(4.1.3) If L

T

> 0, T = tt.

(4.1.4) Then Lj = L

T

- (LT/|L

T

|) D-/R, B = 90 - (LT/|LT | ) 90, =

R

; go "to (4.3).

(4.2) If cos L

T f 0.

(4.3) T = cos"1 [(sin L-n - cos D sin L )/(cos L<p sin D)] from the

law of cosines again applied to the path triangle of figure 1.

(5) Calculate the latitudes at the jth segment along the p~i,th:

L . = 90 - cos" 1 (sin Lm cos D^ + cos L sin D. cos T) . This

is another variation of the law of cosines applied to the triangle of the pole, point j and the transmitter. The use of

90 - cos"1 instead of sin"1 removes an ambiguity in the result,

(5.1) If cos L. =0, A0 = 0, since point j is at a pole. Here A0

J is the difference in longitude between the transmitter and

point j. Go to (6.3)

.

(5.2) If cos L. / 0.

J

(6) Calculate the longitudes.

(6.l) Cos A0 = (cos D. - sin L^ sin L.)/(cos LTM cos L.). This also

comes from the law of cosines.
(6.2) Now, using the law of sines for sin A0: sin A0 = (sin T sin

Dp/Ccos Lj).

(6.2.1) If sin A0 ^ 0, A0 = cos"1 (cos A0), go to (6.3).

(6.2.2) If sin A0 < 0, A0 = 36O - cos"1 (cos A0).

(6.3) If (LT

- L ) £ 0, continue; if < 0, go to (6.4).

(6.3.1) If flagl = 1 [see (2.l)-(2 .2)] , 0j = T

- A0, go to (7).

(6.3.2) If flagl = 2, = T

+ A0, go to (7).

(6.1+) If (LT

- LR

) < and

(6.4.1) If flagl = 1 use (6.3.2), go to (7).

(6.4.2) If flagl = 2 use (6.3.1).

(7)

(7.1)

Check longitude for proper hemisphere

If |0.| - 180 £ 0, continue; if > 0, go to (8).

J (7.1.1) If 0- < 0, 0. = 0. + 360°, gO to (8).

(7.2) If |0,j ^ 0, = 360° + 0..

(8).

(8.1)

(8.2)

(8.3)

Find transmitter -bearing B..

If cos \ = t, B = 90 + (I|

R)/C1L

R

!0

• 90, go to (9).

Or if cos L_ ^ 0, B. = cos

-1 [(sin Lm - sin L. cos D.)/

(cos L. sin D.)] again from the law of cosines.

3 3

If (^ -

<f> ) < 0, continue; if > 0, go to (8.4).

(8.3.1) Flagl = 1 [see (2.l)-(2.2)] B. = 36O - B., go to (9)

J J (8.3.2) Or if flagl =2, B. = B., go to (9).

J J (8.4) if (^ - ^) * 0.

(8.4.1) And flagl = 1, use (8.3.2), go to (9).

(8.4.2) Or if flag 1=2, use (8.3.1).

8
(9) Find receiver "bearing from transmitter.

(9.1) If (0T

- R

) < 0, continue; if S 0, go to (9.2).

(9.1.1) And flagl = 1, BR

= T, go to (10).

(9.1.2) Or if flagl =2, B

R

= 36O - T, go to (10).

(9.2) If (0T

- R

) £ 0.

(9.2.1) And flagl = 1 use (9.1.2), go to (10).

(9.2.2) Or if flagl = 2 use (9.I.1).

(10) Since the transmitter "bearing is undefined at the transmitter

we use:

(10.1) Bj. = Bg + 180.

(10.2) If (B1 - 360) £ 0, B-l = B1 - 360, go to (ll),

(10.3) Or if (Bx - 360) < 0, B1 = Bj_.

(11) The subroutine returns to the main program here.

2.5. Solar Angle Subroutine (SOLANG)

This subroutine finds the solar bearing angles S. and the zenith

angles \±

.

The inputs are an array of N sets of path latitude (XLAT) and longitude (XLON) , solar declination (DECL) , and Greenwich hour angle (GHA),

where N ^ 1000. The variable NTIMES may be used to specify whether there

is only one pair of values (NTIMES = l) for GHA and DECL, or whether there

is a pair for each point (NTIMES > l).

The logic statements given here are used to find the proper solution

to equations involving the arc cosine function. It is necessary to pick
the solution corresponding to the desired triangle. It is also necessary

to make decisions so that the solar bearing angle will follow the convention of being measured clockwise from North and the solar zenith angle

will lie between -l80° and +l80°

The calculation then proceeds as follows, using the solar triangle

of figure 2

(1) Find the solar local hour angle SLHA = GHA + $..

(1.1) If SEHA < 0, then SLHA = SLHA + 36O, go to (l.3).

(1.2) If SLHA > then

(1.3) If (SLHA - 360) > 0, continue; if < 0, go to (1.4).

(1.3.1) SLHA = SLHA - 36O, go to (2).

(1.4) If (SLHA - 360) < 0, SLHA = SLHA.

(2) Calculate v. = cos

-1 (sin L. sin DECL + cos L. cos DECL

M. 1 1 cos SLHA) from the law of cosines. Next find the bearing

angle.

(2.1) If I L.| - 90 = 0, continue; if ^ 0, gc to (2.2)

(2.1.1) S. = SLHA, go to (2.5).

(2.2) If |L.| - 90 ^ 0.

(2.3) If (xlxl " l8°) = °> continue; if ^ 0, go to (2.4).

(2.3.1) And if L. - DECL < 0, S. = l80° , go to (2.5).

(2.3.2) Or if L. - DECL > 0, S. = 0°, go to (2.5).

(2.4) If x (| x | - 180) + 0.

(2.4.1) S. = cos

-1 [(sin DECL - sin L. cos x)/(cos L. sin y)] from the

law of cosines.

(2.5) If (SLHA - 180) < 0, S. = 360 - S., go to (3).

10
(2.6) Or if (SLHA - 180) ^ 0, S^^ = S±

.

(3) Return to the main program.

3. Some Illustrations of the Results for Long and Short Paths

Recordings of VLF signals propagating over long VLF paths (> 15,000

km or so) often show diurnal variations in phase and amplitude which are

indicative of at least some energy coming over the long great circle path

[Crombie, 1958] , particularly when the short path is completely or almost

completely sunlit. This makes the variation of daylight or darkness over

both the long and short paths especially interesting.

Figures 3 and k show typical sets of the computer output for sunrise

and sunset calculations for the short and long path from GBR, Rugby, England, to Byrd Station, Antarctica, for April 15, 1962 . The upper set in

each figure is taken for an ionospheric height of 90 km, and a screen

height of 27 km, or x ~ 98^ • The lower set in each figure is for ground

level or x

= 90 . These figures illustrate several interesting points.

First is the essentially complementary nature of the two paths, i.e., one

is becoming light when the other is becoming dark, and vice versa. An<r

other is the extremely rapid transition, especially at ground level, around

1900 UT. This could be expected to cause a very rapid transition of the

predominant component of the received signal from the long path to the

short path. This effect occurs when the great circle joining the transmitter and receiver nearly coincides with the terminator. It is indicated

by a solar and transmitter bearing difference close to 90° or 270° at sunrise or sunset.

11
These effects are also illustrated in figures 5 and 6. These show

the output plot of percentage of path in darkness, distance along path of

terminator location, the solar bearing angle, and the cosine of the solartransmitter bearing angle difference at sunrise and sunset. The rapid

transition occurs when the bearing angle at which the path crosses the

equator (9.63° ) is nearly equal in magnitude to the solar declination

(9.6^°). The bearing at which the path crosses the equator is also equal

to the complement of the maximum latitude (80.37°) on the great circle

path. The bearing difference angle is interesting, not only because of

its influence on the time interval during which the night-day transition

is made, but also because experimental data indicate that the angle at

which the terminator crosses the path appears to influence the amount of

mode conversion at the discontinuity in the waveguide height caused by

the terminator [Crombie, I96U]

Figure 6 illustrates an interesting characteristic of some paths at

certain times . Note that the path at ground level remains about 80$ dark

until I85O UT. The percent of darkness then drops to 18$ in about fifteen

minutes. The percentage in darkness at the ionospheric height changes

gradually from about 1^30 UT. It begins to drop rapidly after 1700 and

reaches zero at 1825 UT. It then rises from zero to 10$ between 1950 and

2120 UT.

A corresponding phase retardation at about the proper time (1958 UT

in this case) on a path showing such behavior may enable one to determine

the critical value of x, at which the ionizing radiation first affects

the D layer.

12
k. Acknowledgment s

The author is indebted to D. D. Crombie for his suggestions which

added considerably to the usefulness of these programs. The work was

sponsored by the Advanced Research Projects Agency, Washington, D. C,

under Order No. 183.

5. References

Brady, A. H., and D. D. Crombie (196U), Calculation of sunrise and sunset

times at ionospheric heights along a great circle path, NBS Technical Note No. 209.

Crombie, D. D. (1958), Differences between the EW and WE propagation of

VLF signals over long distances, J. Atmos. Terrest. Phys. 12, 110-

117.

Crombie, D. D. (196U), Periodic fading of VLF signals received over long

paths during sunrise and sunset, Radio Sci. J. Res. NBS 68D, No. 1,

27-3^.

13
Figure

Geometry of a Great Circle Path

1U
-^^IIAL^EQUATOR

Figure 2

The relationship of the solar local hour angle, zenith angle, declination and bearing angle,

15
i noo* ui. ioim

III i

HEi

[TOOife

i i

11 11 ll

I 1 IJ:

l!SS: : :188:

ill*: : :ip:

hi- . ?:: IJ1& E jI:::::::::::::::::::::::;:::::::::::::::::::::::::::::::::::::::: :

j|:

III: I

1115- "III'IIC""- - «bI life , :i::::!::::E E Hi

il»l

: !ii;iiiiiE I

If; Ill*- - - «. Hp =iiiiiiiE :1

:??:

iii- •••- - ill IHt •:::: : «:

IS!?: :::: :»:

Figure 3. Light and Darkness Plot vs. UT from the Short Path from GBR, Rugby, England, to Byrd Station, Antarctica for 15 April (1962), for (a) ground or y = 90 , ang

(b) ionospheric reflection height .ss 90 km, screen height = 27 km, or \ = 9& .

16
iD*T.*n.t, noon ui. ioi»i OEC lONl p»Tm i - CHi

1 OtC, L»T. S2,« D£G »1fl |

:«:

:??.-:•

liii::

ilii::

Figure 4. Light and Darkness Plots vs. UT for Long Path from GBR, Rugby, Engl

Byrd Station, Antarctica for 15 April (1962), for (a) ground or x

(b) ionospheric reflection height - 90 km, screen height = 27 km, o

17
Figure
i !i1

•| -

&&&&&a&&o.a

1

"

-o

3

3
"ill* !

"' i ;: i c| • I; . . -i 1 1 ' 1 i L- ; 1 ill! i ;= 

")-' ! i ! "••-,! ! I.--"1 1 .'..:: I ! i \ j!:fn i i i ! 1 i "

""'! j"-"" fill

n

[l|sfF§i" ' i ;

-. >"-' ! ! ' i ! ' .i

I'll . * : ! I 1 ! ! ; ?.' i \ ;

Si i i ! i i : \f ' ;:,::; :

ii ! i i ! s"1 ' :"*•.! ! i i ! ! ! : 1 I j i i ! ' / ' ';"""!' "! T--...:I...... Ill 1

./'" | 1 :-j--.....| 1 I ! ? 1 ] ! ! ! ill! l . . . 

1H ' i 1 7 -. •*

:;;;[

IP. . ! II!!.!;

.... 1 ! 1 i - s> . . ..! 1 I ,..*" i \i 

;ir 1 | _-~ Mi^H, 1 "•"'''. ._._. 1 /" 1 :=r Hill ="?* j

.1 1 f" 1

:S;s* | ; * 

i ! ; : 1

liii ::•?; 11 Is I \i i i j | | ! 1 i

(ill* 1 1 i

....... i i i !1 1 -|; = = > i j

i':: i': . . . i 1 - i _.-

~ : : • j [ \ !11

;['A i > i i

:' *ti~K

%

£

S

\ j o% "

'& ' ' « i J , H,
I'M I

!i "Jl

4\

jijifi i Ml

mmh

iiil

L | .1

-_^j.. r..,-<" 1 I : : | ! i ; ! 1 i i i J

,» ] • ! • '""- >i j

..• Hi i ii • !, 1

V-I.-I- •• "• | •'";

\

—'- : < • t I .1 I '1 i . . . i ' ! ' l t . . . t

liill

mil

in: » E ! 

-ii I 1 \ ! ! ! 1 i 1 ! ! 3 i i j l

i!i'Mi

ii3>LSiL(i

Ii': - -

lid

r 11 1 5 <rf=i- -'

....

••-...

...-••••" X.\. 1 \ i ' "--, • i 1

J'

""•H.

\ 1 v„ . . 1- •: • :•'" ^ " - i t: 1 •

;nl

V'i

Ji!

if

iM'i

!' • r :; 8 .

'"i

ifI ' ! I I'