/***************************************************************************
**
**   SUN.C   Version 1.0  Michael Schwartz  December 25, 1984
**
** This program uses the USNO algorithm to calculate sunrise and sunset 
** times from standard lattitudes and longitudes.  It has been set up as
** system independently as possible, with the exception being a function
** which gets the current (default) date and day of the year.  It is meant
** to be particularly useful for Jewish persons in out-of-the-way locations
** who might need to know the time of candle_lighting or havdalah.  With
** little difficulty, z'man k'riath shema can be calculated and added to the
** routines.
**
** Inquiries can reach me at PO BOX 24536, Denver, Colorado 80224
** (a place no longer 'out of the way')
**
** Ref: Almanac for Computers, pub. by US NAVAL OBSERVATORY
** Usage:
** sun [-c] [-l[lat]] [-L[long]] [-M] [-m[1-12|a]] [-d[day]]
**     [-t<S|C|N|A>]  [-<12|24>]
**  Latitude in deg,min'sec" N; Longitude in deg,min'sec" W
**
** Accuracy:  The USNO claims accuracy withing 2 minutes except at extreme
**  northern or southern lattitudes.  Comparison to local NWS charts for
**  sunrise and sunset (which are cheap and easy to come by) shows that with
**  the double precision calculations, the charts produced by this program
**  are no more than 1 minute removed from those charts in lattitudes lower
**  than 41 degrees.  Candle lighting times agree with those on popular 
**  calendars also to the 1 minute accuracy.
**
**  The program, unlike its fortran predecessor has a number of 
** important options and defaults.  It is capable of getting today's
** sunrise and sunset at a default location (now Denver, of course),
** producing a calendar-like table of a month or a year, and allowing the
** user to produce reams of data without reinvoking the program.
**
** Permission is granted to distribute and use this program as desired.
**
** Bugs:  Automatic calculation of non-default time zones will not produce
**  correct local times (using the APPROXIMATE_TZ define).  Otherwise, 
**  when lattitudes and longitudes are specified, the result is in Zulu (UT).
**	  No method of calculating DST or aligning years to Fridays/Saturdays
**  or holidays has been provided in the MS-DOS version:  However,
**  as more MS-DOS compilers support the ANSI time standards, the MS-DOS
**  option may not be necessary for your compiler.
**************************************************************************/
#include <stdio.h>

#ifdef UNIX
#define LONG int
#ifdef MSDOS
#undef MSDOS
#endif
#else
#define LONG long
#define MSDOS
#define FUNC_8086 sysint21
#define INT_8086  bdos
/*
** The last two are currently set for the CI-86 compiler.  They are the
** machine interrupt for function calls (function 21) and the general 
** interrupt with return from AL respectively.	Change for other compilers.
** In fact, modern compilers should all be able to use the UNIX def.
*/
#endif

#define TWELVE 1
#define TWENTYFOUR 0

#define SUNSET 0
#define CIVIL_TWILIGHT 1
#define NAUTICAL_TWILIGHT 2
#define ASTRONOMICAL_TWILIGHT 3

int myclock = TWELVE;
/* Chooses a 12 hour clock */
static char *location = "Denver, Colorado (East Side)";
/* Defines the location for the report header */
int longdeg=104, longmin=54, longsec=16, latdeg=39, latmin=41, latsec=45;
/* Lattitude and longitude of east Denver,  Colorado			*/
int TZ = -7;
/* This is the difference between MST and GMT */

#define BAD_ARG -1
#define NOT_IMPLEMENTED -2
#define BAD_MONTH -3
#define BAD_DATE -4

#define NORMAL	   0
#define NO_SUNSET  1
#define NO_SUNRISE 2

static char *months[] = {
	"Nonesuch",
	"January", "February", "March", "April", "May", "June",
	"July", "August", "September", "October", "November", "December"
			};
static int count[] = {
	1,32,60,91,121,152,182,213,244,274,305,335,366
		     };
/*
** Note -- this is for a 365 day year.	The best possible approximation.
** However, if one asks for 2/29, the sheet will still show it.
*/

