/*  File drawmap.c  */

#include "intuition/intuition.h"
#include "graphics/gfxmacros.h"
#include "exec/memory.h"
#include <stdio.h>
#include <fcntl.h>
#include <math.h>
#include <popmenu.h>
#include <drawmap.h>
#include <drawmap-menu.h>

short *map, *map_trig;                 /* workspaces for map  */

struct Screen *s;                      /* pointer to screen   */
struct Window *w;                      /* pointer to Window   */
struct RastPort *rp;                   /* pointer to RastPort */
struct ViewPort *vp;                   /* pointer to ViewPort */
struct AreaInfo mapAreaInfo;
struct TmpRas mapTmpRas;
struct Library *GfxBase;
struct Library *IntuitionBase;
short areaArray[5*NUMPTS];             /* 5 words per point for drawing    */
short area[2*NUMPTS];                  /* storage for each individual area */
unsigned short *arrow, *cross;         /* storage for mouse pointers       */

short rim_in[MAXRIM], rim_out[MAXRIM]; /* storge for rim point numbers */

/* ============================================================= */

main ()

{

   struct IntuiMessage *msg;
   PLANEPTR workspace;
   extern void free_svector();
   extern int readmap();
   extern int HandleEvent();
   extern short *svector();
   int ix;

   if ((GfxBase = OpenLibrary ("graphics.library",0)) == NULL)  {
      printf ("Can't open graphics library\n");
      exit (10);
   }
   if ((IntuitionBase = OpenLibrary ("intuition.library",0)) == NULL)  {
      printf ("Can't open intuition library\n");
      goto end1;
   }
   if ((s = OpenScreen (&mapscreen)) == NULL)  {
      printf ("Can't open screen\n");
      goto end2;
   }
   mapWindow.Screen = s;
   if ((w = OpenWindow (&mapWindow)) == NULL)  {
      printf ("Can't open window\n");
      goto end3;
   }
   vp = &(s->ViewPort);                /* pointer to viewport */
   LoadRGB4 (vp, &mapcolors[0], 4);    /* init. color values  */
                                       /* init. mouse pointers */
   arrow = (UWORD *) AllocMem (arrow_size, MEMF_CHIP);
   for (ix=0; ix<arrow_size/2; ++ix)
      arrow[ix] = arrow_data[ix];
   SetPointer (w, arrow, arrow_size/4-2, 16, arrow_x_offset, arrow_y_offset);
   
   cross = (UWORD *) AllocMem (cross_size, MEMF_CHIP);
   for (ix=0; ix<cross_size/2; ++ix)
      cross[ix] = cross_data[ix];
   
   rp = w->RPort;
   if ((workspace = (PLANEPTR) AllocRaster (WWIDTH,WHEIGHT)) == NULL)  {
      printf ("No space for Temporary Raster\n");
      goto end4;
   }
   InitTmpRas (&mapTmpRas, workspace, RASSIZE(WWIDTH,WHEIGHT));
   rp->TmpRas = &mapTmpRas;            /* link it to the RastPort */
   InitArea (&mapAreaInfo, &areaArray[0], NUMPTS);
   rp->AreaInfo = &mapAreaInfo;        /* link it to the RastPort */

   if ((map = svector (0, MAXVAL-1)) == NULL)  {
      printf ("Can't get space for map\n");
      goto end5;
   }
   if ((ix = readmap (map, 0)) != OK)  { /* read map file */
      printf ("Error reading map file\n");
      goto end6;
   }
   if ((map_trig = svector (0, MAXVAL-1)) == NULL)  {
      printf ("Can't get space for map_trig\n");
      goto end6;
   }
   if ((ix = readmap (map_trig, 1)) != OK)  { /* read map_trig file */
      printf ("Error reading map-trig file\n");
      goto end7;
   }

   SetAPen (rp,ORANGE);                /* land  */
   SetBPen (rp,BLUE);                  /* water */
   SetDrMd (rp,JAM2);

   pi2 = pi/2.;                        /* initialize constants */
   twopi = 2.*pi;
   rad = pi/180.;
   
   view_height = VIEW_HEIGHT;          /* constants for globe view */
   eta = view_height/RE;
   facp = 1. + eta;
   etap = 1./facp;
   
   WaitPort ( w->UserPort );           /* wait for message from */
                                       /*   Intuition           */
   while (1)  {
      msg = (struct IntuiMessage *)GetMsg (w->UserPort);
      if (msg==NULL)
         WaitPort (w->UserPort);
      else
         if (HandleEvent (msg->Class, msg->Code)!=OK)
            break;
   }
   
end7:
   free_svector (map_trig, 0, MAXVAL-1); /* done, so cleanup */
end6:
   free_svector (map, 0, MAXVAL-1);
end5:
   FreeRaster (workspace, WWIDTH, WHEIGHT);
end4:
   FreeMem (cross, cross_size);
   FreeMem (arrow, arrow_size);
   ClearPointer (w);
   CloseWindow (w);
end3:
   CloseScreen (s);
end2:
   CloseLibrary (IntuitionBase);
end1:
   CloseLibrary (GfxBase);
}

/* ============================================================= */

int readmap (ws, iflag)                /* reads map files into memory */

short *ws;
int iflag;

{

   int ibin, disp, numrd;
   
   if (iflag==0)  {
      if ((ibin = open ("map.bin", O_RDONLY)) == -1)
         return (NOT_OK);
   }
   else {
      if ((ibin = open ("map-trig.bin", O_RDONLY)) == -1)
         return (NOT_OK);
   }
   disp = 0;
   while (1)  {
      if ((numrd = read (ibin, (char*) &ws[disp], 4000)) <= 0)
         break;
      disp += numrd / sizeof(short);
   }
   close (ibin);
   if (disp!=MAXVAL)
      return (NOT_OK);
   return (OK);
}

/* ============================================================= */

int HandleEvent (class, code)          /* processes main Intuition events */

long class;
unsigned short code;

