jevent.cpp

/* jevent.cpp: computes Galilean satellite events

Copyright (C) 2010, Project Pluto

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.    */

/* The following code is an early and somewhat embarrassing effort.
Please do not judge my programming abilities by it.

   It computes the date/times of "classical" Jovian satellite events
(transits,  shadows,  occultations,  eclipses) over a given time
span,  and sorts them out to produce the sort of event list one
sees in _Sky and Telescope_ and similar publications.  It can
also store them in a binary form that is used within Guide,  so
that that program can show Galilean events without having to
do all the math itself.   */

#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include "watdefs.h"
#include "lunar.h"
#include "afuncs.h"
#include "date.h"

#define PI 3.1415926535897932384626433832795028841971693993751058209749445923
#define JUPITER_DIAMETER_IN_KM (88000. * 1.609)
#define AU_PER_JRAD (JUPITER_DIAMETER_IN_KM / AU_IN_KM)
#define EVENT struct event
#define VSOP_CHUNK 2767U

EVENT
   {
   double t;
   unsigned sat, event_type;
   };

/* We don't actually explicitly use these four macros :
#define OCCULTA  0x00
#define TRANSIT  0x01
#define ECLIPSE  0x02
#define SHADOW   0x03
*/
#define IN_FRONT  0x01
#define FROM_SUN 0x02

#define EVENT_START   0x04
#define EVENT_UNSEEN  0x08

const double speeds[4] = {203.488955432, 101.374724550, 50.317609110, 21.571071314};
static char FAR *vsop_data;
static int quiet = 0;

static void show_event( FILE *ofile, EVENT *e)
{
   const char *event_str[4] = {"Occ", "Tra", "Ecl", "Sha"};
   char buff[80];
   const double ut_jd = e->t - td_minus_ut( e->t) / seconds_per_day;

   if( e->event_type & EVENT_UNSEEN)
      return;
   full_ctime( buff, ut_jd, FULL_CTIME_FORMAT_HH_MM);
   fprintf( ofile, "sat %u: %s %s: %s\n", e->sat, event_str[e->event_type & 3],
              (e->event_type & EVENT_START) ? "start" : "end  ", buff);
}