int cal_opt, month_opt, day_opt, multi_opt, lat_opt, long_opt, ang_opt=0;
int j_opt;

int month_num, month_day, first_day, last_day;

main(argc,argv)
int argc;
char *argv[];
{
  double xsunrise, xsunset;
  int j=0, day, day_adj, status;
/*
** day_adj and status are used for those embarrasing places where the sun
** doesn't rise and set once a day as in the middle latitudes.
** day_adj will be -1 if today's event really happened yesterday, and
**		    1 if today's event won't happen until tomorrow.
** e.g, if today's sunrise occurred at 11:00PM yesterday.
** status will be 0 if NORMAL, NO_SUNRISE or NO_SUNSET if that event doesn't
** happen at all on that day (the range of an acos has been violated)
*/
  char *sunrise, *sunset, *candle_light, *havdalah, *timeadj();
  char *get_date();  /* Puts current date in a string of fixed length */

#define LIGHT -18
#define HAV    42
/* The offsets from sunset for candle lighting and havdalah */

  /* Default is today */
  first_day = last_day = get_current_day();

  /* Handle args */
  while (++j < argc)
  {
    if (argv[j][0] != '-')
    {
	prerror(BAD_ARG,argv[j]);
	prusage();
    }
    else
    {
      switch (argv[j][1])
      {
      case 'c':
	cal_opt++; break;
      case 'M':
	multi_opt++;
	day_opt++; first_day = last_day = get_day("\0");
	break;
      case 'm':
	month_opt++; month_num = get_month(argv[j]+2); break;
      case 'd':
	day_opt++;
	first_day = last_day = get_day(argv[j]+2); break;
      case 'l':
	lat_opt++;
	get_angle("Latitude",&latdeg,&latmin,&latsec,argv[j]+2); break;
      case 'L':
	long_opt++;
	get_angle("Longitude",&longdeg,&longmin,&longsec,argv[j]+2);
#ifdef APPROXIMATE_TZ
	TZ = -((longdeg*60 + longmin) + 450)/900; /* Approximation */
#else
	TZ = 0;	  /* Use GMT */
#endif
	break;
      case 'j':
	j_opt++; break;
      case '1':   /* Assume arg is -12 */
	myclock = TWELVE; break;
      case '2':   /* Assume arg is -24 */
	myclock = TWENTYFOUR; break;
      case 't':   /* Alternate sun angle for sunrise/set calculation */
	if (argv[j][2] == 'S') /* Normal */
	  ang_opt=SUNSET;
	else if (argv[j][2] == 'C') /* Civil Twilight */
	  ang_opt=CIVIL_TWILIGHT;
	else if (argv[j][2] == 'N')
	  ang_opt=NAUTICAL_TWILIGHT;
	else if (argv[j][2] == 'A')
	  ang_opt=ASTRONOMICAL_TWILIGHT;
	else
	  prerror(BAD_ARG,argv[j]);
	break;
      default:
	prerror(BAD_ARG,argv[j]);
      }
    }
  }

  if (cal_opt && (lat_opt || long_opt) && !multi_opt) /* New location */
    get_new_location_name();

  if (cal_opt && !multi_opt) prheader();

  for (day=first_day;day<=last_day;day++)
  {
    status = suntime(&xsunrise,&xsunset,day);
/* Note: an "adjusted hourlength" is (xsunset - xsunrise)/12, unless this is
** less than 0; in that case, add 2.  (24/12)
*/

    if (status & NO_SUNRISE)
       sunrise = "No sunrise today";
    else
       sunrise	    = timeadj("Sunrise at",	   xsunrise,0	 ,&day_adj);
    if (status & NO_SUNSET)
    {
      sunset = "No sunset today";
      if (j_opt)
      {
	candle_light = "Candle lighting at **:**";
	havdalah = "Havdalah at **:**";
      }
    }
    else
    {
      sunset	   = timeadj("Sunset at",	  xsunset ,0	,&day_adj);
      if (j_opt)
      {
	candle_light = timeadj("Candle lighting at",xsunset ,LIGHT,&day_adj);
	havdalah     = timeadj("Havdalah at",	    xsunset ,HAV  ,&day_adj);
      }
    }
    if (cal_opt)
    {
      if (j_opt)
	printf("%-12s  %s  %s  %s  %s\n",get_date(),
	    sunrise, candle_light, sunset, havdalah);
      else
	printf("%-12s   %s  %s\n",get_date(),sunrise,sunset);
    }
    else
    {
      if (j_opt)
	printf("%s %d\n\t%s              %s\n\t%s      %s\n",
	  months[month_num],month_day,sunrise, sunset, candle_light, havdalah);
      else
	printf("%s %2d\n\t%s\t\t%s\n",
	  months[month_num],month_day,sunrise,sunset);
    }
    if (day == count[month_num]-1)  /* We did last day of a month */
    {
      month_num++; month_day=1;
      if (cal_opt && !multi_opt)
      {
	putchar('\f');
	prheader();
      }
    }
    else
    {
      if (cal_opt && !multi_opt)
      {
	if (!(month_day%10))
	{
	  for (j=0;j<80;j++)
	    putchar('-');
	  putchar('\n');
	}
      }
      month_day++;
    }

    if (!(status & NO_SUNRISE))
      free(sunrise);

    if (!(status & NO_SUNSET))
    {
      free(sunset);
      if (j_opt)
      {
	free(candle_light); free(havdalah);
      }
    }
  } 
  if (multi_opt)
	exit(main(argc,argv));
}   