{

   struct IntuiMessage *msgf;
   extern long PopChoose();
   extern void fullmap(), globe(), stars(), floodfill(), getcoord();
   extern void box(), grid();
   extern int getbox();
   static char title_FLAT[]     = "Flat Map";
   static char title_MERCATOR[] = "Mercator";
   static char title_GLOBE[]    = "Globe...view from infinitely far away";
   static char title_ORBITAL[]  = "                                        "
                                  "                                        ";
   static char title_zoom[]     = "                                        "
                                  "                                        ";
   static char title_GRID[]     = "Grid";
   static char title_FLOOD[]    = "Flood Fill";
   static char title_BOX[]      = "Box";
   static char title_COLORS[]   = "Colors";
   static char title_CLEARS[]   = "Clear Screen";
   static char title_wait[]     = "Wait for drawing to complete";
   static char press_prompt[]   = "Press left button to select center point";
   static char drag_prompt[]    = "Press and drag left button to select box";
   static char stars_wait[]     = "Wait for stars to appear";
   static char flood_wait[]     = "Flood fill...press left button to select "
                                  "area to fill";
   static char box_error[]      = "Box of zero size not allowed";
   static char grid_error[]     = "Invalid map displayed for Grid option";
   static char fmt_ORBITAL[]    = "Orbital...view from %.2lf kilometers";
   static char fmt_ZOOM_IN[]    = "Zoom In...view from %.2lf kilometers";
   static char fmt_ZOOM_OUT[]   = "Zoom Out...view from %.2lf kilometers";
   long val;
   static long oldval = -1L;
   static int fill = FILL;
   static int color_offset = 0;
   int i, x, y, x0, y0, result;
   double lat[2], lam[2];
   static double lamg, latg;
   
   switch (class) {
   case CLOSEWINDOW:
      return (NOT_OK);
      break;
   case MOUSEBUTTONS:
      switch (code)  {
      case MENUDOWN:
         val = PopChoose (&map_menu, NULL);
         switch (val)  {
         case COLOR_F:                 /* toggle color-fill flag */
            if (map_COLOR_F.Flags & CHECKED)
               fill = FILL;
            else
               fill = NOFILL;
            break;
         case FLAT:                    /* flat map */
            SetWindowTitles (w, title_wait, -1);
            fullmap (map, FLAT, fill);
            SetWindowTitles (w, title_FLAT, -1);
            oldval = val;
            break;
         case MERCATOR:                /* Mercator map */
            SetWindowTitles (w, title_wait, -1);
            fullmap (map, MERCATOR, fill);
            SetWindowTitles (w, title_MERCATOR, -1);
            oldval = val;
            break;
         case GRID:
            if (oldval!=FLAT     && oldval!=MERCATOR &&
                oldval!=GLOBE    && oldval!=ORBITAL  &&
                oldval!=ZOOM_IN  && oldval!=ZOOM_OUT)
               SetWindowTitles (w, grid_error, -1);
            else  {
               SetWindowTitles (w, title_wait, -1);
               grid (oldval, latg, lamg);
               SetWindowTitles (w, title_GRID, -1);
            }
            break;
         case ZOOM_IN:                 /* zoom views */
         case ZOOM_OUT:
            if (val==ZOOM_IN)  {
               view_height /= 2.;
               if (view_height<MIN_HEIGHT)
                  view_height = MIN_HEIGHT;
               sprintf (title_zoom, fmt_ZOOM_IN, view_height);
            }
            else  {
               view_height *= 2.;
               sprintf (title_zoom, fmt_ZOOM_OUT, view_height);
            }
            eta = view_height/RE;
            facp = 1. + eta;
            etap = 1./facp;
         case GLOBE:                 /* globe views */
         case ORBITAL:
            if ( (val!=ZOOM_IN && val!=ZOOM_OUT) ||
             (oldval!=ZOOM_IN && oldval!=ZOOM_OUT &&
              oldval!=GLOBE   && oldval!=ORBITAL) )  {
               if (oldval!=FLAT && oldval!=MERCATOR)  {
                  SetWindowTitles (w, title_wait, -1);
                  fullmap (map, FLAT, FILL);
                  oldval = FLAT;
               }
               SetWindowTitles (w, press_prompt, -1); /* wait for mouse button */
               SetPointer (w, cross, cross_size/4-2, 16, 
                           cross_x_offset, cross_y_offset);
               WaitPort (w->UserPort);
               while (1)  {
                  msgf = (struct IntuiMessage *) GetMsg (w->UserPort);
                  if (msgf==NULL)
                     WaitPort (w->UserPort);
                  else
                     if (msgf->Code==SELECTDOWN)
                        break;
               }
               x = msgf->MouseX;          /* get mouse position */
               y = msgf->MouseY;
               getcoord (x, y, oldval, &latg, &lamg);
               SetPointer (w, arrow, arrow_size/4-2, 16,
                           arrow_x_offset, arrow_y_offset);
            }
            SetWindowTitles (w, stars_wait, -1);
            Move (rp, 0, 0);
            SetBPen (rp, BLACK);       /* black background */
            ClearScreen (rp);
            SetAPen (rp, BLUE);        /* blue globe */
            DrawEllipse (rp, CENTERX, CENTERY, HRADIUS, VRADIUS);
            Flood (rp, 1, CENTERX, CENTERY);
            SetAPen (rp, ORANGE);      /* reset colors */
            SetBPen (rp, BLUE);
            if (val==ORBITAL)  {       /* re-initialize values for */
               view_height = VIEW_HEIGHT; /* orbital view          */
               eta = view_height/RE;
               facp = 1. + eta;
               etap = 1./facp;
               sprintf (title_ORBITAL, fmt_ORBITAL, view_height);
            }
            globe (map, map_trig, latg, lamg, val, fill); /* draw map */
            stars();
            if (val==GLOBE)
               SetWindowTitles (w, title_GLOBE, -1);
            else if (val==ORBITAL)
               SetWindowTitles (w, title_ORBITAL, -1);
            else
               SetWindowTitles (w, title_zoom, -1);
            oldval = val;
            break;
         case FLOOD:                   /* flood fill */
            SetWindowTitles (w, flood_wait, -1);
            SetPointer (w, cross, cross_size/4-2, 16,
                        cross_x_offset, cross_y_offset);
            floodfill ();
            SetPointer (w, arrow, arrow_size/4-2, 16,
                        arrow_x_offset, arrow_y_offset);
            SetWindowTitles (w, title_FLOOD, -1);
            break;
         case BOX:                     /* box */
            if (oldval!=FLAT && oldval!=MERCATOR)  {
               SetWindowTitles (w, title_wait, -1);
               fullmap (map, FLAT, FILL);
               oldval = FLAT;
            }
            SetWindowTitles (w, drag_prompt, -1);
            SetPointer (w, cross, cross_size/4-2, 16,
                        cross_x_offset, cross_y_offset);
            result = getbox (&x0, &y0, &x, &y);
            SetPointer (w, arrow, arrow_size/4-2, 16,
                        arrow_x_offset, arrow_y_offset);
            if (result!=OK)  {
               SetWindowTitles (w, box_error, -1);
               break;
            }
            getcoord (x0, y0, oldval, &lat[0], &lam[0]);
            getcoord (x, y, oldval, &lat[1], &lam[1]);
            SetWindowTitles (w, title_wait, -1);
            Move (rp, 0, 0);
            ClearScreen (rp);
            box (map, lat, lam, fill);
            SetWindowTitles (w, title_BOX, -1);
            oldval = val;
            break;
         case COLORS:                  /* modify color table */
            color_offset += 4;
            if (color_offset>=num_colors)
               color_offset = 0;
            LoadRGB4 (vp, &mapcolors[color_offset], 4);
            SetWindowTitles (w, title_COLORS, -1);
            break;
         case CLEARS:                  /* clear screen */
            SetWindowTitles (w, title_CLEARS, -1);
            Move (rp, 0, 0);
            ClearScreen (rp);
            oldval = val;
            break;
         default:
            if (map_COLOR_F.Flags & CHECKED)
               fill = FILL;
            else
               fill = NOFILL;
            break;
         }
         break;
      default:
         break;
      }
   default:
      break;
   }
   return (OK);
}