static unsigned find_events( unsigned sat_no, double t1, double t2, int viewpoint, EVENT *e)
{
   double t, lon_j, lat_j, rad_j, lon_e, lat_e, rad_e;
   double loc[9], tloc[18], *tptr, step, delta_lat, prev_delta_lat = 0.;
   double tc, lon_s, lat_s, rad_s, prev_delta = 0., delta;
   int i;
   unsigned rval = 0;

   t = t1;
   step = 10. / speeds[sat_no - 1];
   while( t < t2)
      {
      tc = (t - 2451545.) / 36525.;     /* re-cvt to julian centuries */

      lon_j = calc_vsop_loc( vsop_data, 5, 0, tc, 0.);
      lat_j = calc_vsop_loc( vsop_data, 5, 1, tc, 0.);
      rad_j = calc_vsop_loc( vsop_data, 5, 2, tc, 0.);

      lon_e = calc_vsop_loc( vsop_data, 3, 0, tc, 0.);
      lat_e = calc_vsop_loc( vsop_data, 3, 1, tc, 0.);
      rad_e = calc_vsop_loc( vsop_data, 3, 2, tc, 0.);

      loc[0] = rad_j * cos( lat_j) * cos( lon_j);
      loc[1] = rad_j * cos( lat_j) * sin( lon_j);
      loc[2] = rad_j * sin( lat_j);

      loc[6] = rad_e * cos( lat_e) * cos( lon_e);
      loc[7] = rad_e * cos( lat_e) * sin( lon_e);
      loc[8] = rad_e * sin( lat_e);

      calc_jsat_loc( t, tloc, 1 << (sat_no - 1), 0L);
      tptr = tloc + (sat_no - 1) * 3;
      loc[3] = loc[0] + tptr[0] * AU_PER_JRAD;
      loc[4] = loc[1] + tptr[1] * AU_PER_JRAD;
      loc[5] = loc[2] + tptr[2] * AU_PER_JRAD;

      if( viewpoint)
         {
         for( i = 0; i < 9; i++)
            loc[i] -= loc[6 + i % 3];
         lon_j = atan2( loc[1], loc[0]);
         lat_j = atan( loc[2] / sqrt( loc[0] * loc[0] + loc[1] * loc[1]));
         }
      rad_s = sqrt( loc[3] * loc[3] + loc[4] * loc[4] + loc[5] * loc[5]);
      lon_s = atan2( loc[4], loc[3]);
      while( lon_s - lon_j > PI)
         lon_s -= PI + PI;
      while( lon_s - lon_j <-PI)
         lon_s += PI + PI;
      lat_s = asin( loc[5] / rad_s);
      delta = (lon_s - lon_j) * rad_s / AU_PER_JRAD;
      delta_lat = (lat_s - lat_j) * rad_s / AU_PER_JRAD;
      delta_lat *= 1.071374;     /* stretch for jup's oblateness */
      if( delta * prev_delta < 0.)     /* zero crossed */
         {
         double t_crossing, diff, a, b, c, dx, dy, dist = 0.;

         dx = prev_delta - delta;
         dy = prev_delta_lat - delta_lat;
         a = dx * dx + dy * dy;              /* quadratic for intercept */
         b = 2. * (dx * delta + dy * delta_lat);
         c = delta_lat * delta_lat + delta * delta - 1.;
         diff = b * b - 4. * a * c;
         t_crossing = t + b * step / (2. * a);
         t = t_crossing + (180. / speeds[sat_no - 1]) * .9;
         if( diff > 0.)       /* if real solution, ie, sat doesn't miss */
            {
            diff = sqrt( diff) * step / (2. * a);
            diff = fabs( diff);
            for( i = 0; i < 3; i++)
               dist += (loc[i + 6] - loc[i]) * (loc[i + 6] - loc[i]);
            t_crossing += sqrt( dist) / AU_PER_DAY;
            e[0].t = t_crossing - diff;
            e[0].sat = sat_no;
            e[0].event_type = (prev_delta > 0.);
            if( !viewpoint)
               e[0].event_type |= FROM_SUN;
            e[1].t = t_crossing + diff;
            e[1].sat = sat_no;
            e[1].event_type = e[0].event_type;
            e[0].event_type |= EVENT_START;
            if( !quiet)
               {
               show_event( stdout, e);
               show_event( stdout, e + 1);
               }
            e += 2;
            rval += 2;
            }
         delta = 0.;
         }
      prev_delta = delta;
      prev_delta_lat = delta_lat;
      t += step;
      }
   return( rval);
}

/* see 'shellsor.cpp' in the 'find_orb' repository
for comments on this ShellSort implementation */

static void _sort_events( EVENT *e, const size_t n_events)
{
   size_t gap = 250104703;

   while( gap < n_events)
      gap = gap * 8 / 3 + 1;
   while( (gap = gap * 3 / 8) != 0)
      {
      size_t i, j;

      for( i = 0; i < gap; i++)
         for( j = i; j + gap < n_events; j += gap)
            if( e[j].t > e[j + gap].t)
               {
               EVENT temp;

               memcpy( &temp, e + j + gap, sizeof( EVENT));
               memcpy( e + j + gap, e + j, sizeof( EVENT));
               memcpy( e + j, &temp, sizeof( EVENT));
               if( j >= gap)
                  j -= gap + gap;
               }
      }
}

/* Sort events by date,  then look for the start of an eclipse or
occultation.  Anything from then to the end of that eclipse or
occultation is "unseen" (occurs while the satellite is behind
Jupiter or in eclipse). */

static void _mark_hidden_events( EVENT *e, const int n_events)
{
   int i;

   _sort_events( e, n_events);
   for( i = 0; i < n_events - 1; i++)
      if( !(e[i].event_type & IN_FRONT) && (e[i].event_type & EVENT_START))
         {
         int j = i + 1;

         while( j < n_events && ((e[j].event_type ^ e[i].event_type) & 7) != EVENT_START)
            {
            e[j].event_type |= EVENT_UNSEEN;
            j++;
            }
         }
}