/******************************************************************
C SUBROUTINES
C *****************************************************************
c     this program will calCULATE SUNSET TIMES FOR EAST DENVER.
C	 LAMBDA = 104,54'30" W      Longitude
C	 PHI = 39,42'0"N            Latitude
C    FOR 1978, M=.985600t -3.251 -- M is the Sun's mean anomaly
C	L=M + 1.916 SIN(M) + .020SIN(2M) + 282.565 -- L is Sun's TRUE long.
C    TANa = .91746 TAN L	-- a is Sun's right ascension
C    SIN DELTA = .39782 SIN L	-- Sun's declination
C    COS h = (COSz - SIN DELTA*SIN PHI)/(COS DELTA*COS PHI)
C    T = h + a - 0.065710t -6.620  -- h is local hour angle, a now in hours
C    UT = T + LAMBDA --UT is GMT.  LAMBDA is now in hours (=deg/15)
C
C    WHERE Z is the zenith distance at phenomena in question. Samples:
c    Z=90,50'  REGULAR SUNRISE/SET  (Atmospheric refraction + disk diameter)
C    Z=96, CIVIL TWILIGHT
C    Z=102 NAUTICAL TWILIGHT
C    Z=106 ASTRONOMICAL TWILIGHT
*************************************************************************/
 
#define PI  3.14159265
#define TODEC(a,b,c) ((double)a + (double)b/60.0 + (double)c/3600.0)
#define DEGRAD (PI/180.0)
#define RADDEG (180.0/PI)
#define D 0.91746
#define E 0.39782
#define M(x) (0.985600*x - 3.251)
#define L(x) (x + 1.916*sin(DEGRAD*x) + 0.020*sin(2*DEGRAD*x) + 282.565)
#define ADJ(x) (-0.065710*x - 6.620)