/* ============================================================= */

void fullmap (ws, type, fill)          /* draws flat map projections */

short *ws;
int type, fill;

{

   extern void afill(), adraw(), plotpoint();
   int i, j, k, h1, h2, narea;
   double t;
   long pen;
   short np, xs, ys;
   
   Move (rp, 0, 0);                    /* clear screen */
   ClearScreen (rp);
   j = 0;
   np = 0;
   narea = 0;
   pen = ORANGE;
   SetAPen (rp, pen);
   for (i=0; i<MAXVAL; i+=2)  {
      if (ws[i]==0 && ws[i+1]==0)  {   /* check for end of area */
         ++narea;
         pen = ORANGE;
         for (k=0; k<num_lakes; ++k)  {
            if (narea==lakes[k])  {
               pen = BLUE;
               break;
            }
         }
         SetAPen (rp, pen);
         if (np<=1)  {
            xs = CENTERX + area[j-2];
            ys = CENTERY + area[j-1];
            plotpoint (xs, ys);
         }
         else  {
            if (fill==FILL)
               afill (np, pen);
            else
               adraw (np);
         }
         j = 0;
         np = 0;
         area[0] = 0;
         area[1] = 0;
         pen = ORANGE;
         SetAPen (rp, pen);
      }
      else  {                          /* get coords of new point */
         t = ws[i];                    /* y = latitude */
         if (type==FLAT)
            t = (t/100.) * VFACTOR;
         else if (type==MERCATOR)  {
            t = (t/200. + 45.) * rad;
            t = log (tan (t)) * M_VFACTOR;
         }
         if (t<0.)
            t -= 0.5;
         else
            t += 0.5;
         h2 = -t;
         t = ws[i+1];                  /* x = longitude */
         t = (t/100.) * HFACTOR;
         if (t<0.)
            t -= 0.5;
         else
            t += 0.5;
         h1 = t;
         if (np!=0)  {                 /* disallow identical adjacent pts */
            if (h1==area[j-2] && h2==area[j-1])
               continue;
         }
         xs = CENTERX + h1;
         ys = CENTERY + h2;
         plotpoint (xs, ys);
         area[j] = h1;
         area[j+1] = h2;
         j += 2;
         ++np;
      }
   }
   if (np>0)  {                        /* check for last area */
      pen = ORANGE;
      SetAPen (rp, pen);
      if (np==1)  {
         xs = CENTERX + area[j-2];
         ys = CENTERY + area[j-1];
         plotpoint (xs, ys);
      }
      else  {
         if (fill==FILL)
            afill (np, pen);
         else
            adraw (np);
      }
   }
   SetAPen (rp, ORANGE);
}

/* ============================================================= */

void grid (val, lat0, lam0)            /* controls drawing of grids */

long val;
double lat0, lam0;

{

   extern void globe_grid();
   long oldpen, h1;
   double lam, lat, temp;
   static double lat_interval = 20.;
   static double lam_interval = 30.;

   oldpen = rp->FgPen;                 /* save original pen color */
   if (val!=FLAT && val!=MERCATOR)     /* draw grid for globe     */
      globe_grid (val, lat0, lam0);
   else  {                             /* otherwise grid for flat */
                                       /*   or Mercator map       */
      SetAPen (rp, BLACK);             /* set grid color to black */
      for (lam=-180.; lam<=+180.; lam += lam_interval)  {
         h1 = lam*HFACTOR + CENTERX;
         Move (rp, h1, 0);
         Draw (rp, h1, WHEIGHT-1);
      }
      for (lat=80.; lat>=-80.; lat -= lat_interval)  {
         if (val==FLAT)
            h1 = -lat*VFACTOR;
         else
            h1 = -log (tan((lat/2.+45.)*rad)) * M_VFACTOR;
         h1 += CENTERY;
         Move (rp, 0, h1);
         Draw (rp, WWIDTH-1, h1);
      }
      SetAPen (rp, WHITE);             /* draw coordinate axes in white */
      Move (rp, CENTERX, 0);
      Draw (rp, CENTERX, WHEIGHT-1);
      Move (rp, 0, CENTERY);
      Draw (rp, WWIDTH-1, CENTERY);
   }
   SetAPen (rp, oldpen);               /* restore pen color */
}

/* ============================================================= */

void globe_grid (val, lat0, lam0)      /* controls drawing globe grid */

double lat0, lam0;
long val;

{

   extern void globe_grid_plot();
   extern int limit_lam(), limit_lat();
   double lat, lam, c0, s0, c1, s1, c2, s2, rlam, rlat;
   double hp, scale, fac3, dlat, dlam, ddelt;
   double lamp[2], latp[2];
   int k;
   long oldpen;
   char first;
   static double lat_interval = 20.;
   static double lam_interval = 30.;
   static double delta = 5.;
   
   lat0 *= rad;
   c0 = cos (lat0);
   s0 = sin (lat0);
   dlat = lat_interval;
   dlam = lam_interval;
   ddelt = delta;
   if (val==GLOBE)  {                  /* ordinary globe view */
      hp = 0.;
      scale = 1.;
      fac3 = 1.;
   }
   else  {                             /* orbital view */
      hp = etap;
      scale = sqrt (1.-etap*etap);
      fac3 = (facp/(facp-hp)) * scale;
      if (view_height<=1200.)  {
         dlat /= 4.;
         dlam /= 4.;
         ddelt /= 5.;
      }
   }
   SetAPen (rp, BLACK);                /* grid lines in black */
   for (lat=80.; lat>=-80.; lat-=dlat)  {  /* lines of equal latitude */
      rlat = lat*rad;
      if ((k=limit_lam (&lamp[0], rlat, lat0, hp))!=OK)
         continue;                     /* skip if entirely out of view */
      lamp[0] += lam0;
      lamp[1] += lam0;
      k = lamp[0]/ddelt - 1.;          /* express limits as multiple */
      lamp[0] = k*ddelt;               /*   of ddelt                 */
      k = lamp[1]/ddelt + 1.;
      lamp[1] = k*ddelt;
      c1 = cos (rlat);
      s1 = sin (rlat);
      first = TRUE;
      if (lat==0.)
         SetAPen (rp, WHITE);          /* draw equator in white */
      for (lam=lamp[0]; lam<=lamp[1]; lam+=ddelt)  {
         rlam = (lam-lam0)*rad;
         if (rlam<-pi)
            rlam += twopi;
         if (rlam>+pi)
            rlam -= twopi;
         c2 = cos (rlam);
         s2 = sin (rlam);
         globe_grid_plot (rlat, rlam, c0, s0, c1, s1, c2, s2,
                          val, hp, scale, fac3, &first);
      }
      if (lat==0.)
         SetAPen (rp, BLACK);          /* reset pen to black */
   }
   for (lam=-180.; lam<+180.; lam+=dlam)  {  /* meridian circles */
      rlam = (lam-lam0)*rad;
      if (rlam<-pi)
         rlam += twopi;
      if (rlam>+pi)
         rlam -= twopi;
      if ((k=limit_lat (&latp[0], lat0, rlam, hp))!=OK)
         continue;                     /* skip if entirely out of view */
      k = latp[0]/ddelt + 1.;          /* express limits as multiple */
      latp[0] = k*ddelt;               /*   of ddelt                 */
      k = latp[1]/ddelt - 1;
      latp[1] = k*ddelt;
      if (latp[0]>=90.)                /* exclude North polar point */
         latp[0] = 90. - ddelt;
      if (latp[1]<=-90.)               /* exclude South polar point */
         latp[1] = -90. + ddelt;
      c2 = cos (rlam);
      s2 = sin (rlam);
      first = TRUE;
      if (lam==0. || lam==-180.)       /* prime meridian in white */
         SetAPen (rp, WHITE);
      for (lat=latp[0]; lat>=latp[1]; lat-=ddelt)  {
         rlat = lat*rad;
         c1 = cos (rlat);
         s1 = sin (rlat);
         globe_grid_plot (rlat, rlam, c0, s0, c1, s1, c2, s2,
                          val, hp, scale, fac3, &first);
      }
      if (lam==0. || lam==-180.)
         SetAPen (rp, BLACK);          /* reset pen to black */
   }
}

