/*
 * xearth.c
 * kirk johnson
 * july 1993
 *
 * RCS $Id: xearth.c,v 1.5 1993/07/23 05:25:02 tuna Exp $
 *
 * Copyright 1993 by Kirk Lauritz Johnson (see the included file
 * "kljcpyrt.h" for complete copyright information)
 */

#include "xearth.h"
#include "kljcpyrt.h"

#define JUST_OVER_ONE   (1.01)
#define MAP_DATA_SCALE  (30000)

#define XPROJECT(x)     ((proj_xofs+(x))*proj_scale)
#define YPROJECT(y)     ((proj_yofs-(y))*proj_scale)
#define INV_XPROJECT(x) (((x)*inv_proj_scale)-proj_xofs)
#define INV_YPROJECT(y) (proj_yofs-((y)*inv_proj_scale))

#define XingTypeEntry (0)
#define XingTypeExit  (1)

#define ModeX   (0)             /* possible output_mode values */
#define ModePPM (1)
#define ModeGIF (2)

typedef struct
{
  short y;
  short lo_x;
  short hi_x;
  short val;
} ScanBit;

typedef struct
{
  int    type;
  double x, y;
  double angle;
} EdgeXing;

typedef unsigned long Pixel;

void    scan_map();
void    scan_outline();
void    scan_curves();
double *extract_curve();
void    scan_along_curve();
void    handle_xings();
void    scan_arc();
void    scan();
void    get_scanbits();
int     double_comp();
int     scanbit_comp();
int     edgexing_comp();
void    output();
void    render_no_shading();
void    render_with_shading();
void    set_defaults();
void    command_line();
void    usage();
void    fatal();

char    *progname;              /* program name               */
int      output_mode;           /* output mode (X, PPM, ...)  */
double   view_lon;              /* viewing position longitude */
double   view_lat;              /* viewing position latitude  */
int      shading;               /* render with shading?       */
double   sun_lon;               /* sun position longitude     */
double   sun_lat;               /* sun position latitude      */
int      compute_sun_pos;       /* compute sun's position?    */
double   add_lon;               /* additions when nogeostatic */
double   add_lat;               /* additions when nogeostatic */
int      wdth;                  /* image width (pixels)       */
int      hght;                  /* image height (pixels)      */
int      label;                 /* label image                */
int      wait;                  /* wait time between redraw   */
int      day;                   /* day side brightness (%)    */
int      night;                 /* night side brightness (%)  */
int      nofork;		/* don't fork child process?  */
int      geostatic;		/* don't move sun, move earth */
char    *dpyname;               /* X display name             */

double cos_view_lat;
double sin_view_lat;
double cos_view_lon;
double sin_view_lon;

double proj_scale;
double proj_xofs;
double proj_yofs;
double inv_proj_scale;

ExtArr  scanbits;
ExtArr  edgexings;
int     min_y, max_y;
ExtArr *scanbuf;


int main(argc, argv)
     int   argc;
     char *argv[];
{
  command_line(argc, argv);

  if (!geostatic) {
    if (compute_sun_pos)
      sun_position(time(NULL), &sun_lat, &sun_lon);
    view_lon = sun_lon + add_lon;
    view_lat = sun_lat + add_lat;
    }

  scan_map();
  output();

  exit(0);
}


void scan_map()
{
  int i;

  cos_view_lat = cos(view_lat * (M_PI/180));
  sin_view_lat = sin(view_lat * (M_PI/180));
  cos_view_lon = cos(view_lon * (M_PI/180));
  sin_view_lon = sin(view_lon * (M_PI/180));

  if (hght < wdth)
  {
    proj_scale = hght / (2 * JUST_OVER_ONE);
    proj_xofs  = (JUST_OVER_ONE * wdth) / hght;
    proj_yofs  = JUST_OVER_ONE;
  }
  else
  {
    proj_scale = wdth / (2 * JUST_OVER_ONE);
    proj_xofs  = JUST_OVER_ONE;
    proj_yofs  = (JUST_OVER_ONE * hght) / wdth;
  }
  inv_proj_scale = 1 / proj_scale;

  scanbits  = extarr_alloc(sizeof(ScanBit));
  edgexings = extarr_alloc(sizeof(EdgeXing));

  scanbuf = (ExtArr *) malloc(sizeof(ExtArr) * hght);
  assert(scanbuf != NULL);
  for (i=0; i<hght; i++)
    scanbuf[i] = extarr_alloc(sizeof(double));

  scan_outline();
  scan_curves();

  extarr_free(edgexings);
  for (i=0; i<hght; i++)
    extarr_free(scanbuf[i]);
  free(scanbuf);

  qsort(scanbits->body, scanbits->count, sizeof(ScanBit), scanbit_comp);
}