suntime(sunrise,sunset,day)
double *sunrise, *sunset;
int day;
/******************************************************************
**	  --- This is where the REAL work is done ---
** This routine will seem FORTRANy.  An attempt has been made to
** mirror the procedure outlined in the USGS publication; very little
** optimization has been done, so that the connection between the
** publication and the software will be extremely clear.
** Thus, all angles are preserved in degrees (degrees are converted to
** hours by the observation that 360 degrees = 24 hours -- so 15.0
** degrees are 1 hour).
******************************************************************/
{
  double lohr, lambda, phi, xm, xl, zm, zl, a, aa, ahr;
  double sind, sindel, cosd, cosdel, h, hh, hhr, hphr, t, tprime;
  double hprime, tt, ttt, ut, uut, aahr;
  double cosphi, sinphi, cosz;
  double sin(),cos(),acos(),tan(),atan(),sqrt();
  double fabs();
  int retval=NORMAL;
/*********************************************************************
** lambda - longitude in degrees	phi - lattitude in degrees
** [sin|cos][d|del] - functions of pseudo-angle delta
**********************************************************************/
  
  switch(ang_opt)
  {
     case SUNSET:
	cosz = cos(DEGRAD*TODEC(90,50,0)); break;
     case CIVIL_TWILIGHT:
	cosz = cos(DEGRAD*TODEC(96,0,0)); break;
     case NAUTICAL_TWILIGHT:
	cosz = cos(DEGRAD*TODEC(102,0,0)); break;
     case ASTRONOMICAL_TWILIGHT:
	cosz = cos(DEGRAD*TODEC(106,0,0)); break;
  }

  lambda = TODEC(longdeg,longmin,longsec);
  lohr = lambda/15.0;
  phi	 = TODEC(latdeg,latmin,latsec);
  cosphi = cos(DEGRAD*phi);
  sinphi = sin(DEGRAD*phi);

  t = (double)day + (18.0 + lohr)/24.0;
  tprime = (double)day + (6.0 + lohr)/24.0;

  xm = M(t);
  zm = M(tprime);
  xl = L(xm);
  zl = L(zm);
/*   XM,XL ARE BOTH IN DEGREES. */
	
  a  = RADDEG*atan( D*tan(DEGRAD*xl) ); /* Conversion to and from radians */
  aa = RADDEG*atan( D*tan(DEGRAD*zl) );

/******* readjustment of angles a and aa *****************/
  if ( fabs(a+360.0-xl) > 90.0 )  a += 180.0;
  if (a > 360.0) a -= 360.0;

  if ( fabs(aa+360.0-zl) > 90.0 )  aa += 180.0;
  if (aa > 360.0) aa -= 360.0;

  ahr = a/15.0;
  sindel = E*sin(DEGRAD*xl);
  cosdel = sqrt(1.0 - sindel*sindel);  /* cos delta must ALWAYS be >0 */
  h = (cosz - sindel*sinphi)/( cosdel*cosphi );

  aahr= aa/15.0;
  sind	 = E*sin(DEGRAD*zl);
  cosd	 = sqrt(1.0 - sind*sind);
  hh= (cosz - sind  *sinphi)/( cosd  *cosphi );

  if (fabs(h) <= 1.0)
    h = RADDEG*acos(h);
  else
    retval |= NO_SUNSET;

  if (fabs(hh) <= 1.0)
    hh= RADDEG*acos(hh);
  else
    retval |= NO_SUNRISE;

  hhr = h/15.0;
  tt  = hhr + ahr + ADJ(t);
  ut = tt + lohr;

  hprime = 360.0 - hh;	/* Puts sunrise in correct quadrant */
  hphr= hprime/15.0;
  ttt = hphr+ aahr + ADJ(tprime);
  uut = ttt + lohr;

#ifdef DEBUG
  fprintf(stderr,"DEBUG: phi %f  lambda %f\n",phi,lambda);
  fprintf(stderr,"DEBUG: t   %f  tprime %f\n",t,tprime);
  fprintf(stderr,"DEBUG: xl  %f  zl     %f\n",xl,zl);
  fprintf(stderr,"DEBUG: a   %f  aa     %f\n",a,aa);
  fprintf(stderr,"DEBUG: h   %f  hh     %f  hprime  %f\n",h, hh, hprime );
  fprintf(stderr,"DEBUG: hhr %f  hphr   %f\n",hhr,hphr);
  fprintf(stderr,"DEBUG: tt  %f  ttt    %f\n",tt,ttt);
  fprintf(stderr,"DEBUG: ut  %f  uut    %f\n",ut,uut);
#endif

  *sunset  = ut + TZ;
  *sunrise = uut+ TZ;
  return(retval);
}