/* ============================================================= */

void globe_grid_plot (lat, lam, c0, s0, c1, s1, c2, s2,
                      val, hp, scale, fac3, first)

double lat, lam, c0, s0, c1, s1, c2, s2;   /* draws globe grids */
long val;
double hp, scale, fac3;
char *first;

{

   extern void plotpoint(), globe_rim_fix();
   char in_view;
   short h1, h1c;                      /* x-dist. (pix) from center */
   short h2, h2c;                      /* y-dist. (pix) from center */
   short xs, ys;
   long x, y;
   double lamc, dlam;                  /* longitude */
   double latc, dlat;                  /* latitude  */
   double zp, facz, h1d, h2d, fac2;
   static char prev_in_view;
   static short h1prev, h2prev;
   static double latprev, lamprev, zpprev;
   
   in_view = FALSE;                    /* get status of current point */
   zp = s1*s0 + c1*c0*c2;
   if (zp>=hp)  {                      /* zp > hp => in view */
      in_view = TRUE;
      h1d = HRADIUS * c1 * s2;
      h2d = -VRADIUS * (s1*c0 - c1*s0*c2);
      if (val!=GLOBE)  {
         fac2 = (facp/(facp-zp))*scale;
         h1d *= fac2;
         h2d *= fac2;
      }
      h1 = h1d;
      h2 = h2d;
   }
   if (*first==TRUE)  {                /* get status of first point */
      *first = FALSE;
      if (in_view==TRUE)  {            /* if first point is in view, */
         x = h1 + CENTERX;             /*   move pen to first point  */
         y = h2 + CENTERY;             /*   and plot it              */
         Move (rp, x, y);
         xs = x;
         ys = y;
         plotpoint (xs, ys);
         h1prev = h1;
         h2prev = h2;
      }
      prev_in_view = in_view;          /* save status of first point */
      latprev = lat;
      lamprev = lam;
      zpprev = zp;
      return;
   }
   if (in_view==TRUE)  {               /* if current point is in view, */
      if (prev_in_view==FALSE)  {      /*   but previous point was not */
         facz = zp / (zpprev-zp);      /*   in view, get rim point by  */
         latc = lat - (latprev-lat)*facz; /* linear interpolation      */
         lamc = lam - (lamprev-lam)*facz;
         dlat = fabs (lat-latprev);
         dlam = fabs (lam-lamprev);
         if ( fabs (latc-latprev)> dlat || /* if rim point not between */
              fabs (latc-lat)    > dlat)   /*   current and previous   */
            latc = (lat+latprev)/2.;       /*   points, use midpoint   */
         if ( fabs (lamc-lamprev)> dlam ||
              fabs (lamc-lam)    > dlam )
            lamc = (lam+lamprev)/2.;
         c1 = cos (latc);
         h1d = HRADIUS * c1 * sin (lamc);
         h2d = -VRADIUS * (sin (latc)*c0 - c1*s0*cos (lamc));
         if (val!=GLOBE)  {
            h1d *= fac3;
            h2d *= fac3;
         }
         h1c = h1d;
         h2c = h2d;
         globe_rim_fix (&h1c, &h2c);
         x = h1c + CENTERX;            /* move to rim point & plot it */
         y = h2c + CENTERY;
         Move (rp, x, y);
         xs = x;
         ys = y;
         plotpoint (xs, ys);
         h1prev = h1c;
         h2prev = h2c;
      }
      if (h1!=h1prev || h2!=h2prev)  {
         x = h1 + CENTERX;             /* draw to current point */
         y = h2 + CENTERY;
         Draw (rp, x, y);
         Move (rp, x, y);
      }
      h1prev = h1;
      h2prev = h2;
   }
   else  {                             /* current point is out of view   */
      if (prev_in_view==TRUE)  {       /* if previous point was in view, */
         facz = zp / (zpprev-zp);      /*   get rim point by linear      */
         latc = lat - (latprev-lat)*facz; /* interpolation               */
         lamc = lam - (lamprev-lam)*facz;
         dlat = fabs (lat-latprev);
         dlam = fabs (lam-lamprev);
         if ( fabs (latc-latprev)> dlat || /* if rim point not between */
              fabs (latc-lat)    > dlat)   /*   current and previous   */
            latc = (lat+latprev)/2.;       /*   points, use midpoint   */
         if ( fabs (lamc-lamprev)> dlam ||
              fabs (lamc-lam)    > dlam )
            lamc = (lam+lamprev)/2.;
         c1 = cos (latc);         
         h1d = HRADIUS * c1 * sin (lamc);
         h2d = -VRADIUS * (c0*sin (latc) - c1*s0 * cos (lamc));
         if (val!=GLOBE)  {
            h1d *= fac3;
            h2d *= fac3;
         }
         h1c = h1d;
         h2c = h2d;
         globe_rim_fix (&h1c, &h2c);
         x = h1c + CENTERX;            /* draw to rim point */
         y = h2c + CENTERY;
         Draw (rp, x, y);
         Move (rp, x, y);
      }
   }
   prev_in_view = in_view;             /* save status of current point */
   latprev = lat;
   lamprev = lam;
   zpprev = zp;
}

/* ============================================================= */

int limit_lam (lam, lat, lat0, etap)   /* computes limits on longitude */
                                       /*   for constant latitude      */
double *lam, lat, lat0, etap;

{
   double alpha;
   
   alpha = (etap-sin(lat)*sin(lat0)) / (cos(lat)*cos(lat0));
   if (alpha<=-1.)  {                  /* negative => lamda covers */
      lam[0] = -180.;                  /*   entire hemisphere      */
      lam[1] = +180.;
      return (OK);
   }
   else if (alpha>=+1.)                /* positive => nothing in view */
      return (NOT_OK);
   else  {                             /* otherwise, compute limits */
      lam[0] = -acos (alpha)/rad;
      lam[1] = -lam[0];
      return (OK);
   }
}

