#include <math.h>
#include <stdio.h>
#include <string.h>

int find_millennium_atlas_page( const double ra, const double dec);
int find_atlas_pages( const double ra, const double dec,
                                         int *sky_atlas, int *urano);
int get_new_uranometria_page( const double ra_hrs, const double dec_deg);

/* Given an RA in decimal hours and a dec in decimal degrees,  the
   following routine returns the number of the Millennium Star Atlas
   page that best shows that location. */

int find_millennium_atlas_page( const double ra, const double dec)
{
   int rval;

   if( dec >= 87.)       /* polar cap pages */
      rval = ( ra < 4. || ra > 16. ? 2 : 1);
   else if( dec <= -87.)       /* polar cap pages */
      rval = ( ra < 4. || ra > 16. ? 516 : 515);
   else
      {
      int gore = (int)( ra / 8.), zone = (int)(( 93. - dec) / 6.);
      double remains = ceil( ra / 8.) * 8. - ra;
      const char per_zone[31] = { 2, 4, 8, 10, 12, 14, 16, 20, 20, 22, 22,
                                     24, 24, 24, 24, 24, 24, 24, 24, 24,
                                 22, 22, 20, 20, 16, 14, 12, 10, 8, 4, 2 };

      rval = (int)( remains * (double)per_zone[zone] / 8.) + 1 + gore * 516;
      while( zone--)
         rval += per_zone[zone];
      }
   return( rval);
}

#define N_DIVIDES 17

/* Given an RA in decimal hours and a dec in decimal degrees,  the
   following routine returns the numbers of the _Sky Atlas 2000_ and
   (original) _Uranometria_ pages that best show that location.  See
   below for code for the "new" Uranometria. */

int find_atlas_pages( const double ra, const double dec,
                                         int *sky_atlas, int *urano)
{
   static const short dec_limits[N_DIVIDES + 1] = { -900, -845, -725, -610,
        -500, -390, -280, -170, -55, 55, 170,
         280, 390, 500, 610, 725, 845, 900 };
   static const char n_divides[N_DIVIDES] = { 2, 12, 20, 24, 30, 36, 45, 45,
        45, 45, 45, 36, 30, 24, 20, 12, 2 };
   int start_value = 472;
   int divide;
   double angle;

   if( fabs( dec) < 18.5)           /* between -18.5 and 18.5 */
      {
      *sky_atlas = 9 + (int)( ra / 3. + 5. / 6.);
      if( *sky_atlas == 9)
         *sky_atlas = 17;
      }
   else if( fabs( dec) < 52.)       /* between 18.5 and 52, N and S */
      {
      *sky_atlas = 4 + (int)( ra / 4.);
      if( dec < 0.)
         *sky_atlas += 14;
      }
   else                             /* above 52,  N and S */
      {
      *sky_atlas = 1 + (int)( ra / 8.);
      if( dec < 0.)
         *sky_atlas += 23;
      }

   for( divide = 0; (double)dec_limits[divide + 1] < dec * 10.; divide++)
      start_value -= (int)n_divides[divide + 1];
   angle = ra * (int)n_divides[divide] / 24.;
   if( (int)n_divides[divide] >= 20)
      angle += .5;
   else if( (int)n_divides[divide] == 12)
      angle += 5. / 12.;
   *urano = (int)angle % (int)n_divides[divide] + start_value;
   if( *urano >= 472)         /* charts 472 and 473 are "flipped" */
      *urano = (472 + 473) - *urano;
   return( 0);
}

/* 17 November 2003:  first cut at code for the "new Uranometria"
page layout,  which is quite a bit different from the "old" layout. */