int main( int argc, char **argv)
{
   unsigned i, julian = 0;
   unsigned n_days = 30, sat_no = 15;
   unsigned max_events;
   unsigned n_events = 0, n_sun, n_earth;
   double t1, t2;
   long jd;
   EVENT *e;
   FILE *ofile = NULL, *data_file = NULL;
   FILE *vsop_file;
   char *vsop_tbuff;

   vsop_file = fopen( "vsop.bin", "rb");
   if( !vsop_file)
      {
      printf( "Couldn't open vsop.bin");
      return( -3);
      }
   vsop_tbuff = (char *)malloc( VSOP_CHUNK);
   vsop_data = (char FAR *)FMALLOC( VSOP_CHUNK * 22U);
   for( i = 0; i < 22; i++)
      {
      if( !fread( vsop_tbuff, VSOP_CHUNK, 1, vsop_file))
         {
         printf( "Couldn't read VSOP data\n");
         free( vsop_tbuff);
         free( vsop_data);
         return( -1);
         }
      memcpy( vsop_data + (unsigned)i * VSOP_CHUNK, vsop_tbuff, VSOP_CHUNK);
      }
   fclose( vsop_file);
   free( vsop_tbuff);

   if( argc < 4)
      {
      printf( "JEVENT calculates Jovian satellite events (transits, eclipses,\n");
      printf( "shadows,  occultations) as seen from Earth.  It requires, as a\n");
      printf( "minimum,  a day,  month,  and year.  Given,  say,  the command\n");
      printf( "\nJEVENT 18 4 1993\n\n");
      printf( "JEVENT will calculate all events from 18 Apr 1993 for the next\n");
      printf( "thirty days.  You can add on the following parameters:\n\n");
      printf( "   -j        Use Julian calendar\n");
      printf( "   -d(#)     Calculate for (#) days instead of 30\n");
      printf( "   -f(name)  Put results in ASCII file (name) as well as on screen\n");
      return( -2);
      }
   for( i = 0; i < (unsigned)argc; i++)
      if( argv[i][0] == '-')
         switch( argv[i][1])
            {
            case 'j': case 'J':
               julian = 1;
               break;
            case 'q': case 'Q':
               quiet = 1;
               break;
            case 's': case 'S':
               sat_no = (unsigned)atoi( argv[i] + 2);
               break;
            case 'd': case 'D':
               n_days = (unsigned)atoi( argv[i] + 2);
               break;
            case 'f': case 'F':
               ofile = fopen( argv[i] + 2, "wb");
               break;
            case 'r': case 'R':
               data_file = fopen( argv[i] + 2, "ab");
               printf( "Appending data to %s\n", argv[i] + 2);
               break;
            default:
               break;
            }
                     /* We seem to average about eight events/day,  before */
                     /* sorting and removing 'hidden' events.  But let's   */
                     /* allow a little margin :                            */
   max_events = n_days * 10 + 50;
   e = (EVENT *)calloc( max_events, sizeof( EVENT));
   if( !e)
      return( -1);
   jd = dmy_to_day( 0, atoi( argv[2]), atol( argv[3]), (int)julian);
   t1 = (double)jd - .5 + atof( argv[1]);
   t2 = t1 + (double)n_days;
   printf( "JD %f to %f\n", t1, t2);
   for( i = 0; i < 4; i++)
      if( sat_no & (1 << i))
         {
         printf( "Sat %u from sun\n", i + 1);
         n_sun = find_events( i + 1, t1 - 1., t2 + 1., 0, e + n_events);
         printf( "Sat %u from earth\n", i + 1);
         n_earth = find_events( i + 1, t1 - 1., t2 + 1., 1, e + n_events + n_sun);
         printf( "Finding hidden events\n");
         _mark_hidden_events( e + n_events, n_sun + n_earth);
         n_events += n_earth + n_sun;
         }
   printf( "Sorting %u events\n", n_events);
   _sort_events( e, n_events);

   if( !quiet)
      {
      printf( "Final results:\n");
      for( i = 0; i < n_events; i++)
         if( e[i].t > t1 && e[i].t < t2)
            show_event( stdout, e + i);
      }

   if( ofile)
      {
      for( i = 0; i < n_events; i++)
         if( e[i].t > t1 && e[i].t < t2)
            show_event( ofile, e + i);
      fclose( ofile);
      }

   if( data_file)
      {
      for( i = 0; i < n_events; i++)
         if( e[i].t > t1 && e[i].t < t2)
            if( !( e[i].event_type & EVENT_UNSEEN))
               {
               char buff[5];
               int32_t tval;

                        /* store minutes from 2000.0 */

               tval = (int32_t)( (e[i].t - 2451545.0) * 1440.);
               memcpy( buff, &tval, sizeof( int32_t));
               buff[4] = (char)( (e[i].event_type & 7) | (e[i].sat << 3));
               fwrite( buff, 1, 5, data_file);
               }
      fclose( data_file);
      }
   return( 0);
}