/* ============================================================= */

int limit_lat (lat, lat0, lam, etap)   /* computes limits on latitude */
                                       /*   for constant longitude    */
double *lat, lat0, lam, etap;

{
   double radical, a, b, sum, fac1;
   
   a = sin (lat0);
   b = cos (lat0) * cos (lam);
   sum = a*a + b*b;
   radical = sum - etap*etap;
   if (radical<=0.)                    /* no real solutions */
      return (NOT_OK);
   else  {                             /* two real solutions      */
      radical = sqrt (radical);        /* solve quadrtic equation */
      fac1 = (a*etap + b*radical)/sum;
      lat[0] = asin (fac1)/rad;
      fac1 = (a*etap - b*radical)/sum;
      lat[1] = asin (fac1)/rad;
      if (lat[0]<lat[1])  {            /* put in correct order */
         b = lat[0];
         lat[0] = lat[1];
         lat[1] = b;
      }
      if (a>etap)                      /* check North pole */
         lat[0] = 90.;
      if (a<-etap)                     /* check South pole */
         lat[1] = -90.;
      return (OK);
   }
}

/* ============================================================= */

void getcoord (x, y, val, lat, lam)    /* converts screen coordinates   */
                                       /*   into latitude and longitude */
int x, y;
long val;
double *lat, *lam;

{

   (*lam) = (x - CENTERX) / HFACTOR;
   if (val==FLAT)
      (*lat) = (CENTERY - y) / VFACTOR;
   else
      (*lat) = -90. + 2.*atan(exp((CENTERY-y)/M_VFACTOR))/rad;
}

/* ============================================================= */

void globe (ws, ws_trig, lat0, lam0, val, fill) /* draws globe projections */

short *ws, *ws_trig;
double lat0, lam0;
long val;
int fill;

{

   extern void plotpoint(), globe_fill(), globe_rim_fix();
   char first, prev_in_view, in_view, latzero;
   short h1, h1c, h1prev;              /* x-dist. (pix) from center */
   short h2, h2c, h2prev;              /* y-dist. (pix) from center */
   short np, narea, nrim_in, nrim_out;
   short xs, ys;
   int i, j;
   long x, y, pen, first_color;
   double lam, lamc, lamprev;          /* longitude */
   double lat, latc, latprev;          /* latitude  */
   double c0, s0, c1, s1, c2, zp, zpprev;
   double hp, fac2, fac3, facz, scale;
   double h1d, h2d, dlat, dlam, lam0p;
   static double fac = 1./32767.0;

   latzero = FALSE;
   if (val==GLOBE)  {                  /* ordinary globe view */
      hp = 0.;
      scale = 1.;
      fac3 = 1.;
   }
   else  {                             /* orbital globe view */
      hp = etap;
      scale = sqrt (1.-etap*etap);
      fac3 = (facp/(facp-hp)) * scale;
   }
   if (lat0==0.)  {
      c0 = 1.;
      s0 = 0.;
      if (val==GLOBE)
         latzero = TRUE;               /* equatorial, ordinary globe view */
   }
   else  {
      lat0 *= rad;
      c0 = cos (lat0);
      s0 = sin (lat0);
   }
   first = TRUE;
   pen = ORANGE;
   SetAPen (rp, pen);
   first_color = BLACK;
   j = 0;
   np = 0;
   narea = 0;
   nrim_in = 0;
   nrim_out = 0;
   lam0p = lam0 * rad;                 /* convert lam0 to radians */
   for (i=0; i<MAXVAL; i+=2)  {
      if (ws[i]==0 && ws[i+1]==0)  {   /* if end of area, then     */
         ++narea;                      /*   color-fill and skip to */
         first = TRUE;                 /*    next area             */
         if (np>0 && fill==FILL)
            globe_fill (np, narea, nrim_in, nrim_out, lat0, lam0p,
                        first_color);
         j = 0;
         np = 0;
         nrim_in = 0;
         nrim_out = 0;
         first_color = BLACK;
         continue;
      }
      lat = ws[i];                     /* latitude */
      lat = (lat/100.) * rad;
      lam = ws[i+1];                   /* longitude */
      lam = (lam/100. - lam0) * rad;
      if (lam<-pi)
         lam += twopi;
      if (lam>pi)
         lam -= twopi;
      c1 = ws_trig[i];                 /* cosine of latitude */
      c1 *= fac;
      s1 = ws_trig[i+1];               /* sine of latitude */
      s1 *= fac;
      in_view = FALSE;                 /* get status of current point */
      if (latzero==TRUE)  {            /* equatorial globe view */
         zp = c1*cos (lam);
         if (lam>=-pi2 && lam<=+pi2)  {
            in_view = TRUE;
            h1 = HRADIUS * c1 * sin (lam);
            h2 = -VRADIUS * s1;
         }
      }
      else  {                          /* oblique earth view */
         c2 = cos (lam);
         zp = s1*s0 + c1*c0*c2;
         if (zp>=hp)  {                /* zp > hp => in view */
            in_view = TRUE;
            h1d = HRADIUS * c1 * sin (lam);
            h2d = -VRADIUS * (s1*c0 -c1*s0*c2);
            if (val!=GLOBE)  {
               fac2 = (facp/(facp-zp)) * scale;
               h1d *= fac2;
               h2d *= fac2;
            }
            h1 = h1d;
            h2 = h2d;
         }
      }
      if (first==TRUE)  {              /* get status of first point */
         first = FALSE;
         if (in_view==TRUE)  {         /* if first point is in view, */
            x = h1 + CENTERX;          /*   move pen to first point  */
            y = h2 + CENTERY;          /*   and plot it              */
            Move (rp, x, y);
            first_color = ReadPixel (rp, x, y);
            xs = x;
            ys = y;
            plotpoint (xs, ys);
            h1prev = h1;
            h2prev = h2;
            area[0] = h1;              /* save first point for later use */
            area[1] = h2;
            j = 2;
            np = 1;
         }
         prev_in_view = in_view;       /* save status of first point */
         latprev = lat;
         lamprev = lam;
         zpprev = zp;
         continue;
      }
      if (in_view==TRUE)  {            /* if current point is in view, */
         if (prev_in_view==FALSE)  {   /*   but previous point was not */
            facz = zp / (zpprev-zp);   /*   in view, get rim point by  */
            latc = lat - (latprev-lat)*facz; /* linear interpolation   */
            lamc = lam - (lamprev-lam)*facz;
            dlat = fabs (lat-latprev);
            dlam = fabs (lam-lamprev);
            if ( fabs (latc-latprev)> dlat || /* if rim point not between */
                 fabs (latc-lat)    > dlat)   /*   current and previous   */
               latc = (lat+latprev)/2.;       /*   point, use midpoint    */
            if ( fabs (lamc-lamprev)> dlam ||
                 fabs (lamc-lam)    > dlam )
               lamc = (lam+lamprev)/2.;
            if (latzero==TRUE)  {
               h1c = HRADIUS * cos (latc) * sin (lamc);
               h2c = -VRADIUS * sin (latc);
            }
            else  {
               c1 = cos (latc);
               h1d = HRADIUS * c1 * sin (lamc);
               h2d = -VRADIUS * (sin (latc)*c0 - c1*s0*cos (lamc));
               if (val!=GLOBE)  {
                  h1d *= fac3;
                  h2d *= fac3;
               }
               h1c = h1d;
               h2c = h2d;
            }
            globe_rim_fix (&h1c, &h2c); /* correct the computed rim point */
            x = h1c + CENTERX;         /* move to rim point & plot it */
            y = h2c + CENTERY;
            Move (rp, x, y);
            xs = x;
            ys = y;
            plotpoint (xs, ys);
            h1prev = h1c;
            h2prev = h2c;
            area[j] = h1;              /* save coords of rim point */
            area[j+1] = h2;
            j += 2;
            if (nrim_in<MAXRIM)  {
               rim_in[nrim_in] = np;
               ++nrim_in;
            }
            ++np;
         }
         if (h1!=h1prev || h2!=h2prev)  {
            x = h1 + CENTERX;          /* draw to current point */
            y = h2 + CENTERY;
            Draw (rp, x, y);
            Move (rp, x, y);
            area[j] = h1;              /* save coords of current point */
            area[j+1] = h2;
            j += 2;
            ++np;
         }
         h1prev = h1;
         h2prev = h2;
      }
      else  {                          /* current point is out of view   */
         if (prev_in_view==TRUE)  {    /* if previous point was in view, */
            facz = zp / (zpprev-zp);   /*   get rim point by linear      */
            latc = lat - (latprev-lat)*facz; /* interpolation            */
            lamc = lam - (lamprev-lam)*facz;
            dlat = fabs (lat-latprev);
            dlam = fabs (lam-lamprev);
            if ( fabs (latc-latprev)> dlat || /* if rim point not between */
                 fabs (latc-lat)    > dlat)   /*   current and previous   */
               latc = (lat+latprev)/2.;       /*   point, use midpoint    */
            if ( fabs (lamc-lamprev)> dlam ||
                 fabs (lamc-lam)    > dlam )
               lamc = (lam+lamprev)/2.;
            if (latzero==TRUE)  {
               h1c = HRADIUS * cos (latc) * sin (lamc);
               h2c = -VRADIUS * sin (latc);
            }
            else  {
               c1 = cos (latc);
               h1d = HRADIUS * c1 * sin (lamc);
               h2d = -VRADIUS * (c0*sin(latc) - c1*s0 * cos(lamc));
               if (val!=GLOBE)  {
                  h1d *= fac3;
                  h2d *= fac3;
               }
               h1c = h1d;
               h2c = h2d;
            }
            globe_rim_fix (&h1c, &h2c); /* correct the computed rim point */ 
            x = h1c + CENTERX;         /* draw to rim point */
            y = h2c + CENTERY;
            Draw (rp, x, y);
            Move (rp, x, y);
            area[j] = h1;              /* save coords of rim point */
            area[j+1] = h2;
            j += 2;
            if (nrim_out<MAXRIM)  {
               rim_out[nrim_out] = np;
               ++nrim_out;
            }
            ++np;
         }
      }
      prev_in_view = in_view;          /* save status of current point */
      latprev = lat;
      lamprev = lam;
      zpprev = zp;
   }
}