void scan_outline()
{
  min_y = hght;
  max_y = -1;
  scan_arc(1.0, 0.0, 0.0, 1.0, 0.0, 2*M_PI);
  get_scanbits(64);
}


void scan_curves()
{
  int     i;
  int     npts;
  int     val;
  short  *raw;
  double *pos;
  double *prev;
  double *curr;

  raw = map_data;
  while (1)
  {
    npts = *(raw++);
    if (npts == 0) break;
    val = *(raw++);

    pos   = extract_curve(npts, raw);
    prev  = pos + (npts-1)*3;
    curr  = pos;
    min_y = hght;
    max_y = -1;

    for (i=0; i<npts; i++)
    {
      scan_along_curve(prev, curr);
      prev  = curr;
      curr += 3;
    }

    free(pos);
    if (edgexings->count > 0)
      handle_xings();
    if (min_y <= max_y)
      get_scanbits(val);

    raw += 3*npts;
  }
}


double *extract_curve(npts, data)
     int    npts;
     short *data;
{
  int     i;
  double  tmp1;
  double  tmp2;
  double *pos;
  double *rslt;

  rslt = (double *) malloc(sizeof(double) * 3 * npts);
  assert(rslt != NULL);

  pos = rslt;
  for (i=0; i<npts; i++)
  {
    pos[0] = data[0] * (1.0/MAP_DATA_SCALE);
    pos[1] = data[1] * (1.0/MAP_DATA_SCALE);
    pos[2] = data[2] * (1.0/MAP_DATA_SCALE);

    tmp1 = (cos_view_lon * pos[0]) - (sin_view_lon * pos[2]);
    tmp2 = (sin_view_lon * pos[0]) + (cos_view_lon * pos[2]);
    pos[0] = tmp1;
    pos[2] = tmp2;

    tmp1 = (cos_view_lat * pos[1]) - (sin_view_lat * pos[2]);
    tmp2 = (sin_view_lat * pos[1]) + (cos_view_lat * pos[2]);
    pos[1] = tmp1;
    pos[2] = tmp2;

    data += 3;
    pos  += 3;
  }

  return rslt;
}


void scan_along_curve(prev, curr)
     double *prev;
     double *curr;
{
  double    tmp;
  double    extra[3];
  EdgeXing *xing;

  if (prev[2] <= 0)             /* prev not visible */
  {
    if (curr[2] <= 0)
      return;                   /* neither point visible */

    tmp = curr[2] / (curr[2] - prev[2]);
    extra[0] = curr[0] - tmp * (curr[0] - prev[0]);
    extra[1] = curr[1] - tmp * (curr[1] - prev[1]);
    extra[2] = 0;

    tmp = sqrt(extra[0]*extra[0] + extra[1]*extra[1]);
    extra[0] /= tmp;
    extra[1] /= tmp;

    /* extra[] is an edge crossing (entry point) */
    xing = extarr_next(edgexings);
    xing->type  = XingTypeEntry;
    xing->x     = extra[0];
    xing->y     = extra[1];
    xing->angle = atan2(extra[1], extra[0]);

    prev = extra;
  }
  else if (curr[2] <= 0)        /* curr not visible */
  {
    tmp = curr[2] / (curr[2] - prev[2]);
    extra[0] = curr[0] - tmp * (curr[0] - prev[0]);
    extra[1] = curr[1] - tmp * (curr[1] - prev[1]);
    extra[2] = 0;

    tmp = sqrt(extra[0]*extra[0] + extra[1]*extra[1]);
    extra[0] /= tmp;
    extra[1] /= tmp;

    /* extra[] is an edge crossing (exit point) */
    xing = extarr_next(edgexings);
    xing->type  = XingTypeExit;
    xing->x     = extra[0];
    xing->y     = extra[1];
    xing->angle = atan2(extra[1], extra[0]);

    curr = extra;
  }

  scan(XPROJECT(prev[0]), YPROJECT(prev[1]),
       XPROJECT(curr[0]), YPROJECT(curr[1]));
}