/*
zone page dec   RA
 0     1 +89.9   0 00  One polar zone
 1     2 +79.0   0 00  Six four-hour zones
 2     8 +68.0   0 00  Ten 2.4 hour zones
 3    18 +57.0   0 00  Twelve two-hour zones
 4    30 +46.0   0 00  Fifteen 1.6-hour zones
 5    45 +35.0   0 00  Eighteen 1 1/3 hour zones
 6    63 +23.5   0 00  Eighteen 1 1/3 hour zones
 7    81 +11.5   0 00  Twenty 1.2 hour zones
 8   101   0.0   0 00  Twenty 1.2 hour zones
 9   121 -11.5   0 00  Twenty 1.2 hour zones
10   141 -23.5   0 00  Eighteen 1 1/3 hour zones
11   159 -35.0   0 00  Eighteen 1 1/3 hour zones
12   177 -46.0   0 00  Fifteen 1.6-hour zones
13   192 -57.0   0 00  Twelve two-hour zones
14   204 -68.0   0 00  Ten 2.4 hour zones
15   214 -79.0   0 00  Six four-hour zones
16   220 -89.9   0 00  One polar zone   */

/* Given an RA in decimal hours and a dec in decimal degrees,  the
   following routine returns the number of the new Uranometria 2000
   page that best shows that location.  See 'new_uran.dat' for a
   list of page centers.  Basically,  the sky is divided into sixteen
   zones in declination,  with each zone divided into one to twenty
   pages (fewer pages closer to the poles,  with the zones at the
   celestial equator and just above and below rendered as twenty pages.)
   Within a zone,  pages start at RA=0,  then work their way _west_
   around the sky (hence the "24-ra" in this code.)
 */

int get_new_uranometria_page( const double ra_hrs, const double dec_deg)
{
   int zone_no = (int)( (90. - dec_deg) * 16. / 180. + .5);
   const char zone_splits[17] = { 1, 6, 10, 12, 15, 18, 18, 20, 20,
               20, 18, 18, 15, 12, 10, 6, 1 };
   int rval = (int)( (24. - ra_hrs) * (double)zone_splits[zone_no] / 24. + .5);

   rval %= zone_splits[zone_no];
   while( zone_no--)
      rval += zone_splits[zone_no];
   return( rval + 1);
}

#define PI 3.141592652589793

/* Given a _lunar_ latitude/longitude,  this figures out the best page to
   look for in Antonin Rukl's _Atlas of the Moon_.  The result is stored
   as a character string rather than a page number,  since the libration
   zones are stored on Roman-numeral pages.  */

int get_rukl_page( const double lon, const double lat, char *buff)
{
   double x = cos( lat) * cos( lon);
   double y = cos( lat) * sin( lon);
   double z = sin( lat);
   int rval = -1, ix = (int)( y * 5.5 + 5.5), iy = (int)( 4. - 4. * z);
   const char page_starts[8] = { -1, 7, 17, 28, 39, 50, 60, 68 };

   strcpy( buff, "Rukl ");
   if( x > 0.)
      {
      if( iy <= 1 || iy >= 6)
         if( !ix)
            ix = 1;
         else if( ix == 10)
            ix = 9;
      if( !iy || iy == 7)
         if( ix == 1)
            ix = 2;
         else if( ix == 9)
            ix = 8;
      rval = ix + page_starts[iy];
      sprintf( buff + 5, "%d", rval);
      }
               /* Allow a basically eight-degree libration zone.  This */
               /* isn't _perfect_,  but it's "not bad" either. */
   if( x < PI * 8. / 180. && x > -PI * 8. / 180.)
      {
      int zone_no = (int)( atan2( z, y) * 4. / PI + 4.);
      const char *libration_zone_pages[8] = { "VII", "VI", "V", "IV",
                        "III", "II", "I", "VIII" };

      if( rval > -1)
         strcat( buff, "/");
      strcat( buff, libration_zone_pages[zone_no]);
      rval = zone_no + 256;
      }
   strcat( buff, "   ");
   if( rval == -1)
      *buff = '\0';
   return( rval);
}

/* A little test routine I used to check out everything except the Rukl code.
*/

#ifdef TEST_ROUTINE
#include <stdlib.h>

void main( int argc, char **argv)
{
   int urano, sky_atlas;
   double ra = atof( argv[1]);
   double dec = atof( argv[2]);

   find_atlas_pages( ra, dec, &sky_atlas, &urano);
   printf( "Sky atlas: %d\n", sky_atlas);
   printf( "Uranometria: %d\n", urano);
   printf( "Millennium: %d\n", find_millennium_atlas_page( ra, dec));
}
#endif
