Package: astronomical-almanac / 5.6-7

trnsit.c.patch Patch series | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
Description: Fix rise and set time convergence failure at high southern latitudes.
Author: Stephen Moshier
Forwarded: patch is from upstream
--- astronomical-almanac-5.6/trnsit.c	2005-11-08 19:46:37.000000000 -0500
+++ new/trnsit.c	2020-03-29 12:56:44.668203604 -0400
@@ -6,6 +6,7 @@
 
 #include "kep.h"
 extern double glat;
+extern double azimuth, elevation; /* degrees */
 
 /* Earth radii per au */
 #define DISFAC 2.3454780e4
@@ -19,7 +20,16 @@
 double r_trnsit;
 double r_rise;
 double r_set;
+double elevation_threshold, semidiameter;
 int f_trnsit;
+int found_rise, found_set;
+int no_rise_set (double, int (*)());
+
+double search_halve (double tl, double yl, double th, double yh, int (* func)());
+int iter_trnsit( int (* func)() );
+static void iter_func( double t, int (* func)() );
+
+
 
 int trnsit(J, lha, dec)
 double J; /* Julian date */
@@ -29,9 +39,9 @@
 double x, y, z, N, D, TPI;
 double lhay, cosdec, sindec, coslat, sinlat;
 
-f_trnsit = 0;
 TPI = 2.0*PI;
 /* Initialize to no-event flag value. */
+f_trnsit = 0;
 r_rise = -10.0;
 r_set = -10.0;
 /* observer's geodetic latitude, in radians */
@@ -42,12 +52,12 @@
 cosdec = cos(dec);
 sindec = sin(dec);
 
-/*  x = (J - floor(J-0.5) - 0.5) * TPI; */
- x = floor(J) - 0.5;
+/* Refer to same start of date as iter_trnsit,
+   so r_trnsit means the same thing in both programs.  */
+ x = floor(J - 0.5) + 0.5;
  x = (J - x) * TPI;
 /* adjust local hour angle */
 y = lha;
-/* printf ("%.7f,", lha); */
 while( y < -PI )
 	y += TPI;
 while( y > PI )
@@ -55,8 +65,7 @@
 lhay = y;
 y =  y/( -dradt/TPI + 1.00273790934);
 r_trnsit = x - y;