void handle_xings()
{
  int       i;
  int       nxings;
  EdgeXing *xings;
  EdgeXing *from;
  EdgeXing *to;

  xings  = edgexings->body;
  nxings = edgexings->count;

  xings  = edgexings->body;
  assert((nxings % 2) == 0);
  qsort(xings, nxings, sizeof(EdgeXing), edgexing_comp);

  if (xings[0].type == XingTypeExit)
  {
    for (i=0; i<nxings; i+=2)
    {
      from = &(xings[i]);
      to   = &(xings[i+1]);
      assert(from->type == XingTypeExit);
      assert(to->type == XingTypeEntry);
      scan_arc(from->x, from->y, from->angle, to->x, to->y, to->angle);
    }
  }
  else
  {
    from = &(xings[nxings-1]);
    to   = &(xings[0]);
    assert(from->type == XingTypeExit);
    assert(to->type == XingTypeEntry);
    assert(from->angle >= to->angle);
    scan_arc(from->x, from->y, from->angle, to->x, to->y, to->angle+2*M_PI);

    for (i=1; i<(nxings-1); i+=2)
    {
      from = &(xings[i]);
      to   = &(xings[i+1]);
      assert(from->type == XingTypeExit);
      assert(to->type == XingTypeEntry);
      scan_arc(from->x, from->y, from->angle, to->x, to->y, to->angle);
    }
  }

  edgexings->count = 0;
}


void scan_arc(x0, y0, a0, x1, y1, a1)
     double x0, y0, a0;
     double x1, y1, a1;
{
  int    i;
  int    lo, hi;
  double angle, step;
  double prev_x, prev_y;
  double curr_x, curr_y;

  assert(a0 < a1);

  step = M_PI / 50;
  lo   = ceil(a0 / step);
  hi   = floor(a1 / step);

  prev_x = XPROJECT(x0);
  prev_y = YPROJECT(y0);

  for (i=lo; i<=hi; i++)
  {
    angle  = i * step;
    curr_x = XPROJECT(cos(angle));
    curr_y = YPROJECT(sin(angle));
    scan(prev_x, prev_y, curr_x, curr_y);

    prev_x = curr_x;
    prev_y = curr_y;
  }

  curr_x = XPROJECT(x1);
  curr_y = YPROJECT(y1);
  scan(prev_x, prev_y, curr_x, curr_y);
}


void scan(x0, y0, x1, y1)
     double x0, y0;
     double x1, y1;
{
  int    i;
  int    lo_y, hi_y;
  double x_value;
  double x_delta;

  if (y0 < y1)
  {
    lo_y = ceil(y0 - 0.5);
    hi_y = floor(y1 - 0.5);

    if (hi_y == (y1 - 0.5))
      hi_y -= 1;
  }
  else
  {
    lo_y = ceil(y1 - 0.5);
    hi_y = floor(y0 - 0.5);

    if (hi_y == (y0 - 0.5))
      hi_y -= 1;
  }

  if (lo_y > hi_y)
    return;                     /* no scan lines crossed */

  if (lo_y < min_y) min_y = lo_y;
  if (hi_y > max_y) max_y = hi_y;

  x_value = x0 - (x0 - x1) * (y0 - (lo_y + 0.5)) / (y0 - y1);
  x_delta = (x0 - x1) / (y0 - y1);

  for (i=lo_y; i<=hi_y; i++)
  {
    *((double *) extarr_next(scanbuf[i])) = x_value;
    x_value += x_delta;
  }
}


