1 /*----------------------------------------------------------------------------+
2  | Sunrise/Sunset functions for Heyu.                                         |
3  |                                                                            |
4  | The solar computation functions below were programmed by Charles W.        |
5  | Sullivan utilizing the techniques and astronomical constants published by  |
6  | Roger W. Sinnott in the August 1994 issue of Sky & Telescope Magazine,     |
7  | Page 84.                                                                   |
8  +----------------------------------------------------------------------------*/
9 
10 /*
11  *   This program is free software: you can redistribute it and/or modify
12  *   it under the terms of the GNU General Public License as published by
13  *   the Free Software Foundation, either version 3 of the License, or
14  *   (at your option) any later version.
15  *
16  *   This program is distributed in the hope that it will be useful,
17  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
18  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19  *   GNU General Public License for more details.
20  *
21  *   You should have received a copy of the GNU General Public License
22  *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
23  *
24  */
25 
26 #include <stdlib.h>
27 #include <stdio.h>
28 #include <string.h>
29 #include <math.h>
30 #include <time.h>
31 
32 #include "sun.h"
33 #include "process.h"
34 
35 #define PI    3.14159265
36 #define D2R   (PI / 180.)
37 
38 /* Ratio, Mean Solar Day / Mean Siderial Day */
39 #define MSSR   1.0027379
40 
41 #define KS15   (15. * D2R * MSSR)
42 
43 /* Local functions */
44 static void sun_position ( double, double *, double *);
45 static double local_sidereal_time (double, long int, double);
46 
47 /*---------------------------------------------------------+
48  | Calculate local times and Azimuths of Sunrise and       |
49  | Sunset on a specified date.                             |
50  | Input:                                                  |
51  |   latitude and longitude in degrees are positive for    |
52  |     North and East of Greenwich respectively.           |
53  |   JulianDay is the Julian Day number at Greenwich Noon. |
54  |   timezone in seconds from Greenwich, positive for      |
55  |     localities west of Greenwich.                       |
56  | Return:                                                 |
57  |   sunrise and sunset times are in minutes after local   |
58  |     midnight.                                           |
59  |   azimuths at rise and set are in degrees (0 - 360).    |
60  |   return code (defined in sun.h):                       |
61  |     NORMAL_SUN    Sun rises and sets on this day.       |
62  |     DOWN_ALL_DAY  Sun below horizon all day.            |
63  |     UP_ALL_DAY    Sun above horizon all day.            |
64  |     NO_SUNRISE    Sun does not rise on this day.        |
65  |     NO_SUNSET     Sun does not set on this day.         |
66  |                                                         |
67  +---------------------------------------------------------*/
68 
suntimes(double latitude,double longitude,long int timezone,long int JulianDay,int sunmode,int offset,int * sunrise,int * sunset,double * azrise,double * azset)69 int suntimes ( double latitude, double longitude, long int timezone,
70                long int JulianDay,
71 	       int sunmode, int offset,
72                int *sunrise, int *sunset,
73                double *azrise, double *azset)
74 {
75    double    tdays, lst, zendist;
76    double    slat, clat,cozend ;
77    double    rasc[2], decl[2];
78    double    hj ;
79    int       flag_rise, flag_set;
80    int       j ;
81    int       retcode = NORMAL_SUN;
82    int       *srssptr = NULL;
83    double    *azptr = NULL;
84    double    a, b, a0, a2, d0, d1, d2, dela, deld, p,
85              el0, el2, h0, h1, h2, v0, v1, v2,
86              d, e, t3, h7, n7, d7, az;
87 
88 /* Zenith angles for Sunrise/set and Twilights */
89 #define ANG_RISESET  (90. + 50./60.)
90 #define ANG_CIVILTWI (96.)
91 #define ANG_NAUTITWI (102.)
92 #define ANG_ASTROTWI (108.)
93 #define ANG_ANGOFFS  (90. + offset/60.)
94 
95    double zenangle[] = {
96       ANG_RISESET, ANG_CIVILTWI, ANG_NAUTITWI, ANG_ASTROTWI, ANG_ANGOFFS, };
97 
98 
99    *sunrise = *sunset = 0 ;
100 
101    /* Elapsed days from 1 Jan 2000 at 00:00 hours UTC0 */
102    tdays = (double)(JulianDay - 2451545L) - 0.5 ;
103 
104    /* Determine Local Siderial Time. */
105    lst = local_sidereal_time( tdays, timezone, longitude );
106 
107    tdays += (double)timezone / (3600. * 24.);
108 
109    /* Get sun's position */
110    for ( j = 0; j < 2; j++ ) {
111       sun_position ( tdays, &rasc[j], &decl[j] );
112       tdays += 1.0;
113    }
114 
115    if ( rasc[1] < rasc[0] )
116       rasc[1] += 2.*PI ;
117 
118 #if 0
119    /* sunrise and sunset are defined when the sun is
120       50 minutes below the horizon (at sea level).     */
121    zendist = D2R * (90. + 50./60.) ;
122 #endif
123    zendist = D2R * zenangle[sunmode];
124 
125    slat = sin(D2R * latitude);
126    clat = cos(D2R * latitude);
127    cozend = cos(zendist);
128 
129    flag_rise = flag_set = 0;
130 
131    a0 = rasc[0];
132    d0 = decl[0];
133    v0 = 0.0;
134    dela = rasc[1] - rasc[0];
135    deld = decl[1] - decl[0];
136 
137    for ( j = 0; j < 24; j++ ) {
138       hj = (double) j ;
139       p = (1.0 + hj) / 24. ;
140 
141       a2 = rasc[0] + p * dela ;
142       d2 = decl[0] + p * deld ;
143 
144       /* Test an hour for an event */
145       el0 = lst + hj * KS15 ;
146       el2 = el0 + KS15 ;
147       h0 = el0 - a0;
148       h2 = el2 - a2 ;
149       h1 = 0.5 * (h2 + h0);
150       d1 = 0.5 * (d2 + d0) ;
151 
152       if ( j == 0 ) {
153          v0 = slat * sin(d0) + clat * cos(d0) * cos(h0) - cozend ;
154       }
155       v2 = slat * sin(d2) + clat * cos(d2) * cos(h2) - cozend ;
156       if ( v0 * v2 < 0. ) {
157          v1 = slat * sin(d1) + clat * cos(d1) * cos(h1) - cozend ;
158          a = 2. * v2 - 4. * v1 + 2. * v0 ;
159          b = 4. * v1 - 3. * v0 - v2 ;
160          d = b * b - 4. * a * v0 ;
161          if ( d >= 0. ) {
162             d = sqrt(d);
163             if ( v0 < 0. && v2 > 0. ) {
164                /* Event is Sunrise */
165                srssptr = sunrise ;
166                azptr = azrise ;
167                flag_rise = 1 ;
168             }
169             if ( v0 > 0. && v2 < 0. ) {
170                /* Event is Sunset */
171                srssptr = sunset ;
172                azptr = azset ;
173                flag_set = 1;
174             }
175             e = (-b + d) / (2. * a);
176             if ( e > 1. || e < 0. )
177                e = (-b -d) / (2. * a) ;
178             t3 = hj + e ;
179             *srssptr = (int)(60.* t3 + 0.5) ;  /* Round off*/
180 
181             h7 = h0 + e * (h2 -h0);
182             n7 = - cos(d1) * sin(h7);
183             d7 = clat * sin(d1) - slat * cos(d1) * cos(h7);
184 
185             az = fmod(((atan2(n7, d7) / D2R) + 360.), 360.) ;
186             if ( azptr != NULL )
187                *azptr = az ;
188          }
189       }
190       if ( flag_rise & flag_set )
191          break;
192       a0 = a2;
193       d0 = d2;
194       v0 = v2;
195    }
196    if ( !(flag_rise | flag_set) ) {
197       if ( v2 < 0. ) {
198          /* Sun down all day */
199          retcode = DOWN_ALL_DAY ;
200       }
201       else if ( v2 >= 0. ) {
202          /* Sun up all day */
203          retcode = UP_ALL_DAY ;
204       }
205    }
206    else if ( !flag_rise ) {
207       /* No sunrise this date */
208       retcode = NO_SUNRISE ;
209    }
210    else if ( !flag_set ) {
211       /* No sunset this date */
212       retcode = NO_SUNSET ;
213    }
214 
215    return retcode;
216 }
217 
218 /*---------------------------------------------------------+
219  | Calculate local sidereal time.                          |
220  | Input:                                                  |
221  |   tday - number of days from 2000 Jan 1 at              |
222  |     00:00 hours at Greenwich (UTC0).                    |
223  +---------------------------------------------------------*/
224 
local_sidereal_time(double tday,long int timezone,double longitude)225 double local_sidereal_time ( double tday, long int timezone, double longitude )
226 {
227    double s ;
228 
229    s = 24110.5 + 8640184.812999999*tday/36525. ;
230    s += 86636.6*(double)timezone/(3600.*24.) + 86400.*longitude/360. ;
231    s = s/86400. ;
232    s = s - floor(s) ;
233 
234    return  (s * 360. * D2R) ;
235 }
236 
237 
238 /*---------------------------------------------------------+
239  | Calculate sun's Right Ascention and Declination.        |
240  +---------------------------------------------------------*/
241 
sun_position(double tday,double * rasc,double * decl)242 void sun_position ( double tday, double *rasc, double *decl )
243 {
244    double l, g, s, u, v, w ;
245    double tcent ;
246 
247    /* Julian centuries from 1900.0 */
248    tcent = tday / 36525. + 1.0 ;
249 
250    /* Fundamental arguments (Van Flandern & Pulkinnen, 1979) */
251 
252    l = .779072 + .00273790931 * tday ;
253    g = .993126 + .0027377785 * tday ;
254    l = 2. * PI * (l - floor(l)) ;
255    g = 2. * PI * (g - floor(g)) ;
256 
257    v = .39785 * sin(l) - .01 * sin(l - g) +
258        .00333 * sin(l + g) - .00021 * sin(l) * tcent ;
259    u = 1.0 - .03349 * cos(g) - .00014 * cos(2. * l) +
260        .00008 * cos (l) ;
261    w = - .0001 - .04129 * sin(2. * l) +
262        .03211 * sin(g) + .00104 * sin(2. * l - g) -
263        .00035 * sin( 2. * l + g) - .00008 * sin(g) * tcent ;
264    s = w / sqrt(u - v * v);
265    *rasc = l + atan(s / sqrt(1.0 - s * s));
266    s = v / sqrt(u);
267    *decl = atan(s / sqrt(1.0 - s * s));
268 
269    return ;
270 }
271 
272 /*---------------------------------------------------------+
273  | Gregorian calendar day to Julian Day number.            |
274  | Returns the Julian Day number at Greenwich Noon         |
275  | (UTC0 12:00) for the year, month, and day arguments.    |
276  |                                                         |
277  | The Julian Day number is the count of whole days from   |
278  | Noon on 1 Jan 4713 B.C.E. in the Julian Proleptic       |
279  | Calendar.  This calculation is historically valid from  |
280  | 15 Oct 1582 onward, however any year, month, and day    |
281  | greater than zero are acceptable as arguments and will  |
282  | yield the logically correct result.                     |
283  +---------------------------------------------------------*/
284 
greg2jd(int year,int month,int day)285 long int greg2jd( int year, int month, int day )
286 {
287    int count ;
288 
289    if ( month > 12 ) {
290       year += ( month - 1 ) / 12 ;
291       month = ( month - 1 ) % 12 + 1 ;
292    }
293 
294    if ( month > 2 )
295       count = - ( 4 * month + 23 ) / 10 ;
296    else  {
297       count = 365 ;
298       year-- ;
299    }
300 
301    count = count + year / 4
302                  - 3 * ( year / 100 + 1 ) / 4
303                  + 31 * ( month - 1 ) + day ;
304 
305    return ( (long)count + 365L*(long)year + 1721060L ) ;
306 }
307 
308 /*---------------------------------------------------------------------+
309  | Adjust the times of Dawn and Dusk for abnormal sun conditions, as   |
310  | in polar regions, by creating and artificial Dawn and/or Dusk at    |
311  | 00:01 or 23:58 as appropriate.                                      |
312  +---------------------------------------------------------------------*/
abnormal_sun_adjust(int scode,int * dawnp,int * duskp)313 int abnormal_sun_adjust ( int scode, int *dawnp, int *duskp )
314 {
315    switch ( scode ) {
316       case NORMAL_SUN :
317          break;
318       case UP_ALL_DAY :
319          *dawnp = 1;
320          *duskp = 23 * 60 + 58;
321          break;
322       case DOWN_ALL_DAY :
323          *dawnp = 23 * 60 + 58;
324          *duskp = 1;
325          break;
326       case NO_SUNRISE :
327          *dawnp = 1;
328          break;
329       case NO_SUNSET :
330          *duskp = 23 * 60 + 58;
331          break;
332       default :
333          break;
334    }
335    return scode;
336 }
337 
338 /*---------------------------------------------------------------------+
339  | Return the Julian Day corresponding to UTC0 Noon for the UTC0 time  |
340  | argument (in seconds from 1/1/1970 at 00:00:00 UTC0).               |
341  +---------------------------------------------------------------------*/
utc2jd(long lutc0)342 long int utc2jd ( long lutc0 )
343 {
344    extern CONFIG *configp;
345 
346    return (lutc0 - ((lutc0 - configp->tzone) % 86400L)) / 86400L + 2440588L;
347 }
348 
349 /*---------------------------------------------------------------------+
350  | Compute the UTC0 times of local Dawn and Dusk for the day which     |
351  | includes the argument UTC0 time.  UTC0 times are all expressed as   |
352  | seconds elapsed from 1/1/1970 at 00:00:00 UTC0.                     |
353  +---------------------------------------------------------------------*/
local_dawndusk(time_t utc0,time_t * utc0_dawn,time_t * utc0_dusk)354 int local_dawndusk( time_t utc0, time_t *utc0_dawn, time_t *utc0_dusk )
355 {
356    long int jd;
357    int      dawn, dusk, scode;
358    long     lutc0, midnight;
359 
360    extern CONFIG   *configp;
361 
362    if ( configp->loc_flag != (LATITUDE | LONGITUDE) ) {
363       *utc0_dawn = *utc0_dusk = 0;
364       return -1;
365    }
366 
367    lutc0 = (long)utc0;
368 
369    midnight = lutc0 - ((lutc0 - configp->tzone) % 86400L);
370 
371    jd = utc2jd(lutc0);
372 
373    scode = suntimes(configp->latitude, configp->longitude, configp->tzone, jd,
374                configp->sunmode, configp->sunmode_offset, &dawn, &dusk, NULL, NULL);
375 
376    /* Adjust for abnormal sun conditions */
377    abnormal_sun_adjust(scode, &dawn, &dusk);
378 
379    *utc0_dawn = (time_t)(midnight + 60L * (long)dawn);
380    *utc0_dusk = (time_t)(midnight + 60L * (long)dusk);
381 
382    return scode;
383 }
384 
385 
386 /*-------------------------------------------------------------+
387  | Display a table of Dawn and Dusk for the year in the        |
388  | format used by the US Naval Observatory website at the time |
389  | of this writing at:                                         |
390  |    http://aa.usno.navy.mil/data/docs/RS_OneYear.html        |
391  | (Standard time only is displayed.)                          |
392  | Argument sunmode indicates the definition of Dawn and Dusk: |
393  |   0  -> Rise and Set                                        |
394  |   1  -> Civil Twilight                                      |
395  |   2  -> Nautical Twilight                                   |
396  |   3  -> Astronomical Twilight                               |
397  | Printing on letter size or A4 paper requires 8 point fixed  |
398  | font and landscape mode.                                    |
399  +-------------------------------------------------------------*/
display_sun_table_wide(FILE * fd_sun,int year,long timezone,int sunmode,int offset,int timemode,int lat_d,int lat_m,int lon_d,int lon_m)400 int display_sun_table_wide ( FILE *fd_sun, int year, long timezone,
401     int sunmode, int offset, int timemode, int lat_d, int lat_m, int lon_d, int lon_m )
402 {
403    static struct tzones {
404       char *name;
405       long seconds;
406    } us_tzones[] = {
407      {"Atlantic",         14400 },
408      {"Eastern",          18000 },
409      {"Central",          21600 },
410      {"Mountain",         25200 },
411      {"Pacific",          28800 },
412      {"Alaska",           32400 },
413      {"Hawaii-Aleutian",  36000 },
414      {"Samoa",            39600 },
415      {"Wake Island",     -43200 },
416      {"Guam",            -39600 },
417    };
418    static int n_tzones = ( sizeof(us_tzones) / sizeof(struct tzones) );
419 
420    static int mdays[] = {0,31,28,31,30,31,30,31,31,30,31,30,31};
421 
422    int time_adjust ( int, int, unsigned char );
423 
424    long   julianday, julianday0;
425 
426    int    j, month, day, yday ;
427    int    rise, set, scode;
428    int    retcode;
429    char   tmark = ' ';
430    char   *timename;
431 
432    double latitude, longitude;
433 
434    get_dst_info(year);
435 
436    timename = ( timemode == TIMEMODE_CIVIL ) ? "Civil" : "Standard";
437 
438    latitude  = (lat_d < 0) ? (double)lat_d - (double)lat_m / 60.  :
439                              (double)lat_d + (double)lat_m / 60. ;
440    longitude = (lon_d < 0) ? (double)lon_d - (double)lon_m / 60.  :
441                              (double)lon_d + (double)lon_m / 60. ;
442 
443    julianday0 = greg2jd(year, 1, 1);
444 
445    if ( greg2jd(year + 1, 1, 1) - julianday0 == 366 )
446       mdays[2] = 29;
447    else
448       mdays[2] = 28;
449 
450    (void) fprintf(fd_sun, "\n              o  ,    o  ,  \n");
451    (void) fprintf(fd_sun, "Location: ");
452    (void) fprintf(fd_sun, "%s%03d %02d,",
453         longitude < 0 ? "W" : "E", abs(lon_d), abs(lon_m));
454    (void) fprintf(fd_sun, " %s%02d %02d",
455         latitude  < 0 ? "S" : "N", abs(lat_d), abs(lat_m));
456 
457    switch (sunmode) {
458       case RiseSet :
459          (void) fprintf(fd_sun, "%*sSunrise and Sunset for %d", 29, " ", year);
460          break;
461       case CivilTwi :
462          (void) fprintf(fd_sun, "%*sCivil Twilight for %d", 29, " ", year);
463          break;
464       case NautTwi :
465          (void) fprintf(fd_sun, "%*sNautical Twilight for %d", 29, " ", year);
466          break;
467       case AstroTwi :
468          (void) fprintf(fd_sun, "%*sAstronomical Twilight for %d", 26, " ", year);
469          break;
470       case AngleOffset :
471          (void) fprintf(fd_sun, "%*sSun centre at %d angle minutes below horizon for %d", 26, " ", offset, year);
472          break;
473    };
474    (void) fprintf(fd_sun, "%*sHEYU ver 2.0 \n\n", 39, " ");
475 
476 
477    for ( j = 0; j < n_tzones; j++ ) {
478       if ( us_tzones[j].seconds == timezone )
479          break;
480    }
481    if ( j < n_tzones )
482       (void) fprintf(fd_sun, "%*sUS/%s %s Time\n\n",
483                                             54, " ", us_tzones[j].name, timename);
484    else
485       (void) fprintf(fd_sun, "%*sTimezone: %.1fh %s of Greenwich - %s Time\n\n", 51, " ",
486             (double)abs(timezone)/3600., (timezone < 0 ? "East" : "West"), timename );
487 
488    (void) fprintf(fd_sun, "       Jan        Feb        Mar        Apr ");
489    (void) fprintf(fd_sun, "       May        Jun        Jul        Aug ");
490    (void) fprintf(fd_sun, "       Sep        Oct        Nov        Dec \n");
491    (void) fprintf(fd_sun, "Day Rise  Set");
492 
493    for ( month = 2; month <= 12; month++ )
494       (void) fprintf(fd_sun, "  Rise  Set");
495    (void) fprintf(fd_sun, "\n  ");
496 
497    for ( month = 1; month <= 12; month++ )
498       (void) fprintf(fd_sun, "   h m  h m");
499 
500    scode = NORMAL_SUN;
501    for ( day = 1; day < 32; day++ ) {
502       (void) fprintf(fd_sun, "\n%02d", day);
503       for ( month = 1; month <= 12; month++ ) {
504          if ( day > mdays[month] ) {
505             (void) fprintf(fd_sun, "           ");
506             continue;
507          }
508          julianday = greg2jd( year, month, day );
509 
510          retcode = suntimes(latitude, longitude, timezone, julianday,
511                  sunmode, offset, &rise, &set, NULL, NULL );
512 
513          yday = (int)(julianday - julianday0);
514 
515          if ( timemode == TIMEMODE_CIVIL ) {
516             /* Adjust times for Daylight Time */
517             if ( retcode == NORMAL_SUN || retcode == NO_SUNSET ) {
518                rise += time_adjust(yday, rise, LGL2STD);
519                if ( rise < 0 || rise > 1439 )
520                   retcode = NO_SUNRISE;
521             }
522             if ( retcode == NORMAL_SUN || retcode == NO_SUNRISE ) {
523                set += time_adjust(yday, set, LGL2STD);
524                if ( set < 0 || set > 1439 )
525                   retcode = NO_SUNSET;
526             }
527 
528             /* Mark day of time change */
529             tmark = (yday == 0 ||
530                   time_adjust(yday, 720, LGL2STD) == time_adjust(yday - 1, 720, LGL2STD)) ? ' ' : '*';
531          }
532 
533          scode |= retcode;
534 
535          switch ( retcode ) {
536             case DOWN_ALL_DAY :
537                (void) fprintf(fd_sun, " %c---- ----", tmark);
538                break ;
539             case UP_ALL_DAY :
540                (void) fprintf(fd_sun, " %c**** ****", tmark);
541                break ;
542             case NO_SUNRISE :
543                (void) fprintf(fd_sun, " %c     %02d%02d", tmark, set/60, set % 60);
544                break ;
545             case NO_SUNSET :
546                (void) fprintf(fd_sun, " %c%02d%02d     ", tmark, rise/60, rise%60);
547                break ;
548             default :
549                (void) fprintf(fd_sun, " %c%02d%02d %02d%02d",
550                             tmark, rise/60, rise%60, set/60,set%60 );
551          }
552       }
553    }
554 
555    if ( timemode == TIMEMODE_CIVIL ) {
556       fprintf(fd_sun, "\n\n(*) Denotes time change");
557    }
558    else {
559       fprintf(fd_sun,
560          "\n\n%*sAdd offset for Daylight Time (usually +60 minutes), if and when in effect.", 34, " ");
561    }
562 
563    if ( scode & (DOWN_ALL_DAY | UP_ALL_DAY) ) {
564       fprintf(fd_sun, "\n\n(****) Sun continuously above horizon");
565       fprintf(fd_sun,  "%*s(----) Sun continuously below horizon", 60, " ");
566    }
567 
568    fprintf(fd_sun, "\n");
569 
570    return 0;
571 }
572 
573 
574 /*-------------------------------------------------------------+
575  | Display a table of Dawn and Dusk for the year.              |
576  | Argument sunmode indicates the definition of Dawn and Dusk: |
577  |   0  -> Rise and Set                                        |
578  |   1  -> Civil Twilight                                      |
579  |   2  -> Nautical Twilight                                   |
580  |   3  -> Astronomical Twilight                               |
581  +-------------------------------------------------------------*/
display_sun_table(FILE * fd_sun,int year,long timezone,int sunmode,int offset,int timemode,int lat_d,int lat_m,int lon_d,int lon_m)582 int display_sun_table ( FILE *fd_sun, int year, long timezone,
583     int sunmode, int offset, int timemode, int lat_d, int lat_m, int lon_d, int lon_m )
584 {
585    static struct tzones {
586       char *name;
587       long seconds;
588    } us_tzones[] = {
589      {"Atlantic",         14400 },
590      {"Eastern",          18000 },
591      {"Central",          21600 },
592      {"Mountain",         25200 },
593      {"Pacific",          28800 },
594      {"Alaska",           32400 },
595      {"Hawaii-Aleutian",  36000 },
596      {"Samoa",            39600 },
597      {"Wake Island",     -43200 },
598      {"Guam",            -39600 },
599    };
600    static int n_tzones = ( sizeof(us_tzones) / sizeof(struct tzones) );
601 
602    static int mdays[] = {0,31,28,31,30,31,30,31,31,30,31,30,31};
603 
604    int time_adjust ( int, int, unsigned char );
605 
606    long   julianday, julianday0;
607 
608    int    j, month, day, yday, half ;
609    int    rise, set, scode;
610    int    retcode;
611    char   tmark = ' ';
612    char   minibuf[64];
613    char   *timename;
614 
615    double latitude, longitude;
616 
617    get_dst_info(year);
618 
619    timename = ( timemode == TIMEMODE_CIVIL ) ? "Civil" : "Standard";
620 
621    latitude  = (lat_d < 0) ? (double)lat_d - (double)lat_m / 60.  :
622                              (double)lat_d + (double)lat_m / 60. ;
623    longitude = (lon_d < 0) ? (double)lon_d - (double)lon_m / 60.  :
624                              (double)lon_d + (double)lon_m / 60. ;
625 
626    julianday0 = greg2jd(year, 1, 1);
627 
628    if ( greg2jd(year + 1, 1, 1) - julianday0 == 366 )
629       mdays[2] = 29;
630    else
631       mdays[2] = 28;
632 
633    switch (sunmode) {
634       case RiseSet :
635          sprintf(minibuf, "Sunrise and Sunset for %d\n", year);
636          break;
637       case CivilTwi :
638          sprintf(minibuf, "Civil Twilight for %d\n", year);
639          break;
640       case NautTwi :
641          sprintf(minibuf, "Nautical Twilight for %d\n", year);
642          break;
643       case AstroTwi :
644          sprintf(minibuf, "Astronomical Twilight for %d\n", year);
645          break;
646       case AngleOffset :
647          sprintf(minibuf, "Sun centre at %d angle minutes below horizon for %d\n", offset, year);
648          break;
649    };
650    fprintf(fd_sun, "%*s%s", (80 - (int)strlen(minibuf))/2, " ", minibuf);
651    (void) fprintf(fd_sun, "\nLocation: ");
652    (void) fprintf(fd_sun, "%s%03d:%02d,",
653         longitude < 0 ? "W" : "E", abs(lon_d), abs(lon_m));
654    (void) fprintf(fd_sun, " %s%02d:%02d",
655         latitude  < 0 ? "S" : "N", abs(lat_d), abs(lat_m));
656 
657    for ( j = 0; j < n_tzones; j++ ) {
658       if ( us_tzones[j].seconds == timezone )
659          break;
660    }
661    if ( j < n_tzones )
662       (void) fprintf(fd_sun, "   US/%s %s Time\n\n", us_tzones[j].name, timename);
663    else
664       (void) fprintf(fd_sun, "   Timezone: %.1fh %s of Greenwich - %s Time\n\n",
665             (double)abs(timezone)/3600., (timezone < 0 ? "East" : "West"), timename );
666 
667    for ( half = 0; half < 2; half++ ) {
668 
669       if ( half == 0 )
670          (void) fprintf(fd_sun, "       Jan          Feb          Mar          Apr          May          Jun\n");
671       else
672          (void) fprintf(fd_sun, "\n\n       Jul          Aug          Sep          Oct          Nov          Dec\n");
673 
674       if ( sunmode == 0 ) {
675          fprintf(fd_sun, "Day Rise   Set");
676          for ( month = 2; month <= 6; month++ )
677             fprintf(fd_sun, "   Rise   Set");
678          fprintf(fd_sun, "\n  ");
679       }
680       else {
681          fprintf(fd_sun, "Day Morn   Eve");
682          for ( month = 2; month <= 6; month++ )
683             fprintf(fd_sun, "   Morn   Eve");
684          fprintf(fd_sun, "\n  ");
685       }
686 
687       for ( month = 1; month <= 6; month++ )
688          (void) fprintf(fd_sun, "  hh:mm hh:mm");
689 
690       scode = NORMAL_SUN;
691       for ( day = 1; day < 32; day++ ) {
692          (void) fprintf(fd_sun, "\n%02d", day);
693          for ( month = 6 * half + 1; month <= 6 * half + 6; month++ ) {
694             if ( day > mdays[month] ) {
695                (void) fprintf(fd_sun, "             ");
696                continue;
697             }
698             julianday = greg2jd( year, month, day );
699 
700             yday = (int)(julianday - julianday0);
701 
702             retcode = suntimes(latitude, longitude, timezone, julianday,
703                     sunmode, offset, &rise, &set, NULL, NULL );
704 
705             if ( timemode == TIMEMODE_CIVIL ) {
706                /* Adjust times for Daylight Time */
707                if ( retcode == NORMAL_SUN || retcode == NO_SUNSET ) {
708                   rise += time_adjust(yday, rise, LGL2STD);
709                   if ( rise < 0 || rise > 1439 )
710                      retcode = NO_SUNRISE;
711                }
712                if ( retcode == NORMAL_SUN || retcode == NO_SUNRISE ) {
713                   set += time_adjust(yday, set, LGL2STD);
714                   if ( set < 0 || set > 1439 )
715                      retcode = NO_SUNSET;
716                }
717 
718                /* Mark day of time change */
719                tmark = (yday == 0 ||
720                    time_adjust(yday, 720, LGL2STD) == time_adjust(yday - 1, 720, LGL2STD)) ? ' ' : '*';
721             }
722 
723             scode |= retcode;
724 
725             switch ( retcode ) {
726                case DOWN_ALL_DAY :
727                   (void) fprintf(fd_sun, " %c----- -----", tmark);
728                   break ;
729                case UP_ALL_DAY :
730                   (void) fprintf(fd_sun, " %c***** *****", tmark);
731                   break ;
732                case NO_SUNRISE :
733                   (void) fprintf(fd_sun, " %c      %02d:%02d", tmark, set/60, set % 60);
734                   break ;
735                case NO_SUNSET :
736                   (void) fprintf(fd_sun, " %c%02d:%02d      ", tmark, rise/60, rise%60);
737                   break ;
738                default :
739                   (void) fprintf(fd_sun, " %c%02d:%02d %02d:%02d",
740                                  tmark, rise/60, rise%60, set/60,set%60 );
741             }
742          }
743       }
744    }
745 
746    if ( timemode == TIMEMODE_CIVIL ) {
747       (void) fprintf(fd_sun, "\n\n(*) Denotes time change");
748    }
749    else {
750       (void) fprintf(fd_sun, "\n\nAdd offset for Daylight Time (usually +60 minutes) if and when in effect.");
751    }
752 
753    if ( scode & (DOWN_ALL_DAY | UP_ALL_DAY) ) {
754       (void) fprintf(fd_sun, "\n\n(*****) Sun continuously above horizon");
755       (void) fprintf(fd_sun,  "%*s(-----) Sun continuously below horizon", 4, " ");
756    }
757 
758    (void) fprintf(fd_sun, "\n");
759 
760 
761    return 0;
762 }
763