/**********************************************************************
	SUBROUTINE TIMEADJ(TIME,HR,MIN,MINADJ,DAYADJ)
**********************************************************************/
char *timeadj(mess,time,minadj,dayadj)
char *mess;
double time;
int minadj, *dayadj;
{
  char *calloc();
  static char *str;
  int hour, min, strlen();
  int num = strlen(mess)+9;

  time += (double)minadj/60.0;
  if (time < 0)
  {
    time += 24.0;
    *dayadj -= 1;
  }
  hour = time; /* Type conversion causes truncation */
  min  = (int)(( time - (double)hour )*60.0 + 0.5);
  if (min >= 60)
  {
    hour += 1;
    min  -= 60;
  }
  if (hour > 24)
  {
    hour -= 24;
    *dayadj += 1;
  }
  if (myclock==TWELVE)	     /* Check for 12 hour clock */
    if ( hour > 12)
      hour -= 12;

  str = calloc(num,sizeof(char));
  sprintf(str,"%s %2d:%02d",mess,hour,min);
  return(str);
}

/***************
New functions
***************/
get_current_day()
/* This is a system dependent function which will get the current day of
** the year, and initialize the month-num and month-day
*/
{
  int day;
#ifdef MSDOS
  struct regval { int ax,bx,cx,dx,si,di,ds,es; } regs;
  regs.ax = 0x2a00; /* MS-DOS get date function request */
  FUNC_8086(&regs,&regs); /* Now  DH has month, DL has day */
  month_num = regs.dx>>8;
  month_day = regs.dx&0xff;
  day = count[month_num-1] - 1 + month_day;
#endif
#ifdef UNIX
#include <time.h>

time_t tp;
struct tm	*localtime(),*t;

time(&tp);
t = localtime(&tp);

month_num = t->tm_mon + 1;
month_day = t->tm_mday;
day = t->tm_yday;
if (!lat_opt && !long_opt)
{
  TZ = -timezone / 3600;	/* Timezone is in seconds */
  if (t->tm_isdst) TZ++;
}
#endif

  return(day);
}

#include <ctype.h>
get_day(argin)
char *argin;
/* Prompt if *argin NULL.  Otherwise parse for mm/dd or num */
/* Gets the day of interest, possibly parsed from argin--if not, prompts */
/* Also used as the final exit function if multi_opt			 */
{
  char temp[6], *cp;
  int tnum,j;
  LONG atoi();

  if (*argin == NULL) /* Prompt for input */
  {
    argin = temp;
    get_string(argin,
      "Enter date in mm/dd format: (<CR> to exit): ",5);
    /* This is the exit for multi_opt */
    if (!*argin) exit(0);
    /* If we got a NULL then exit */
  }
  cp = argin;
  tnum = 0;
  while (isdigit(*cp) && *cp)
    tnum = 10*tnum + *cp++ - '0';
  if (*cp)  /* We have a non-digit literal */
  {
    month_num = tnum;
    month_day = atoi(cp+1);
    if (month_num < 1 || month_num > 12)
	prerror(BAD_DATE,argin);
    if (month_day < 1 || month_day > 366)
	prerror(BAD_DATE,argin);
    return(count[tnum-1] - 1 + month_day);
  }
  else
  {
  /* The date is a single integer, and thus the number of a day of the year */
    month_num = 1;
    for (j=0; j<12; j++)
      if (tnum >= count[j]) month_num = j+1;
    month_day = tnum - count[month_num-1] + 1;
    if (tnum < 1 || tnum > 366)
      prerror(BAD_DATE,argin);
    return (tnum);
  }
}

prerror(errnum,mess)
int errnum;
char *mess;
{
  static  char *messages[] =  {
	"***Not used***",
	"Bad argument |%s|\n",
	"Option |%s| not yet implemented\n",
	"Bad argument |%s| to -m[onth] flag: **Fatal**\n",
	"Bad date |%s|.  Heap may be corrupted.\n"
  };

  fprintf(stderr,messages[-errnum],mess);
  return(errnum);
}