void get_scanbits(val)
     int val;
{
  int      i, j;
  int      lo_x, hi_x;
  int      nvals;
  double  *vals;
  ScanBit *scanbit;

  for (i=min_y; i<=max_y; i++)
  {
    vals  = scanbuf[i]->body;
    nvals = scanbuf[i]->count;
    assert((nvals % 2) == 0);
    qsort(vals, nvals, sizeof(double), double_comp);

    for (j=0; j<nvals; j+=2)
    {
      lo_x = ceil(vals[j] - 0.5);
      hi_x = floor(vals[j+1] - 0.5);

      if (lo_x <= hi_x)
      {
        scanbit = extarr_next(scanbits);
        scanbit->y    = i;
        scanbit->lo_x = lo_x;
        scanbit->hi_x = hi_x;
        scanbit->val  = val;
      }
    }

    scanbuf[i]->count = 0;
  }
}


int double_comp(a, b)
     double *a;
     double *b;
{
  double val_a;
  double val_b;
  int    rslt;

  val_a = *a;
  val_b = *b;

  if (val_a < val_b)
    rslt = -1;
  else if (val_a > val_b)
    rslt = 1;
  else
    rslt = 0;

  return rslt;
}


int scanbit_comp(a, b)
     ScanBit *a;
     ScanBit *b;
{
  return (a->y - b->y);
}


int edgexing_comp(a, b)
     EdgeXing *a;
     EdgeXing *b;
{
  double angle_a;
  double angle_b;
  int    rslt;

  angle_a = a->angle;
  angle_b = b->angle;

  if (angle_a < angle_b)
    rslt = -1;
  else if (angle_a > angle_b)
    rslt = 1;
  else
    rslt = 0;

  return rslt;
}


void output()
{
  void (*render)();

  if (shading)
    render = render_with_shading;
  else
    render = render_no_shading;

  switch (output_mode)
  {
  case ModePPM:
    ppm_setup(stdout);
    render(ppm_row);
    break;

  case ModeGIF:
    gif_setup(stdout);
    render(gif_row);
    gif_cleanup();
    break;

  case ModeX:
    if (nofork || (fork() == 0))
    {
      while (1)
      {
        x11_setup();
        render(x11_row);
        x11_cleanup();
        sleep(wait);
      }
    }
    break;

  default:
    assert(0);
  }
}


void render_no_shading(rowfunc)
     void (*rowfunc)();
{
  int      i, j;
  int      cnt;
  ScanBit *scanbit;
  uchar   *scanbuf;
  uchar   *row;
  uchar   *tmp;

  cnt     = scanbits->count;
  scanbit = (ScanBit *) scanbits->body;

  scanbuf = (uchar *) malloc(wdth);
  assert(scanbuf != NULL);
  row = (uchar *) malloc(wdth*3);
  assert(row != NULL);

  for (i=0; i<hght; i++)
  {
    bzero(scanbuf, wdth);
    while ((cnt > 0) && (scanbit->y == i))
    {
      for (j=scanbit->lo_x; j<=scanbit->hi_x; j++)
        scanbuf[j] += scanbit->val;

      scanbit += 1;
      cnt     -= 1;
    }

    tmp = row;
    for (j=0; j<wdth; j++)
    {
      if (scanbuf[j] == 0)
      {
        /* black */
        tmp[0] = 0;
        tmp[1] = 0;
        tmp[2] = 0;
      }
      else if (scanbuf[j] > 64)
      {
        /* green */
        tmp[0] = 0;
        tmp[1] = 255;
        tmp[2] = 0;
      }
      else
      {
        /* blue */
        tmp[0] = 0;
        tmp[1] = 0;
        tmp[2] = 255;
      }
      tmp += 3;
    }

    rowfunc(row);
  }

  free(scanbuf);
  free(row);
}


