/*
 * FILE:  drawmap.c
 * 
 * Function to generate PostScript instructions for drawing and shading land
 * areas for a given longitude-latitude range.
 * 
 * PARAMETERS:
 * 
 * mapfile   - input file containing longitude-latitude coordinates of land
 * outlines fpps      - pointer to PostScript file for output lon_range -
 * minimum & maximum longitude on x-axis lat_range - minimum & maximum
 * latitude  on y-axis shading - gray level for shading (0.0 = black to 1.0 =
 * white)
 * 
 */
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "map.h"

#define LonToIn(xx)  (( (xx) - origin_d.x) / DegPerInch.x)
#define LatToIn(yy)  (( (yy) - origin_d.y) / DegPerInch.y)

extern XY_TYPE DegPerInch;
static XY_TYPE size_i, origin_d;

#include "mapsub.c"
#if PROTOTYPE_ALLOWED
void drawmap(char *mapfile, 
             FILE * fpps,
             RANGE_TYPE lon_range, RANGE_TYPE lat_range,
             double shading)
#else
void drawmap(mapfile, fpps, lon_range, lat_range, shading)
   char *mapfile;
   FILE *fpps;
   RANGE_TYPE lon_range, lat_range;
   double shading;
#endif
{
   LONLAT_TYPE current, previous, *last_in, first_point;
   FILE *fpmap;
   int path_saved = 0;
   char shade_stroke[20];
   double mid_lon, half_lon_span;
   RANGE_TYPE test_range;

   strcpy(shade_stroke, "stroke ");
   origin_d.x = lon_range.min;
   origin_d.y = lat_range.min;
   size_i.x = (lon_range.max - lon_range.min) / DegPerInch.x;
   size_i.y = (lat_range.max - lat_range.min) / DegPerInch.y;
   half_lon_span = 0.5 * (lon_range.max - lon_range.min);
   /* max. plotting range is 360 deg; so define test_range as 360 deg
      centered on the midpoint of the user's longitude range.  This
      allows us to shift the input map file (0-360) to a range
      appropriate to the user's specifications. */
   mid_lon = lon_range.min + half_lon_span;
   test_range.min = mid_lon - 180.0;
   test_range.max = mid_lon + 180.0;

   if (!strcmp(mapfile, ""))
      return;
   fputs("%---Map---\n", fpps);
   fputs("/shade {gsave setgray eofill grestore} bind def\n", fpps);
   fprintf(fpps,
      "gsave newpath %f %f M %f %f L %f %f L %f %f L clip newpath\n",
      0.0, 0.0,
      size_i.x, 0.0,
      size_i.x, size_i.y,
      0.0, size_i.y);

   if ((fpmap = fopen(mapfile, "r")) == NULL)
   {
      fprintf(stderr, "\n ERROR: Unable to open map file %s.\n", mapfile);
      return;
   }

   current.lon = current.lat = LAND_BOUNDARY;

   while (!feof(fpmap))
   {
      previous.lon = current.lon;
      previous.lat = current.lat;

      if (fscanf(fpmap, " %lf %lf", &current.lon, &current.lat) != 2)
         break;
      /* the lon_range can specify any subset from -180 to 720;
         the map files are assumed to contain longitudes between 0 and 360;
         we need to check whether we need to shift the input longitudes:
         (note: be careful not to change the land & water boundary entries!) */
      if (current.lon >= 0.0 && current.lon < test_range.min)
         current.lon += 360.0;
      if (current.lon > test_range.max)
         current.lon -= 360.0;

      if (in_range(current))
      {
         if (in_range(previous))
         {
            if (fabs(previous.lon - current.lon) < half_lon_span) /* continue existing path */
               fprintf(fpps, "%9.4f %9.4f L\n",
                  LonToIn(current.lon), LatToIn(current.lat));
            else /* must be exiting lon-axis; map continues on opposite side */
            {
               /* extend current path to axis being exited */
               previous.lon = ExtendToAxis(previous.lon);
               fprintf(fpps, "%9.4f %9.4f L\n",
                  LonToIn(previous.lon), LatToIn(previous.lat));
               /* & start new path from reentry axis */
               if (current.lat < -60) /* special case for Antarctica -- need
                                       * axes corners to close the polygon */
               {
                  fprintf(fpps,
                     "%9.4f %9.4f L %9.4f %9.4f L %9.4f %9.4f L %9.4f %9.4f L\n",
                     LonToIn(previous.lon), LatToIn(lat_range.min),
                     LonToIn(ExtendToAxis(current.lon)), LatToIn(lat_range.min),
                     LonToIn(ExtendToAxis(current.lon)), LatToIn(current.lat),
                     LonToIn(current.lon), LatToIn(current.lat));
               }
               else if (path_saved)   /* resuming unfinished business */
               {
                  /* finish up the path on this side */
                  connect(fpps, previous, *pop(), lon_range, lat_range, 'M');
                  fputs(shade_stroke, fpps);
                  fputs("grestore\n", fpps);
                  /* resume path on opposite side */
                  current.lon = ExtendToAxis(current.lon);
                  connect(fpps, *pop(), current, lon_range, lat_range, 'L');
                  path_saved--;
               }
               else  /* save incomplete path & start new one on opposite side */
               {
                  if (push(previous)) return;
                  fprintf(fpps, "%9.4f %9.4f L gsave newpath %9.4f %9.4f M %9.4f %9.4f L\n",
                     LonToIn(previous.lon), LatToIn(previous.lat),
                     LonToIn(ExtendToAxis(current.lon)), LatToIn(current.lat),
                     LonToIn(current.lon), LatToIn(current.lat));
                  current.lon = ExtendToAxis(current.lon);
                  if (push(current)) return;
                  path_saved++;
               }
            }
         }
         else /* !in_range(previous), either new path or resume existing path */
         {
            if (equal(previous.lon, LAND_BOUNDARY) ||
                equal(previous.lon, WATER_BOUNDARY))
            {
               fprintf(fpps, "%9.4f %9.4f M\n", LonToIn(current.lon),
                  LatToIn(current.lat));
               if (push(current))
                  return;       /* remember first point */
            }
            else if ((last_in = pop()) == NULL) /* new path */
            {
               fprintf(fpps, "%9.4f %9.4f M\n",
                  LonToIn(first_point.lon), LatToIn(first_point.lat));
               connect(fpps, first_point, previous, lon_range, lat_range, 'L');
               fprintf(fpps, "%9.4f %9.4f L\n",
                  LonToIn(current.lon), LatToIn(current.lat));
               if (push(first_point)) return;
            }
            else /* reentering range so resume existing path */
            {
               connect(fpps, *last_in, previous, lon_range, lat_range, 'L');
               fprintf(fpps, "%9.4f %9.4f L\n",
                  LonToIn(current.lon), LatToIn(current.lat));
            }
         }
      }
      else /* not in_range(current) -- either boundary or out-of-range */
      {
         if (equal(current.lon, LAND_BOUNDARY))
         {
            if (shading < 1.0)
               sprintf(shade_stroke, "%.3f shade stroke ", shading);
            if (in_range(previous)) /* close preceding land area properly */
            {
               connect(fpps, previous, *pop(), lon_range, lat_range, 'L');
               fputs(shade_stroke, fpps);
            }
            else /* either no path or exited range */
            {
               if ((last_in = pop()) != NULL)   /* exited range */
               {
                  connect(fpps, *last_in, previous, lon_range, lat_range, 'L');
                  /* connect last point to first point; this is usually
                     identical for closed land areas, except those
                     that cross 0 deg longitude */
                  connect(fpps, *last_in, *pop(), lon_range, lat_range, 'M');
                  fputs(shade_stroke, fpps);
               }
            }
         }
         else if (equal(current.lon, WATER_BOUNDARY))
         {
            if (shading < 1.0)
               strcpy(shade_stroke, " 1.0 shade stroke ");
            if (in_range(previous))
            {
               connect(fpps, previous, *pop(), lon_range, lat_range, 'L');
               fprintf(fpps, " %f shade stroke\n", 1.0); /* white-out inland waters */
            }
            else
            {
               if ((last_in = pop()) != NULL)
               {
                  connect(fpps, *last_in, *pop(), lon_range, lat_range, 'L');
                  fprintf(fpps, " %f shade stroke\n", 1.0);
               }
            }
         }
         else /* out-of-range */
         {
            if (in_range(previous)) /* exiting axis range */
            {
               fprintf(fpps, "%9.4f %9.4f L\n",
                  LonToIn(current.lon), LatToIn(current.lat));
               if (push(current))
                  return;   /* remember last point in */
            }
            if (equal(previous.lon, LAND_BOUNDARY) ||
                equal(previous.lon, WATER_BOUNDARY))
            {
               first_point.lon = current.lon;
               first_point.lat = current.lat;
            }
         }
      }
   }				/* end while !feof */

   if (stack_index >= 0)	/* unfinished outlines */
   {
      do
      {
         if (in_range(current))
            connect(fpps, current, (current = *pop()), lon_range, lat_range, 'L');
         else
            connect(fpps, *pop(), (current = *pop()), lon_range, lat_range, 'M');
      } while (stack_index > 0);
      fputs(shade_stroke, fpps);
   }

   fclose(fpmap);
   fputs("grestore\n", fpps);
}