/* ============================================================= */

void globe_fill (np, narea, nrim_in, nrim_out, lat0, lam0, first_color)
                                       /* fills area on globe */

short np, narea, nrim_in, nrim_out;
double lat0, lam0;
long first_color;

{

   extern void plotpoint ();
   int k;
   long x, y;
   short xs, ys;
   char is_lake;
   static long pen = ORANGE;
   static char area1 = FALSE;

   if (nrim_in!=nrim_out)
      return;
   if (nrim_in==0)  {                  /* fill areas with no rim points */
      is_lake = FALSE;
      for (k=0; k<num_lakes; ++k)      /* check for lake */
         if (narea==lakes[k])  {
            if (first_color==ORANGE)  {
               pen = BLUE;             /* make lake blue if drawn */
               is_lake = TRUE;         /*   into orange           */
               break;
            }
            else                       /* don't color-fill if drawn */
               return;                 /*   into blue               */
         }
      xs = CENTERX + area[0];
      ys = CENTERY + area[1];
      if (np==1)
         plotpoint (xs, ys);
      else if (np>=GLOBE_FILL_MIN)
         afill (np, pen);
      if (narea==1)                    /* check south polar region */
         area1 = TRUE;
      else if (narea==2)  {
         if (area1==TRUE)  {
            area1 = FALSE;
            if (lat0<0.)  {
               x = CENTERX;
               y = CENTERY + VRADIUS * cos (lat0);
               Flood (rp, 1, x, y);
            }
         }
      }
      if (is_lake==TRUE)
         pen = ORANGE;
   }
}

/* ============================================================= */

void globe_rim_fix (h1, h2)            /* find nearest actual rim point */

short *h1, *h2;

{

   short inc1, inc2;
   int i;
   long x, y, color;
   static int itmax = 20;
   
   if ((*h1)>=0)                       /* right half */
      inc1 = +1;
   else                                /* left half */
      inc1 = -1;
   if ((*h2)>=0)                       /* bottom half */
      inc2 = +1;
   else                                /* top half */
      inc2 = -1;
   x = (*h1) + CENTERX + 5*inc1;       /* coords of test pixel */
   y = (*h2) + CENTERY + 5*inc2;
   for (i=0; i<itmax; ++i)  {          /* look for nearest black pixel */
      if ((color = ReadPixel (rp, x, y))==BLACK)
         break;
      x += inc1;
      y += inc2;
   }
   x -= inc1;
   y -= inc2;
   for (i=0; i<itmax; ++i)  {          /* look back for previous blue one */
      if ((color = ReadPixel (rp, x, y))==BLUE || color==ORANGE)
         break;
      x -= inc1;
      y -= inc2;
   }
   (*h1) = x - CENTERX;
   (*h2) = y - CENTERY;
}

/* ============================================================= */

void stars ()                          /* draws stars into black background */

{

   extern double ran();
   static int nmax = 150;
   double t;
   short x, y;
   int i, nstars;
   long oldpen, color;
  
   oldpen = rp->FgPen;
   SetAPen (rp, WHITE);
   nstars = 0;
   for (i=0; ; ++i)  {
      t = ran();
      x = t * WWIDTH + 0.5;
      t = ran();
      y = t * WHEIGHT + 0.5;
      if (x<0)
         x = 0;
      if (x>WWIDTH-1)
         x = WWIDTH-1;
      if (y<0)
         y = 0;
      if (y>WHEIGHT-1)
         y = WHEIGHT-1;
      if ((color = ReadPixel (rp, x, y)) == BLACK)  {
         WritePixel (rp, x, y);
         ++nstars;
         if (nstars>=nmax)
            break;
      }
   }
   SetAPen (rp, oldpen);
}

/* ============================================================= */