void render_with_shading(rowfunc)
     void (*rowfunc)();
{
  int      i, j;
  int      cnt;
  int      base_val;
  double   delta_val;
  int      val;
  ScanBit *scanbit;
  char    *scanbuf;
  uchar   *row;
  uchar   *tmp;
  double   tmp1, tmp2;
  double   x, y, z;
  double   sol[3];
  double   scale;
  double  *inv_x;

  cnt     = scanbits->count;
  scanbit = (ScanBit *) scanbits->body;

  scanbuf = (char *) malloc(wdth);
  assert(scanbuf != NULL);
  row = (uchar *) malloc(wdth*3);
  assert(row != NULL);

  /* determine vector pointing at current position of sun,
   * adjusting appropriately for viewing position
   */
  if (compute_sun_pos)
    sun_position(time(NULL), &sun_lat, &sun_lon);

  if (!geostatic) {
    view_lon = sun_lon + add_lon; 
    view_lat = sun_lat + add_lat; 
    scan_map();
    }

  sol[0] = sin(sun_lon * (M_PI/180)) * cos(sun_lat * (M_PI/180));
  sol[1] = sin(sun_lat * (M_PI/180));
  sol[2] = cos(sun_lon * (M_PI/180)) * cos(sun_lat * (M_PI/180));

  tmp1 = (cos_view_lon * sol[0]) - (sin_view_lon * sol[2]);
  tmp2 = (sin_view_lon * sol[0]) + (cos_view_lon * sol[2]);
  sol[0] = tmp1;
  sol[2] = tmp2;

  tmp1 = (cos_view_lat * sol[1]) - (sin_view_lat * sol[2]);
  tmp2 = (sin_view_lat * sol[1]) + (cos_view_lat * sol[2]);
  sol[1] = tmp1;
  sol[2] = tmp2;

  /* precompute INV_XPROJECT() values */
  inv_x = (double *) malloc(sizeof(double) * wdth);
  assert(inv_x != NULL);
  for (i=0; i<wdth; i++)
    inv_x[i] = INV_XPROJECT(i);

  /* precompute shading parameters */
  base_val  = (night * 255.99) / 100;
  delta_val = (day * 255.99) / 100 - base_val;

  for (i=0; i<hght; i++)
  {
    y = INV_YPROJECT(i);

    bzero(scanbuf, wdth);
    while ((cnt > 0) && (scanbit->y == i))
    {
      for (j=scanbit->lo_x; j<=scanbit->hi_x; j++)
        scanbuf[j] += scanbit->val;

      scanbit += 1;
      cnt     -= 1;
    }

    tmp = row;
    for (j=0; j<wdth; j++)
    {
      if (scanbuf[j] == 0)
      {
        /* black */
        tmp[0] = 0;
        tmp[1] = 0;
        tmp[2] = 0;
      }
      else
      {
        x = inv_x[j];
        z = 1 - x*x - y*y;
        z = SQRT(z);

        scale = x*sol[0] + y*sol[1] + z*sol[2];
        if (scale < 0)
        {
          val = base_val;
        }
        else
        {
          val = base_val + scale*delta_val;
          if (val > 255) val = 255;
        }
        assert(val >= 0);

        if (scanbuf[j] > 64)
        {
          /* green */
          tmp[0] = 0;
          tmp[1] = val;
          tmp[2] = 0;
        }
        else
        {
          /* blue */
          tmp[0] = 0;
          tmp[1] = 0;
          tmp[2] = val;
        }
      }
      tmp += 3;
    }

    rowfunc(row);
  }

  free(scanbuf);
  free(row);
  free(inv_x);
}


void set_defaults()
{
  output_mode     = ModeX;
  view_lon        = 0;
  view_lat        = 0;
  shading         = 1;
  compute_sun_pos = 1;
  wdth            = 512;
  hght            = 512;
  label           = 0;
  wait            = 300;
  day             = 100;
  night           = 25;
  nofork          = 0;
  geostatic       = 1;
  dpyname         = NULL;
}