-/* printf ("rt %.7f ", r_trnsit); */
-/* Ordinarily never print here.  */
+/* Ordinarily do not print here.  */
 if( prtflg > 1 )
 	{
 	printf( "approx local meridian transit" );
@@ -72,7 +81,11 @@
 switch( objnum )
 	{
 /* Sun */
-	case 0: N = COSSUN; break;
+	case 0:
+	  N = COSSUN;
+	  semidiameter = 0.2666666666666667;
+	  elevation_threshold = -0.8333333333333333;
+	  break;
 
 /* Moon, elevation = -34' - semidiameter + parallax
  * semidiameter = 0.272453 * parallax + 0.0799"
@@ -80,20 +93,25 @@
 	case 3:
 		N = 1.0/(DISFAC*obpolar[2]);
 		D = asin( N ); /* the parallax */
-		N =  - 9.890199094634534e-3
-			+ (1.0 - 0.2725076)*D
-			- 3.874e-7;
+		semidiameter = 0.2725076*D + 3.874e-7;
+		N =  -9.890199094634534e-3 - semidiameter + D;
+		semidiameter *= RTD;
+		elevation_threshold = -34.0/60.0 - semidiameter;
 		N = sin(N);
 		break;
 
 /* Other object */
-	default: N = COSZEN; break;
+	default:
+	  N = COSZEN;
+	  semidiameter = 0.0;
+	  elevation_threshold = -0.5666666666666666;
+	  break;
 	}
 y = (N - sinlat*sindec)/(coslat*cosdec);
 
 if( (y < 1.0) && (y > -1.0) )
 	{
-	f_trnsit = 1;
+	f_trnsit = 3;
 /* Derivative of y with respect to declination
  * times rate of change of declination:
  */
@@ -105,7 +123,7 @@
 	D = -dradt/TPI + 1.00273790934;
 	r_rise = x - (lhay + y)*(1.0 + z)/D;
 	r_set = x - (lhay - y)*(1.0 - z)/D;
-	/* Ordinarily never print here.  */
+	/* Ordinarily do not print here.  */
 		if( prtflg > 1 )
 		{
 		printf( "rises approx " );
@@ -128,8 +146,8 @@
 /* Julian dates of rise, transet and set times.  */
 double t_rise;
 double t_trnsit;
+double elevation_trnsit;
 double t_set;
-
 /* Compute estimate of lunar rise and set times for iterative solution.  */
 
 static void
@@ -156,8 +174,8 @@
 int (* func)();
 {
   double JDsave, TDTsave, UTsave;
-  double date, date_save, date_trnsit, t0, t1;
-  double rise1, set1, trnsit1, loopctr, retry;
+  double date, date_trnsit, t0, t1;
+  double rise1, set1, trnsit1, loopctr;
   double TPI;
   int prtsave;
 
@@ -167,16 +185,14 @@
   JDsave = JD;
   TDTsave = TDT;
   UTsave = UT;
-  date = floor(UT - 0.5) + 0.5;
-  retry = 0;
-  date_save = date;
-  t1 = date_save;
+  /* Start iteration at time given by the user.  */
+  t1 = UT;
 
   /* Find transit time. */
   do
     {
-      date = date_save;
       t0 = t1;
+      date = floor(t0 - 0.5) + 0.5;
       iter_func(t0, func);
       t1 = date + r_trnsit / TPI;
       if (++loopctr > 10)
@@ -184,35 +200,25 @@
 	  printf ("? Transit time did not converge.\n");
 	  goto no_trnsit;
 	}
-
-  /* Reject transit if it is more than half a day from entered date.  */
-  if (retry == 0)
-    {
-      if ((UTsave - t1) > 0.5)
-	{
-	  date_save += 1.0;
-	  t1 = t_trnsit + 1.0;
-	  retry = 1;
-	  /*	  goto do_retry; */
-	}
-      if ((UTsave - t1) < -0.5)
-	{
-	  date_save -= 1.0;
-	  t1 = t_trnsit - 1.0;
-	  retry = 1;
-	  /*	  goto do_retry; */
-	}
-    }
     }
   while (fabs(t1 - t0) > .0001);
 
   t_trnsit = t1;
+  elevation_trnsit = elevation;
   trnsit1 = r_trnsit;
   set1 = r_set;
 
 
-  if (f_trnsit == 0)
-	goto prtrnsit;
+  if (f_trnsit != 3)
+    {
+      /* Rise or set time not found.  Apply a search technique
+	 if object is above horizon now.  */
+      t_rise = -1.0;
+      t_set = -1.0;
+      if (elevation > elevation_threshold)
+	no_rise_set (t_trnsit, func);
+      goto prtrnsit;
+    }
 
   /* Set current date to be that of the transit just found.  */
   date_trnsit = date;
@@ -228,14 +234,20 @@
     {
       t0 = t1;
       iter_func(t0, func);
-      /* Skip out if no event found.  Don't report rise or set.  */
-      if ((f_trnsit == 0) || (++loopctr > 10))
+      /* Skip out if no event found.  */
+
+      if ((f_trnsit & 1) == 0)
+	{
+	  /* Rise or set time not found.  Apply search technique.  */
+	  t_rise = -1.0;
+	  t_set = -1.0;
+	  no_rise_set (t_trnsit, func);
+	  goto prtrnsit;
+	}
+      if (++loopctr > 10)
 	{
-	  if (loopctr > 10)
-	    {
-	      printf ("? Rise time did not converge.\n");
-	    }
-	  f_trnsit = 0;
+	  printf ("? Rise time did not converge.\n");
+	  f_trnsit &= ~1;
 	  goto prtrnsit;
 	}
       t1 = date + r_rise / TPI;
@@ -264,13 +276,18 @@
     {
       t0 = t1;
       iter_func(t0, func);
-      if ((f_trnsit == 0) || (++loopctr > 10))
+      if ((f_trnsit & 2) == 0)
+	{
+	  /* Rise or set time not found.  Apply search technique.  */
+	  t_rise = -1.0;
+	  t_set = -1.0;
+	  no_rise_set (t_trnsit, func);
+	  goto prtrnsit;
+	}
+      if (++loopctr > 10)
 	{
-	  if (loopctr > 10)
-	    {
-	      printf ("? Set time did not converge.\n");
-	    }
-	  f_trnsit = 0;
+	  printf ("? Set time did not converge.\n");
+	  f_trnsit &= ~2;
 	  goto prtrnsit;
 	}
       t1 = date + r_set / TPI;
@@ -291,21 +308,29 @@
   printf( "local meridian transit " );
   UT = t_trnsit;
   jtocal (t_trnsit);
-  if (f_trnsit != 0)
+  if (f_trnsit & 1)
     {
       printf( "rises " );
       UT = t_rise;
       jtocal (t_rise);
+    }
+  if (f_trnsit & 2)
+    {
       printf( "sets " );
       UT = t_set;
       jtocal (t_set);
-      t0 = t_set - t_rise;
-      if ((t0 > 0.0) && (t0 < 1.0))
-	printf("Visible hours %.4f\n", 24.0 * t0);
-      if ((fabs(JDsave - t_rise) > 0.5) && (fabs(JDsave - t_trnsit) > 0.5)
-	  && (fabs(JDsave - t_set) > 0.5))
-	printf("? wrong event date.\n");
+      if (t_rise != -1.0)
+	{
+	  t0 = t_set - t_rise;
+	  if ((t0 > 0.0) && (t0 < 1.0))
+	    printf("Visible hours %.4f\n", 24.0 * t0);
+	}
     }
+
+  if ((fabs(JDsave - t_rise) > 0.5) && (fabs(JDsave - t_trnsit) > 0.5)
+      && (fabs(JDsave - t_set) > 0.5))
+    printf("? wrong event date.\n");
+
 no_trnsit:
   JD = JDsave;
   TDT = TDTsave;
@@ -314,5 +339,157 @@
   prtflg = 0;
   update();
   prtflg = prtsave;
+  f_trnsit = 3;
   return 0;
 }
+
+
+/* If the initial approximation fails to locate a rise or set time,
+   this function steps between the transit time and the previous
+   or next inferior transits to find an event more reliably.  */
+
+int no_rise_set (t0, func)
+     double t0;
+     int (* func)();
+{
+  double t_trnsit0 = t_trnsit;
+  double el_trnsit0 = elevation_trnsit;
+  double t_above, el_above, t_below, el_below;
+  double el1, el2, t1, t2;
+  
+  /* Step time toward previous inferior transit to find
+     whether a rise event was missed.  The step size is a function
+     of the azimuth and decreases near the transit time.  */
+  t_above = t_trnsit0;
+  el_above = el_trnsit0;
+  t_below = -1.0;
+  el_below = el_above;
+
+  t1 = t_trnsit0;
+  t2 = t1;
+  el1 = el_trnsit0;
+
+  for (;;)
+    {
+      if ((t_trnsit0 - t2) > 1.0)
+	break;
+      t2 = t1 - 0.25;
+      iter_func(t2, func);
+      el2 = elevation;
+
+      if (el2 < elevation_threshold)
+	{
+	  /* Object is below horizon at t2.  Rise event is bracketed.
+	     Proceed to interval halving search.  */
+	  t_below = t2;
+	  el_below = elevation;
+	  t_above = t1;
+	  el_above = el1;
+	  t_rise = search_halve (t_below, el_below,
+				 t_above, el_above, func);
+	  f_trnsit |= 1;
+	  break;
+	}
+      if (el2 < el1)
+	{
+	  /* decreasing elevation, continue */
+	  t1 = t2;
+	  el1 = el2;
+	  continue;
+	}
+      /* Non-decreasing elevation. Stop. */
+      break;
+    }
+  if (elevation > elevation_threshold)
+    {
+      /* No rise event detected.  */
+      /* printf ("Previous inferior transit is above horizon.\n"); */
+      f_trnsit &= ~1;
+    }
+
+
+  /* Step forward in time toward the next inferior transit.  */
+  t1 = t_trnsit0;
+  t2 = t1;
+  el1 = el_trnsit0;
+
+  for (;;)
+    {
+      if ((t2 - t_trnsit0) > 1.0)
+	break;
+      t2 = t1 + 0.25;
+      iter_func(t2, func);
+      el2 = elevation;
+
+      if (el2 < elevation_threshold)
+	{
+	  /* Object is below horizon at t2.  Set event is bracketed.
+	     Proceed to interval halving search.  */
+	  t_below = t2;
+	  el_below = elevation;
+	  t_above = t1;
+	  el_above = el1;
+	  t_set = search_halve (t_below, el_below,
+				 t_above, el_above, func);
+	  f_trnsit |= 2;
+	  break;
+	}
+      if (el2 < el1)
+	{
+	  /* decreasing elevation, continue */
+	  t1 = t2;
+	  el1 = el2;
+	  continue;
+	}
+      /* Non-decreasing elevation. Stop. */
+      break;
+    }
+
+  if (elevation > elevation_threshold)
+    {
+      /* No set event detected.  */
+      /* printf ("Next inferior transit is above horizon.\n"); */
+      f_trnsit &= ~2;
+    }
+  return 0;
+}
+
+
+/* Search rise or set time by simple interval halving
+   after the event has been bracketed in time.  */
+
+double
+search_halve (t1, y1, t2, y2, func)
+     double t1;
+     double y1;
+     double t2;
+     double y2;
+     int (* func)();
+{
+  double e2, em, tm, ym;
+
+  e2 = y2 - elevation_threshold;
+  tm = 0.5 * (t1 + t2);
+
+while( fabs(t2 - t1) > .00001 )
+  {
+    /* Evaluate at middle of current interval.  */
+    tm = 0.5 * (t1 + t2);
+    iter_func(tm, func);
+    ym = elevation;
+    em = ym - elevation_threshold;
+    /* Replace the interval boundary whose error has the same sign as em.  */
+     if( em * e2 > 0 )
+      {
+	y2 = ym;
+	t2 = tm;
+	e2 = em;
+      }
+    else
+      {
+	y1 = ym;
+	t1 = tm;
+      }
+  }
+return (tm);
+}