/* 
    * a program that calculates the amount of dark time in some
    * interval, in a very slow and painful way.
    *
    * Usage:    darktime   startdate  enddate
    *
    * where the dates are expressed in form DD/MM/YYYY.
    */

#include <stdio.h>
#include <time.h>
#include <math.h>
#include "apo.h"
#include "pf.h"

#undef DEBUG
#undef DEBUG2

#define SUNALT   -15.0           /* sun must be below this for darkness */
#define MOONALT    0.0           /* moon must be below this for darkness */
#define TIMESTEP 0.006944444     /* fraction of a day -- this is ten minutes */


static double find_dark( /* double jdstart */ );

char *progname = "darktime";


main(argc, argv)
int argc;
char *argv[];
{
   char dateval[100];
   double sum, ndark;
   double ut, jd, jdstart, jdend;
   int day, month, year;

   /* 
    * starting at UT=0 means that we'll always have nights at APO
    * falling completely within a single 24-hour interval.
    */
   ut = 0.0;

   /* parse command-line args, to get starting and ending dates */
   if (argc != 3) {
      fprintf(stderr, "usage: darktime startdate enddate\n");
      exit(1);
   }
   pf_strtoval(argv[1], dateval);
   read_date(dateval, &day, &month, &year);
   get_jd(day, month, year, ut, &jdstart);
   pf_strtoval(argv[2], dateval);
   read_date(dateval, &day, &month, &year);
   get_jd(day, month, year, ut, &jdend);

#ifdef DEBUG
   printf(" start JD %12.1lf    end JD %12.1lf\n", jdstart, jdend);
#endif

   /* for each day, figure out the number of dark hours */
   sum = 0.0;
   for (jd = jdstart; jd < jdend; jd += 1.0) {
      ndark = find_dark(jd);
      printf("%6.1lf %5.2lf \n", jd - 2450000.0, ndark);
      sum += ndark;
   }
   
   exit(0);
}


   /* 
    * given the JD at the start of some day, figure out how many
    * dark hours there are between JD and (JD + 24 hours).
    * return the number of dark hours.
    */


static double
find_dark(jdstart)
double jdstart;
{
   int dark;
   double jd;
   double sum;
   double lst, ra, dec, ha;
   double rah;
   double alt, az, psecz;

   /* this flag is 0 if it's light, 1 if it's dark. */
   dark = 0;
   
   sum = 0.0;
   for (jd = jdstart; jd < jdstart + 1.0; jd += TIMESTEP) {


      /* find the LST corresponding to this JD */
      get_lst(jd, &lst);

      /* figure out the altitude of the sun. */
      sun_radec(jd, &ra, &dec);
      degtohhh(ra, &rah);
      ha = lst - rah;
      to_altaz(ha, dec, &alt, &az, &psecz);
#ifdef DEBUG2
      printf(" JD is %12.3lf   sun alt %6.1lf \n", jd-2450000.0, alt);
#endif

      /* check to see if sun is "up"; if so, quit and go to next timestep */
      if (alt > SUNALT) {

         /* flip the light/dark flag, if necessary */
         if (dark == 1) {
	    dark = 0;
#ifdef DEBUG2
            printf("     darkness ends   at JD %12.3lf\n", jd-2450000.0);
#endif
         }
	 
	 continue;
      }

      /* okay, the sun is down!  Now, let's check the moon. */
      
      /* figure out moon's altitude. */
      moon_radec(jd, &ra, &dec);
      degtohhh(ra, &rah);
      ha = lst - rah;
      to_altaz(ha, dec, &alt, &az, &psecz);
#ifdef DEBUG2
      printf(" JD is %12.3lf  moon alt %6.1lf \n", jd-2450000.0, alt);
#endif

      /* check to see if moon is "up"; if so, quit and go to next timestep */
      if (alt > MOONALT) {

         /* flip the light/dark flag, if necessary */
         if (dark == 1) {
	    dark = 0;
#ifdef DEBUG
            printf("     darkness ends   at JD %12.3lf\n", jd-2450000.0);
#endif
         }
	 
	 continue;
      }

      /* Hey! If we get here, it must be dark.  Add one TIMESTEP to the sum */
      sum += TIMESTEP;
#ifdef DEBUG2
      printf("     add 1 minute, sum now %12.5lf\n", sum);
#endif

      /* flip the light/dark flag, if necessary */
      if (dark == 0) {
	 dark = 1;
#ifdef DEBUG
      printf("     darkness starts at JD %12.3lf\n", jd-2450000.0);
#endif
      }
	 
   }

#ifdef DEBUG
   printf("     we end with dark = %d\n", dark);
#endif

   /* now convert the sum, which is in fractions of a DAY, to hours */
   sum *= 24.0;
   return(sum);
   
}