void floodfill ()                      /* flood fills an area based */
                                       /*   on mouse position       */
{
   struct IntuiMessage *msgf;
   short x, y;
   long color;
   
   WaitPort (w->UserPort);             /* wait for message */
   while (1)  {
      msgf = (struct IntuiMessage *) GetMsg (w->UserPort);
      if (msgf==NULL)
         WaitPort (w->UserPort);
      else
         if (msgf->Code==SELECTDOWN)
            break;
   }
   x = msgf->MouseX;                   /*  get mouse coordinates  */
   y = msgf->MouseY;                   /*    and flood fill       */
   if ((color = ReadPixel (rp, x, y))==BLUE)
      Flood (rp, 1, x, y);
}

/* ============================================================= */

int getbox (x0, y0, x, y)              /* selects a region to draw to */
                                       /*   larger scale              */
int *x0, *y0, *x, *y;

{

   extern void drawbox();
   struct IntuiMessage *msgf;
   int xd, yd, x1, x2, y1, y2;
   char selectbutton;
   
   selectbutton = FALSE;
   SetDrMd (rp, JAM2 | COMPLEMENT);    /* turn on complement mode and */
   ModifyIDCMP (w, IDCMPFLAGS | MOUSEMOVE); /* enable mouse events    */
   WaitPort (w->UserPort);
   while (1)  {
      msgf = (struct IntuiMessage *) GetMsg (w->UserPort);
      if (msgf==NULL)
         WaitPort (w->UserPort);
      else if (msgf->Code==SELECTDOWN)  { /* if user pressed left button, */
         x1 = msgf->MouseX;               /*   get initial mouse position */
         y1 = msgf->MouseY;
         xd = x1;
         yd = y1;
         selectbutton = TRUE;
      }
      else if (selectbutton==TRUE && msgf->Class==MOUSEMOVE)  {
         x2 = msgf->MouseX;
         y2 = msgf->MouseY;
         drawbox (x1, y1, xd, yd);     /* erase old box */
         xd = x2;
         yd = y2;
         drawbox (x1, y1, xd, yd);     /* draw new box */
      }
      else if (selectbutton==TRUE && msgf->Code==SELECTUP)
         break;
   }
   x2 = msgf->MouseX;                  /* erase current box */
   y2 = msgf->MouseY;
   drawbox (x1, y1, xd, yd);
   SetDrMd (rp, JAM2);                 /* restore original drawing mode */
   ModifyIDCMP (w, IDCMPFLAGS);        /*   and disable mouse events    */
   *x0 = x1;
   *y0 = y1;
   *x = x2;
   *y = y2;
   if (x1==x2 && y1==y2)               /* error if box is of zero size */
      return (NOT_OK);
   if (x1>x2 && y1>y2)  {              /* ensure that the vertices of */
      *x0 = x2;                        /*   the box are in the proper */
      *y0 = y2;                        /*   order                     */
      *x = x1;
      *y = y1;
   }
   else if (x1<x2 && y1>y2)  {
      *x0 = x1;
      *y0 = y2;
      *x = x2;
      *y = y1;
   }
   else if (x1>x2 && y1<y2)  {
      *x0 = x2;
      *y0 = y1;
      *x = x1;
      *y = y2;
   }
   return (OK);
}

/* ============================================================= */

void drawbox (x1, y1, x2, y2)          /* draws a box */

int x1, y1, x2, y2;

{

   Move (rp, x1, y1);
   Draw (rp, x2, y1);
   Move (rp, x2, y1);
   Draw (rp, x2, y2);
   Move (rp, x2, y2);
   Draw (rp, x1, y2);
   Move (rp, x1, y2);
   Draw (rp, x1, y1);
   Move (rp, x1, y1);
}

/* ============================================================= */

void box (ws, latp, lamp, fill)        /* draws areas contained within */
                                       /*   a rectangular region       */
short *ws;
double *latp, *lamp;
int fill;