get_angle(mess,degrees,minutes,seconds,argin)
char *mess, *argin;
int *degrees, *minutes, *seconds;
/* Gets the angle of interest, either parsed from argin or by prompts */
{
  char temp[12],*cc,*degp,*minp,*secp,*strchr();
  static char *zero = "0";
  LONG atoi();

  if (*argin == NULL) /* Must prompt for angle */
  {
    fprintf(stderr,"Enter %s as deg,min'sec\": ",mess);
    fflush(stderr);
    /* Use stderr so as not to mess up a redirection of output */
    get_string(temp," ",11);
    temp[11]='\0'; /* That's right, I don't trust myself */
    argin = temp;
  }
  if ( (cc = strchr(argin,',') ) > 0)
  {
    degp = argin; *cc = '\0'; minp = cc+1;
  }
  else
  {
    minp = argin;
    degp = zero;
  }

  if ( (cc = strchr(minp,'\'') ) > 0)
  {
    *cc = '\0'; secp = cc + 1;
  }
  else
  {
    secp = minp;
    minp = zero;
  }

  if ( (cc = strchr(secp,'"') ) > 0)
  {
    *cc = '\0';
  }
  else
  {
    secp = zero;
  }

  *degrees = atoi(degp);
  *minutes = atoi(minp);
  *seconds = atoi(secp);
}

static char only_date_space[16];

char *get_date()
/* Converts mm/dd into a string of month day of fixed length */
{
  sprintf(only_date_space,"%s %2d",months[month_num],month_day);
  return(only_date_space);
}

prheader()
/* Used to put page headers for the cal_opt */
{
  printf("\n               ");
  printf("Sunrise and Sunset times for %s\n\n\n",location);
}

get_month(argin)
char *argin;
{
  LONG atoi();

  if (*argin == NULL)
  {
    first_day = count[month_num-1];
    last_day  = count[month_num] -1;
    month_day = 1;
    return(month_num);
  }
  if (*argin == 'a')
  {
     first_day = 1;
     last_day = 366;
     month_num = 1;
     month_day = 1;
     return(month_num);
  }
  month_num = atoi(argin);
  if ((month_num < 1) || (month_num > 12) )
  {
    prerror(BAD_MONTH,argin);
    exit(-1);
  }
  first_day = count[month_num-1];
  last_day  = count[month_num] - 1;
  month_day = 1;
  return(month_num);
}

get_new_location_name()
{
  char *calloc();
  location = calloc(51,sizeof(char));
/* Don't use up this space unless we need it (For multi, use realloc?) */
  get_string(location,
    "Enter new location (Title for report): ",50);
}

prusage()
{
  printf("Usage:\n");
  printf("sun [-d[day]] [-c] [-M] [-l[lat]] [-L[long]]\n");
  printf("    [-m[1-12|a]] [-<12|24>] [-t<S|C|N|A>] [-j]\n\n");
  printf("  -m month\t-c calendar output\t-M prompt multiple days\n");
  printf("  -12 - twelve hour clock.  -24 - twenty four hour clock\n");
  printf("  -t: S regular S-nrise/set. C Civil. N Nautical. A Astronomical\n");
  printf("  -j: Jewish - print candle lighting and havdalah times\n");
  printf("If no values are specified, the d, l and L options will prompt\n");
  printf("Default for date is today, default for m is this month.\n");
  printf("Default is for a 12 hour clock, regular sunrise/set.\n");
  exit(-1);
}

get_string(buffer,mess,maxlen)
char *buffer, *mess;
int maxlen;
{
/* Returns 1 if stuff read, 0 if read had error or EOF */
/* Common routine should work for both MS-DOS and UNIX */
  int n_write; /* In case we wish to check status */
  char *start=buffer;
  char c;

  n_write = write(2,mess,strlen(mess)); /* STDERR, so not mess up output */

  while (read(0,&c,1) == 1) /* Until EOF or error */
  {
    if (c == '\b' || c == 127) /* BS or DEL */
    {
	if (buffer-start) buffer--;
    }
    else if (c == '\t') /* HT */
	*buffer++ = c;
    else if (c == '\n' || c == '\r') /* Don't care about mode for NL */
    {
	*buffer = '\0';
	return(1);
    }
    else if ( (buffer - start) >= maxlen || c < 32)
	continue; /* Ignore CTRLs and flush overage to NL */
    else
	*buffer++ = c;
  }
  return(0); /* Get here only if read fails */
}