@@ -22,9 +22,6 @@ import (
2222 "time"
2323)
2424
25- var t , ra , dec , ym , y0 , yp , xe , ye , z1 , z2 float64
26- var nz int
27-
2825/*
2926RiseSet holds the rise and set times as strings in the form hh:mm
3027*/
@@ -57,6 +54,10 @@ Example
5754
5855*/
5956func Riseset (object Object , eventdate time.Time , glong float64 , glat float64 , zone float64 ) (results RiseSet ) {
57+
58+ var ra , dec , ym , y0 , yp , xe , ye , z1 , z2 float64
59+ var nz int
60+
6061 sinho := make ([]float64 , 4 )
6162
6263 day := eventdate .Day ()
@@ -85,18 +86,23 @@ func Riseset(object Object, eventdate time.Time, glong float64, glat float64, zo
8586 rise := 0
8687 sett := 0
8788 hour := 1.
88- ym = sinalt (iobj , date , hour - 1 , glong , cl , sl ) - sinho [iobj ]
89+ ym ,ra ,dec = sinalt (iobj , date , hour - 1 , glong , cl , sl , ra , dec )
90+ ym = ym - sinho [iobj ]
8991
9092 for (hour != 25 ) && (rise * sett != 1 ) {
91- y0 = sinalt (iobj , date , hour , glong , cl , sl ) - sinho [iobj ]
92- yp = sinalt (iobj , date , hour + 1 , glong , cl , sl ) - sinho [iobj ]
93+ y0 ,ra ,dec = sinalt (iobj , date , hour , glong , cl , sl , ra , dec )
94+ y0 = y0 - sinho [iobj ]
95+
96+ yp ,ra ,dec = sinalt (iobj , date , hour + 1 , glong , cl , sl ,ra ,dec )
97+ yp = yp - sinho [iobj ]
98+
9399 xe = 0
94100 ye = 0
95101 z1 = 0
96102 z2 = 0
97103 nz = 0
98104
99- quad () // Note uses and updates package variables
105+ nz , yp , y0 , ym , xe , ye , z1 , z2 = quad (yp , y0 , ym , xe , ye , z1 , z2 )
100106
101107 switch nz {
102108 //cases depend on values of discriminant
@@ -220,9 +226,9 @@ Returns the local siderial time for the mjd and longitude specified
220226func lmst (mjd float64 , glong float64 ) float64 {
221227 mjd0 := ipart (mjd )
222228 ut := (mjd - mjd0 ) * 24
223- t := (mjd0 - 51544.5 ) / 36525
229+ mytim := (mjd0 - 51544.5 ) / 36525
224230 gmst := 6.697374558 + 1.0027379093 * ut
225- gmst = gmst + (8640184.812866 + (.093104 - .0000062 * t ) * t ) * t / 3600
231+ gmst = gmst + (8640184.812866 + (.093104 - .0000062 * mytim ) * mytim ) * mytim / 3600
226232 return 24 * fpart ((gmst - glong / 15 )/ 24 )
227233}
228234
@@ -242,9 +248,9 @@ Finds a parabola through three points and returns values of coordinates of
242248extreme value (xe, ye) and zeros if any (z1, z2)
243249Assumes that the x values are -1, 0, +1
244250*/
245- func quad () {
251+ func quad (yp , y0 , ym , xe , ye , z1 , z2 float64 )( int , float64 , float64 , float64 , float64 , float64 , float64 , float64 ) {
246252 var a , b , c float64
247- nz = 0
253+ nz : = 0
248254 a = 0.5 * (ym + yp ) - y0
249255 b = 0.5 * (yp - ym )
250256 c = y0
@@ -265,7 +271,7 @@ func quad() {
265271 z1 = z2
266272 }
267273 }
268- return
274+ return nz , yp , y0 , ym , xe , ye , z1 , z2
269275}
270276
271277//Returns SIN of x degrees
@@ -283,29 +289,29 @@ Returns sine of the altitude of either the sun or the moon given the modified
283289julian day number at midnight UT and the hour of the UT day, the longitude of
284290the observer, and the sine and cosine of the latitude of the observer.
285291*/
286- func sinalt (iobj Object , mjd0 float64 , hour float64 , glong float64 , cphi float64 , sphi float64 ) float64 {
292+ func sinalt (iobj Object , mjd0 float64 , hour float64 , glong float64 , cphi float64 , sphi float64 , ra float64 , dec float64 ) ( float64 , float64 , float64 ) {
287293 instant := mjd0 + hour / 24.
288- t = (instant - 51544.5 ) / 36525
294+ tim : = (instant - 51544.5 ) / 36525
289295 if iobj == 1 {
290- moonsub ()
296+ ra , dec = moonsub (tim , dec , ra )
291297 } else {
292- sun ()
298+ ra , dec = sun (tim , dec , ra )
293299 }
294300 tau := 15 * (lmst (instant , glong ) - ra ) //hour angle of object
295- return sphi * sn (dec ) + cphi * cn (dec )* cn (tau )
301+ return sphi * sn (dec ) + cphi * cn (dec )* cn (tau ), ra , dec
296302}
297303
298304/*
299305Returns RA and DEC of Sun to roughly 1 arcmin for few hundred years either side
300306of J2000.0
301307*/
302- func sun () {
308+ func sun (tim , dec , ra float64 ) ( float64 , float64 ) {
303309 p2 := 6.283185307
304310 COSEPS := .91748
305311 SINEPS := .39778
306- m := p2 * fpart (.993133 + 99.997361 * t ) //Mean anomaly
312+ m := p2 * fpart (.993133 + 99.997361 * tim ) //Mean anomaly
307313 dL := 6893 * math .Sin (m ) + 72 * math .Sin (2 * m ) //Eq centre
308- L := p2 * fpart (.7859453 + m / p2 + (6191.2 * t + dL )/ 1296000 )
314+ L := p2 * fpart (.7859453 + m / p2 + (6191.2 * tim + dL )/ 1296000 )
309315 // convert to RA and DEC - ecliptic latitude of Sun taken zero
310316 sl := math .Sin (L )
311317 x := math .Cos (L )
@@ -317,7 +323,7 @@ func sun() {
317323 if ra < 0 {
318324 ra = ra + 24
319325 }
320- return
326+ return ra , dec
321327}
322328
323329/*
@@ -326,16 +332,18 @@ centuries either side of J2000.0
326332Predicts rise and set times to within minutes for about 500 years in past
327333TDT and UT time diference may become significant for long times
328334*/
329- func moonsub () {
335+
336+
337+ func moonsub ( tim ,dec , ra float64 ) (float64 ,float64 ){
330338 p2 := 6.283185307
331339 ARC := 206264.8062
332340 COSEPS := .91748
333341 SINEPS := .39778
334- L0 := fpart (.606433 + 1336.855225 * t ) //mean long Moon in revs
335- L := p2 * fpart (.374897 + 1325.55241 * t ) //mean anomaly of Moon
336- LS := p2 * fpart (.993133 + 99.997361 * t ) //mean anomaly of Sun
337- d := p2 * fpart (.827361 + 1236.853086 * t ) //diff longitude sun and moon
338- F := p2 * fpart (.259086 + 1342.227825 * t ) //mean arg latitude
342+ L0 := fpart (.606433 + 1336.855225 * tim ) //mean long Moon in revs
343+ L := p2 * fpart (.374897 + 1325.55241 * tim ) //mean anomaly of Moon
344+ LS := p2 * fpart (.993133 + 99.997361 * tim ) //mean anomaly of Sun
345+ d := p2 * fpart (.827361 + 1236.853086 * tim ) //diff longitude sun and moon
346+ F := p2 * fpart (.259086 + 1342.227825 * tim ) //mean arg latitude
339347 // longitude correction terms
340348 dL := 22640 * math .Sin (L ) - 4586 * math .Sin (L - 2 * d )
341349 dL = dL + 2370 * math .Sin (2 * d ) + 769 * math .Sin (2 * L )
@@ -365,5 +373,5 @@ func moonsub() {
365373 if ra < 0 {
366374 ra = ra + 24
367375 }
368- return
376+ return ra , dec
369377}
0 commit comments