void command_line(argc, argv)
     int   argc;
     char *argv[];
{
  int i;

  progname = argv[0];
  set_defaults();

  for (i=1; i<argc; i++)
  {
    if (strcmp(argv[i], "-ppm") == 0)
    {
      output_mode = ModePPM;
    }
    else if (strcmp(argv[i], "-gif") == 0)
    {
      output_mode = ModeGIF;
    }
    else if (strcmp(argv[i], "-display") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing arg to -display");
      dpyname = argv[i];
    }
  }

  if (output_mode == ModeX)
    x11_resources();

  for (i=1; i<argc; i++)
  {
    if (strcmp(argv[i], "-pos") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing 1st arg to -pos");
      sscanf(argv[i], "%lf", &view_lat);
      if ((view_lat > 90) || (view_lat < -90))
        fatal("viewing latitude must be between -90 and 90");

      i += 1;
      if (i >= argc) usage("missing 2nd arg to -pos");
      sscanf(argv[i], "%lf", &view_lon);
      if ((view_lon > 180) || (view_lon < -180))
        fatal("viewing longitude must be between -180 and 180");
    }
    else if (strcmp(argv[i], "-noshade") == 0)
    {
      shading = 0;
    }
    else if (strcmp(argv[i], "-sunpos") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing 1st arg to -sunpos");
      sscanf(argv[i], "%lf", &sun_lat);
      if ((sun_lat > 90) || (sun_lat < -90))
        fatal("sun latitude must be between -90 and 90");

      i += 1;
      if (i >= argc) usage("missing 2nd arg to -sunpos");
      sscanf(argv[i], "%lf", &sun_lon);
      if ((sun_lon > 180) || (sun_lon < -180))
        fatal("sun longitude must be between -180 and 180");

      compute_sun_pos = 0;
    }
    else if (strcmp(argv[i], "-size") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing 1st arg to -size");
      sscanf(argv[i], "%d", &wdth);

      i += 1;
      if (i >= argc) usage("missing 2nd arg to -size");
      sscanf(argv[i], "%d", &hght);
    }
    else if (strcmp(argv[i], "-label") == 0)
    {
      label = 1;
    }
    else if (strcmp(argv[i], "-wait") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing arg to -wait");
      sscanf(argv[i], "%d", &wait);
      if (wait < 0)
        fatal("arg to -wait must be non-negative");
    }
    else if (strcmp(argv[i], "-day") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing arg to -day");
      sscanf(argv[i], "%d", &day);
      if ((day > 100) || (day < 0))
        fatal("arg to -day must be between 0 and 100");
    }
    else if (strcmp(argv[i], "-night") == 0)
    {
      i += 1;
      if (i >= argc) usage("missing arg to -night");
      sscanf(argv[i], "%d", &night);
      if ((night > 100) || (night < 0))
        fatal("arg to -night must be between 0 and 100");
    }
    else if (strcmp(argv[i], "-nofork") == 0)
    {
      nofork = 1;
    }
    else if (strcmp(argv[i], "-version") == 0)
    {
      printf("xearth version %s\n", VersionString);
      exit(1);
    }
    else if (strcmp(argv[i], "-ppm") == 0)
    {
      /* nothing */
    }
    else if (strcmp(argv[i], "-gif") == 0)
    {
      /* nothing */
    }
    else if (strcmp(argv[i], "-display") == 0)
    {
      i += 1;
      /* nothing */
    }
    else if (strcmp(argv[i], "-nogeostatic") == 0)
    {
      geostatic = 0;

      i += 1;
      if (i >= argc) usage("missing 1st arg to -nogeostatic");
      sscanf(argv[i], "%lf", &add_lat);
      if ((add_lat > 90) || (add_lat < -90))
        fatal("viewing latitude must be between -90 and 90");

      i += 1;
      if (i >= argc) usage("missing 2nd arg to -nogeostatic");
      sscanf(argv[i], "%lf", &add_lon);
      if ((add_lon > 180) || (add_lon < -180))
        fatal("viewing longitude must be between -180 and 180");
    }
    else
    {
      usage(NULL);
    }
  }
}


void usage(msg)
     char *msg;
{
  fprintf(stderr, "\n");
  if (msg != NULL)
    fprintf(stderr, "%s\n", msg);
  fprintf(stderr, "usage: %s", progname);
  fprintf(stderr, " [-pos lat lon] [-noshade] [-sunpos lat lon] [-size w h]\n");
  fprintf(stderr, " [-label] [-wait secs] [-day pct] [-night pct] [-nofork] [-nogeostatic lat lon]\n");
  fprintf(stderr, " [-version] [-ppm] [-gif] [-display dpyname]\n");
  fprintf(stderr, "\n");
  exit(1);
}


void fatal(msg)
     char *msg;
{
  fflush(stdout);
  fprintf(stderr, "\n%s: fatal - %s\n", progname, msg);
  fprintf(stderr, "\n");
  exit(1);
}