{

   extern void plotpoint(), drawbox();
   char first, prev_in_view, in_view, is_lake;
   short h1, h1c, h1prev;              /* x-dist. (pix) from center */
   short h2, h2c, h2prev;              /* y-dist. (pix) from center */
   short np, xs, ys, nrim_in, nrim_out;
   int i, j, k, narea;
   long x, y, pen, first_color;
   double lam, lamc, lamprev;          /* longitude */
   double lat, latc, latprev;          /* latitude  */
   double bwidth, bheight, bcx, bcy, xscale, yscale;
   double lat1, lat2, lam1, lam2;

   lat1 = latp[0];                     /* store values for box corners */
   lat2 = latp[1];                     /*   locally                    */
   lam1 = lamp[0];
   lam2 = lamp[1];
   bwidth = lam2 - lam1;               /* box width (degrees)         */
   bheight = lat1 - lat2;              /* box height (degrees)        */
   bcx = (lam1 + lam2)/2.;             /* x-coord of box center (deg) */
   bcy = (lat1 + lat2)/2.;             /* y-coord of box center (deg) */
   xscale = WWIDTH / bwidth;           /* horizontal scale (pix/deg)  */
   yscale = (WHEIGHT-10) / bheight;    /* vertical scale (pix/deg)    */
   j = 0;
   np= 0;
   narea = 0;
   first = TRUE;
   pen = ORANGE;
   SetAPen (rp, pen);
   first_color = BLACK;

   drawbox (0, 0, WWIDTH-1, WHEIGHT-10); /* outline the box */

   for (i=0; i<MAXVAL; i+=2)  {
      if (ws[i]==0 && ws[i+1]==0)  {   /* if end of area, then */
         ++narea;
         if (np>0 && fill==FILL)  {    /* color-fill if requested */
            if (nrim_in!=nrim_out)
               ;
            else if (nrim_in==0)  {
               is_lake = FALSE;
               for (k=0; k<num_lakes; ++k)
                  if (narea==lakes[k])  {
                     is_lake = TRUE;
                     if (first_color==ORANGE)
                        pen = BLUE;
                     break;
                  }
               if (is_lake==FALSE || (is_lake==TRUE && pen==BLUE))  {
                  xs = CENTERX + area[0];
                  ys = CENTERY + area[1];
                  if (np==1)
                     plotpoint (xs, ys);
                  else
                     afill (np, pen);
               }
            }
         }
         j = 0;
         np = 0;
         nrim_in = 0;
         nrim_out = 0;
         pen = ORANGE;
         first_color = BLACK;
         first = TRUE;
         continue;                     /* skip to next point */
      }
      lat = ws[i];                     /* latitude */
      lat = lat/100.;
      lam = ws[i+1];                   /* longitude */
      lam = lam/100.;
      in_view = FALSE;                 /* get status of current point */
      if ( (lat<=lat1 && lat>=lat2)
        && (lam>=lam1 && lam<=lam2))  {
         in_view = TRUE;
         h1 = (lam-bcx) * xscale;
         h2 = - (lat-bcy) * yscale;
      }
      if (first==TRUE)  {              /* check first point */
         first = FALSE;
         if (in_view==TRUE)  {         /* if first point is in view, */
            x = h1 + CENTERX;          /*   move pen to first point  */
            y = h2 + CENTERY;          /*   and plot it              */
            Move (rp, x, y);
            first_color = ReadPixel (rp, x, y);
            xs = x;
            ys = y;
            plotpoint (xs, ys);
            h1prev = h1;
            h2prev = h2;
            area[0] = h1;
            area[1] = h2;
            j = 2;
            np = 1;
         }
         prev_in_view = in_view;       /* save status of first point */
         latprev = lat;
         lamprev = lam;
         continue;
      }
      if (in_view==TRUE)  {
         if (prev_in_view==FALSE)  {   /* if prev. point was not in view, */
            if (lamprev<=lam1)  {      /*   find rim point                */
               lamc = lam1;
               if (latprev<=lat2)      /* lower left */
                  latc = lat2;
               else if (latprev>=lat1) /* upper left */
                  latc = lat1;
               else                    /* left center */
                  latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
            }
            else if (lamprev>=lam2)  {
               lamc = lam2;
               if (latprev>=lat1)      /* upper right */
                  latc = lat1;
               else if (latprev<=lat2) /* lower right */
                  latc = lat2;
               else                    /* right center */
                  latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
            }
            else  {
               if (latprev>=lat1)      /* top center */
                  latc = lat1;
               else                    /* bottom center */
                  latc = lat2;
               lamc = lam + (lamprev-lam)*(latc-lat)/(latprev-lat);
            }
            h1c = (lamc-bcx) * xscale;
            h2c = - (latc-bcy) * yscale;
            x = h1c + CENTERX;         /* move to rim point & plot it */
            y = h2c + CENTERY;
            Move (rp, x, y);
            xs = x;
            ys = y;
            plotpoint (xs, ys);
            h1prev = h1c;
            h2prev = h2c;
            area[j] = h1;
            area[j+1] = h2;
            j += 2;
            if (nrim_in<MAXRIM)  {
               rim_in[nrim_in] = np;
               ++nrim_in;
            }
            ++np;
         }
         if (h1!=h1prev || h2!=h2prev)  {
            x = h1 + CENTERX;          /* draw to current point */
            y = h2 + CENTERY;
            Draw (rp, x, y);
            Move (rp, x, y);
            area[j] = h1;
            area[j+1] = h2;
            j += 2;
            ++np;
         }
         h1prev = h1;
         h2prev = h2;
      }
      else  {                          /* else out of view */
         if (prev_in_view==TRUE)  {    /* if previous point was in view, */
            if (lam<=lam1)  {          /*   find rim point                */
               lamc = lam1;
               if (lat<=lat2)          /* lower left */
                  latc = lat2;
               else if (lat>=lat1)     /* upper left */
                  latc = lat1;
               else                    /* left center */
                  latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
            }
            else if (lam>=lam2)  {
               lamc = lam2;
               if (lat>=lat1)          /* upper right */
                  latc = lat1;
               else if (lat<=lat2)     /* lower right */
                  latc = lat2;
               else                    /* right center */
                  latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
            }
            else  {
               if (lat>=lat1)          /* top center */
                  latc = lat1;
               else                    /* bottom center */
                  latc = lat2;
               lamc = lam + (lamprev-lam)*(latc-lat)/(latprev-lat);
            }
            h1c = (lamc-bcx) * xscale;
            h2c = - (latc-bcy) * yscale;
            x = h1c + CENTERX;         /* draw to rim point */
            y = h2c + CENTERY;
            Draw (rp, x, y);
            Move (rp, x, y);
            area[j] = h1;
            area[j+1] = h2;
            j += 2;
            if (nrim_out<MAXRIM)  {
               rim_out[nrim_out] = np;
               ++nrim_out;
            }
            ++np;
         }
      }
      prev_in_view = in_view;          /* save status of current point */
      latprev = lat;
      lamprev = lam;
   }
}

/* ============================================================= */

void afill (pairs, pen)                /* draws and fills a region  */
                                       /*   using AreaDraw function */
short pairs;                            /* how many pairs of words */
long pen;                              /* pen number to use       */

{

   short i, j;
   long x, y, oldpen;

   oldpen = rp->FgPen;
   SetAPen (rp, pen);
   x = area[0] + CENTERX;
   y = area[1] + CENTERY;
   if (x<0)
      x = 0;
   if (x>WWIDTH-1)
      x = WWIDTH-1;
   if (y<0)
      y = 0;
   if (y>WHEIGHT-1)
      y = WHEIGHT-1;
   AreaMove (rp, x, y);
   j = 2;
   for (i=1; i<pairs; i++)  {
      if (area[j]==0 && area[j+1]==0)
         break;
      x = area[j] + CENTERX;
      y = area[j+1] + CENTERY;
      if (x<0)
         x = 0;
      if (x>WWIDTH-1)
         x = WWIDTH-1;
      if (y<0)
         y = 0;
      if (y>WHEIGHT-1)
         y = WHEIGHT-1;
      AreaDraw (rp, x, y);
      j += 2;
   }   
   AreaEnd (rp);
   SetAPen (rp, oldpen);
}

/* ============================================================= */

void adraw (pairs)                     /* draws area outlines */

short pairs;

{

   short i, j;
   long x, y, oldpen;
   
   oldpen = rp->FgPen;
   SetAPen (rp, ORANGE);
   x = area[0] + CENTERX;
   y = area[1] + CENTERY;
   if (x<0)
      x = 0;
   if (x>WWIDTH-1)
      x = WWIDTH-1;
   if (y<0)
      y = 0;
   if (y>WHEIGHT-1)
      y = WHEIGHT-1;
   Move (rp, x, y);
   j = 2;
   for (i=1; i<pairs; ++i)  {
      if (area[j]==0 && area[j+1]==0)
         break;
      x = area[j] + CENTERX;
      y = area[j+1] + CENTERY;
      if (x<0)
         x = 0;
      if (x>WWIDTH-1)
         x = WWIDTH-1;
      if (y<0)
         y = 0;
      if (y>WHEIGHT-1)
         y = WHEIGHT-1;
      Draw (rp, x, y);
      Move (rp, x, y);
      j += 2;
   }
   SetAPen (rp, oldpen);
}

/* ============================================================= */

void plotpoint (x, y)                  /* plots a single point */

short x, y;

{

   if (x<0)
      x = 0;
   if (x>WWIDTH-1)
      x = WWIDTH-1;
   if (y<0)
      y = 0;
   if (y>WHEIGHT-1)
      y = WHEIGHT-1;
   WritePixel (rp, x, y);
}

/* ============================================================= */

short *svector (nl,nh)                 /* Allocates a short vector */
                                       /*   with range [nl..nh]    */
int nl, nh;

{
   extern char *calloc();
   extern void errors();
   short *v;

   v = (short*) calloc (1, (unsigned) (nh-nl+1)*sizeof(short));
   if (v==NULL)
      return (NULL);
   else
      return v-nl;
}

/* ============================================================= */

void free_svector (v,nl,nh)            /* Frees a short vector   */
                                       /*   allocated by svector */
short *v;
int nl, nh;

{

   free ((char*) (v+nl));
}
