home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_disks / 200-299 / ff229.lzh / DrawMap / drawmap.c < prev    next >
C/C++ Source or Header  |  1989-07-20  |  54KB  |  1,695 lines
  1. /*  File drawmap.c  */
  2.  
  3. #include "intuition/intuition.h"
  4. #include "graphics/gfxmacros.h"
  5. #include "exec/memory.h"
  6. #include <stdio.h>
  7. #include <fcntl.h>
  8. #include <math.h>
  9. #include <popmenu.h>
  10. #include <drawmap.h>
  11. #include <drawmap-menu.h>
  12.  
  13. short *map, *map_trig;                 /* workspaces for map  */
  14.  
  15. struct Screen *s;                      /* pointer to screen   */
  16. struct Window *w;                      /* pointer to Window   */
  17. struct RastPort *rp;                   /* pointer to RastPort */
  18. struct ViewPort *vp;                   /* pointer to ViewPort */
  19. struct AreaInfo mapAreaInfo;
  20. struct TmpRas mapTmpRas;
  21. struct Library *GfxBase;
  22. struct Library *IntuitionBase;
  23. short areaArray[5*NUMPTS];             /* 5 words per point for drawing    */
  24. short area[2*NUMPTS];                  /* storage for each individual area */
  25. unsigned short *arrow, *cross;         /* storage for mouse pointers       */
  26.  
  27. short rim_in[MAXRIM], rim_out[MAXRIM]; /* storge for rim point numbers */
  28.  
  29. /* ============================================================= */
  30.  
  31. main ()
  32.  
  33. {
  34.  
  35.    struct IntuiMessage *msg;
  36.    PLANEPTR workspace;
  37.    extern void free_svector();
  38.    extern int readmap();
  39.    extern int HandleEvent();
  40.    extern short *svector();
  41.    int ix;
  42.  
  43.    if ((GfxBase = OpenLibrary ("graphics.library",0)) == NULL)  {
  44.       printf ("Can't open graphics library\n");
  45.       exit (10);
  46.    }
  47.    if ((IntuitionBase = OpenLibrary ("intuition.library",0)) == NULL)  {
  48.       printf ("Can't open intuition library\n");
  49.       goto end1;
  50.    }
  51.    if ((s = OpenScreen (&mapscreen)) == NULL)  {
  52.       printf ("Can't open screen\n");
  53.       goto end2;
  54.    }
  55.    mapWindow.Screen = s;
  56.    if ((w = OpenWindow (&mapWindow)) == NULL)  {
  57.       printf ("Can't open window\n");
  58.       goto end3;
  59.    }
  60.    vp = &(s->ViewPort);                /* pointer to viewport */
  61.    LoadRGB4 (vp, &mapcolors[0], 4);    /* init. color values  */
  62.                                        /* init. mouse pointers */
  63.    arrow = (UWORD *) AllocMem (arrow_size, MEMF_CHIP);
  64.    for (ix=0; ix<arrow_size/2; ++ix)
  65.       arrow[ix] = arrow_data[ix];
  66.    SetPointer (w, arrow, arrow_size/4-2, 16, arrow_x_offset, arrow_y_offset);
  67.    
  68.    cross = (UWORD *) AllocMem (cross_size, MEMF_CHIP);
  69.    for (ix=0; ix<cross_size/2; ++ix)
  70.       cross[ix] = cross_data[ix];
  71.    
  72.    rp = w->RPort;
  73.    if ((workspace = (PLANEPTR) AllocRaster (WWIDTH,WHEIGHT)) == NULL)  {
  74.       printf ("No space for Temporary Raster\n");
  75.       goto end4;
  76.    }
  77.    InitTmpRas (&mapTmpRas, workspace, RASSIZE(WWIDTH,WHEIGHT));
  78.    rp->TmpRas = &mapTmpRas;            /* link it to the RastPort */
  79.    InitArea (&mapAreaInfo, &areaArray[0], NUMPTS);
  80.    rp->AreaInfo = &mapAreaInfo;        /* link it to the RastPort */
  81.  
  82.    if ((map = svector (0, MAXVAL-1)) == NULL)  {
  83.       printf ("Can't get space for map\n");
  84.       goto end5;
  85.    }
  86.    if ((ix = readmap (map, 0)) != OK)  { /* read map file */
  87.       printf ("Error reading map file\n");
  88.       goto end6;
  89.    }
  90.    if ((map_trig = svector (0, MAXVAL-1)) == NULL)  {
  91.       printf ("Can't get space for map_trig\n");
  92.       goto end6;
  93.    }
  94.    if ((ix = readmap (map_trig, 1)) != OK)  { /* read map_trig file */
  95.       printf ("Error reading map-trig file\n");
  96.       goto end7;
  97.    }
  98.  
  99.    SetAPen (rp,ORANGE);                /* land  */
  100.    SetBPen (rp,BLUE);                  /* water */
  101.    SetDrMd (rp,JAM2);
  102.  
  103.    pi2 = pi/2.;                        /* initialize constants */
  104.    twopi = 2.*pi;
  105.    rad = pi/180.;
  106.    
  107.    view_height = VIEW_HEIGHT;          /* constants for globe view */
  108.    eta = view_height/RE;
  109.    facp = 1. + eta;
  110.    etap = 1./facp;
  111.    
  112.    WaitPort ( w->UserPort );           /* wait for message from */
  113.                                        /*   Intuition           */
  114.    while (1)  {
  115.       msg = (struct IntuiMessage *)GetMsg (w->UserPort);
  116.       if (msg==NULL)
  117.          WaitPort (w->UserPort);
  118.       else
  119.          if (HandleEvent (msg->Class, msg->Code)!=OK)
  120.             break;
  121.    }
  122.    
  123. end7:
  124.    free_svector (map_trig, 0, MAXVAL-1); /* done, so cleanup */
  125. end6:
  126.    free_svector (map, 0, MAXVAL-1);
  127. end5:
  128.    FreeRaster (workspace, WWIDTH, WHEIGHT);
  129. end4:
  130.    FreeMem (cross, cross_size);
  131.    FreeMem (arrow, arrow_size);
  132.    ClearPointer (w);
  133.    CloseWindow (w);
  134. end3:
  135.    CloseScreen (s);
  136. end2:
  137.    CloseLibrary (IntuitionBase);
  138. end1:
  139.    CloseLibrary (GfxBase);
  140. }
  141.  
  142. /* ============================================================= */
  143.  
  144. int readmap (ws, iflag)                /* reads map files into memory */
  145.  
  146. short *ws;
  147. int iflag;
  148.  
  149. {
  150.  
  151.    int ibin, disp, numrd;
  152.    
  153.    if (iflag==0)  {
  154.       if ((ibin = open ("map.bin", O_RDONLY)) == -1)
  155.          return (NOT_OK);
  156.    }
  157.    else {
  158.       if ((ibin = open ("map-trig.bin", O_RDONLY)) == -1)
  159.          return (NOT_OK);
  160.    }
  161.    disp = 0;
  162.    while (1)  {
  163.       if ((numrd = read (ibin, (char*) &ws[disp], 4000)) <= 0)
  164.          break;
  165.       disp += numrd / sizeof(short);
  166.    }
  167.    close (ibin);
  168.    if (disp!=MAXVAL)
  169.       return (NOT_OK);
  170.    return (OK);
  171. }
  172.  
  173. /* ============================================================= */
  174.  
  175. int HandleEvent (class, code)          /* processes main Intuition events */
  176.  
  177. long class;
  178. unsigned short code;
  179.  
  180. {
  181.  
  182.    struct IntuiMessage *msgf;
  183.    extern long PopChoose();
  184.    extern void fullmap(), globe(), stars(), floodfill(), getcoord();
  185.    extern void box(), grid();
  186.    extern int getbox();
  187.    static char title_FLAT[]     = "Flat Map";
  188.    static char title_MERCATOR[] = "Mercator";
  189.    static char title_GLOBE[]    = "Globe...view from infinitely far away";
  190.    static char title_ORBITAL[]  = "                                        "
  191.                                   "                                        ";
  192.    static char title_zoom[]     = "                                        "
  193.                                   "                                        ";
  194.    static char title_GRID[]     = "Grid";
  195.    static char title_FLOOD[]    = "Flood Fill";
  196.    static char title_BOX[]      = "Box";
  197.    static char title_COLORS[]   = "Colors";
  198.    static char title_CLEARS[]   = "Clear Screen";
  199.    static char title_wait[]     = "Wait for drawing to complete";
  200.    static char press_prompt[]   = "Press left button to select center point";
  201.    static char drag_prompt[]    = "Press and drag left button to select box";
  202.    static char stars_wait[]     = "Wait for stars to appear";
  203.    static char flood_wait[]     = "Flood fill...press left button to select "
  204.                                   "area to fill";
  205.    static char box_error[]      = "Box of zero size not allowed";
  206.    static char grid_error[]     = "Invalid map displayed for Grid option";
  207.    static char fmt_ORBITAL[]    = "Orbital...view from %.2lf kilometers";
  208.    static char fmt_ZOOM_IN[]    = "Zoom In...view from %.2lf kilometers";
  209.    static char fmt_ZOOM_OUT[]   = "Zoom Out...view from %.2lf kilometers";
  210.    long val;
  211.    static long oldval = -1L;
  212.    static int fill = FILL;
  213.    static int color_offset = 0;
  214.    int i, x, y, x0, y0, result;
  215.    double lat[2], lam[2];
  216.    static double lamg, latg;
  217.    
  218.    switch (class) {
  219.    case CLOSEWINDOW:
  220.       return (NOT_OK);
  221.       break;
  222.    case MOUSEBUTTONS:
  223.       switch (code)  {
  224.       case MENUDOWN:
  225.          val = PopChoose (&map_menu, NULL);
  226.          switch (val)  {
  227.          case COLOR_F:                 /* toggle color-fill flag */
  228.             if (map_COLOR_F.Flags & CHECKED)
  229.                fill = FILL;
  230.             else
  231.                fill = NOFILL;
  232.             break;
  233.          case FLAT:                    /* flat map */
  234.             SetWindowTitles (w, title_wait, -1);
  235.             fullmap (map, FLAT, fill);
  236.             SetWindowTitles (w, title_FLAT, -1);
  237.             oldval = val;
  238.             break;
  239.          case MERCATOR:                /* Mercator map */
  240.             SetWindowTitles (w, title_wait, -1);
  241.             fullmap (map, MERCATOR, fill);
  242.             SetWindowTitles (w, title_MERCATOR, -1);
  243.             oldval = val;
  244.             break;
  245.          case GRID:
  246.             if (oldval!=FLAT     && oldval!=MERCATOR &&
  247.                 oldval!=GLOBE    && oldval!=ORBITAL  &&
  248.                 oldval!=ZOOM_IN  && oldval!=ZOOM_OUT)
  249.                SetWindowTitles (w, grid_error, -1);
  250.             else  {
  251.                SetWindowTitles (w, title_wait, -1);
  252.                grid (oldval, latg, lamg);
  253.                SetWindowTitles (w, title_GRID, -1);
  254.             }
  255.             break;
  256.          case ZOOM_IN:                 /* zoom views */
  257.          case ZOOM_OUT:
  258.             if (val==ZOOM_IN)  {
  259.                view_height /= 2.;
  260.                if (view_height<MIN_HEIGHT)
  261.                   view_height = MIN_HEIGHT;
  262.                sprintf (title_zoom, fmt_ZOOM_IN, view_height);
  263.             }
  264.             else  {
  265.                view_height *= 2.;
  266.                sprintf (title_zoom, fmt_ZOOM_OUT, view_height);
  267.             }
  268.             eta = view_height/RE;
  269.             facp = 1. + eta;
  270.             etap = 1./facp;
  271.          case GLOBE:                 /* globe views */
  272.          case ORBITAL:
  273.             if ( (val!=ZOOM_IN && val!=ZOOM_OUT) ||
  274.              (oldval!=ZOOM_IN && oldval!=ZOOM_OUT &&
  275.               oldval!=GLOBE   && oldval!=ORBITAL) )  {
  276.                if (oldval!=FLAT && oldval!=MERCATOR)  {
  277.                   SetWindowTitles (w, title_wait, -1);
  278.                   fullmap (map, FLAT, FILL);
  279.                   oldval = FLAT;
  280.                }
  281.                SetWindowTitles (w, press_prompt, -1); /* wait for mouse button */
  282.                SetPointer (w, cross, cross_size/4-2, 16, 
  283.                            cross_x_offset, cross_y_offset);
  284.                WaitPort (w->UserPort);
  285.                while (1)  {
  286.                   msgf = (struct IntuiMessage *) GetMsg (w->UserPort);
  287.                   if (msgf==NULL)
  288.                      WaitPort (w->UserPort);
  289.                   else
  290.                      if (msgf->Code==SELECTDOWN)
  291.                         break;
  292.                }
  293.                x = msgf->MouseX;          /* get mouse position */
  294.                y = msgf->MouseY;
  295.                getcoord (x, y, oldval, &latg, &lamg);
  296.                SetPointer (w, arrow, arrow_size/4-2, 16,
  297.                            arrow_x_offset, arrow_y_offset);
  298.             }
  299.             SetWindowTitles (w, stars_wait, -1);
  300.             Move (rp, 0, 0);
  301.             SetBPen (rp, BLACK);       /* black background */
  302.             ClearScreen (rp);
  303.             SetAPen (rp, BLUE);        /* blue globe */
  304.             DrawEllipse (rp, CENTERX, CENTERY, HRADIUS, VRADIUS);
  305.             Flood (rp, 1, CENTERX, CENTERY);
  306.             SetAPen (rp, ORANGE);      /* reset colors */
  307.             SetBPen (rp, BLUE);
  308.             if (val==ORBITAL)  {       /* re-initialize values for */
  309.                view_height = VIEW_HEIGHT; /* orbital view          */
  310.                eta = view_height/RE;
  311.                facp = 1. + eta;
  312.                etap = 1./facp;
  313.                sprintf (title_ORBITAL, fmt_ORBITAL, view_height);
  314.             }
  315.             globe (map, map_trig, latg, lamg, val, fill); /* draw map */
  316.             stars();
  317.             if (val==GLOBE)
  318.                SetWindowTitles (w, title_GLOBE, -1);
  319.             else if (val==ORBITAL)
  320.                SetWindowTitles (w, title_ORBITAL, -1);
  321.             else
  322.                SetWindowTitles (w, title_zoom, -1);
  323.             oldval = val;
  324.             break;
  325.          case FLOOD:                   /* flood fill */
  326.             SetWindowTitles (w, flood_wait, -1);
  327.             SetPointer (w, cross, cross_size/4-2, 16,
  328.                         cross_x_offset, cross_y_offset);
  329.             floodfill ();
  330.             SetPointer (w, arrow, arrow_size/4-2, 16,
  331.                         arrow_x_offset, arrow_y_offset);
  332.             SetWindowTitles (w, title_FLOOD, -1);
  333.             break;
  334.          case BOX:                     /* box */
  335.             if (oldval!=FLAT && oldval!=MERCATOR)  {
  336.                SetWindowTitles (w, title_wait, -1);
  337.                fullmap (map, FLAT, FILL);
  338.                oldval = FLAT;
  339.             }
  340.             SetWindowTitles (w, drag_prompt, -1);
  341.             SetPointer (w, cross, cross_size/4-2, 16,
  342.                         cross_x_offset, cross_y_offset);
  343.             result = getbox (&x0, &y0, &x, &y);
  344.             SetPointer (w, arrow, arrow_size/4-2, 16,
  345.                         arrow_x_offset, arrow_y_offset);
  346.             if (result!=OK)  {
  347.                SetWindowTitles (w, box_error, -1);
  348.                break;
  349.             }
  350.             getcoord (x0, y0, oldval, &lat[0], &lam[0]);
  351.             getcoord (x, y, oldval, &lat[1], &lam[1]);
  352.             SetWindowTitles (w, title_wait, -1);
  353.             Move (rp, 0, 0);
  354.             ClearScreen (rp);
  355.             box (map, lat, lam, fill);
  356.             SetWindowTitles (w, title_BOX, -1);
  357.             oldval = val;
  358.             break;
  359.          case COLORS:                  /* modify color table */
  360.             color_offset += 4;
  361.             if (color_offset>=num_colors)
  362.                color_offset = 0;
  363.             LoadRGB4 (vp, &mapcolors[color_offset], 4);
  364.             SetWindowTitles (w, title_COLORS, -1);
  365.             break;
  366.          case CLEARS:                  /* clear screen */
  367.             SetWindowTitles (w, title_CLEARS, -1);
  368.             Move (rp, 0, 0);
  369.             ClearScreen (rp);
  370.             oldval = val;
  371.             break;
  372.          default:
  373.             if (map_COLOR_F.Flags & CHECKED)
  374.                fill = FILL;
  375.             else
  376.                fill = NOFILL;
  377.             break;
  378.          }
  379.          break;
  380.       default:
  381.          break;
  382.       }
  383.    default:
  384.       break;
  385.    }
  386.    return (OK);
  387. }
  388.  
  389. /* ============================================================= */
  390.  
  391. void fullmap (ws, type, fill)          /* draws flat map projections */
  392.  
  393. short *ws;
  394. int type, fill;
  395.  
  396. {
  397.  
  398.    extern void afill(), adraw(), plotpoint();
  399.    int i, j, k, h1, h2, narea;
  400.    double t;
  401.    long pen;
  402.    short np, xs, ys;
  403.    
  404.    Move (rp, 0, 0);                    /* clear screen */
  405.    ClearScreen (rp);
  406.    j = 0;
  407.    np = 0;
  408.    narea = 0;
  409.    pen = ORANGE;
  410.    SetAPen (rp, pen);
  411.    for (i=0; i<MAXVAL; i+=2)  {
  412.       if (ws[i]==0 && ws[i+1]==0)  {   /* check for end of area */
  413.          ++narea;
  414.          pen = ORANGE;
  415.          for (k=0; k<num_lakes; ++k)  {
  416.             if (narea==lakes[k])  {
  417.                pen = BLUE;
  418.                break;
  419.             }
  420.          }
  421.          SetAPen (rp, pen);
  422.          if (np<=1)  {
  423.             xs = CENTERX + area[j-2];
  424.             ys = CENTERY + area[j-1];
  425.             plotpoint (xs, ys);
  426.          }
  427.          else  {
  428.             if (fill==FILL)
  429.                afill (np, pen);
  430.             else
  431.                adraw (np);
  432.          }
  433.          j = 0;
  434.          np = 0;
  435.          area[0] = 0;
  436.          area[1] = 0;
  437.          pen = ORANGE;
  438.          SetAPen (rp, pen);
  439.       }
  440.       else  {                          /* get coords of new point */
  441.          t = ws[i];                    /* y = latitude */
  442.          if (type==FLAT)
  443.             t = (t/100.) * VFACTOR;
  444.          else if (type==MERCATOR)  {
  445.             t = (t/200. + 45.) * rad;
  446.             t = log (tan (t)) * M_VFACTOR;
  447.          }
  448.          if (t<0.)
  449.             t -= 0.5;
  450.          else
  451.             t += 0.5;
  452.          h2 = -t;
  453.          t = ws[i+1];                  /* x = longitude */
  454.          t = (t/100.) * HFACTOR;
  455.          if (t<0.)
  456.             t -= 0.5;
  457.          else
  458.             t += 0.5;
  459.          h1 = t;
  460.          if (np!=0)  {                 /* disallow identical adjacent pts */
  461.             if (h1==area[j-2] && h2==area[j-1])
  462.                continue;
  463.          }
  464.          xs = CENTERX + h1;
  465.          ys = CENTERY + h2;
  466.          plotpoint (xs, ys);
  467.          area[j] = h1;
  468.          area[j+1] = h2;
  469.          j += 2;
  470.          ++np;
  471.       }
  472.    }
  473.    if (np>0)  {                        /* check for last area */
  474.       pen = ORANGE;
  475.       SetAPen (rp, pen);
  476.       if (np==1)  {
  477.          xs = CENTERX + area[j-2];
  478.          ys = CENTERY + area[j-1];
  479.          plotpoint (xs, ys);
  480.       }
  481.       else  {
  482.          if (fill==FILL)
  483.             afill (np, pen);
  484.          else
  485.             adraw (np);
  486.       }
  487.    }
  488.    SetAPen (rp, ORANGE);
  489. }
  490.  
  491. /* ============================================================= */
  492.  
  493. void grid (val, lat0, lam0)            /* controls drawing of grids */
  494.  
  495. long val;
  496. double lat0, lam0;
  497.  
  498. {
  499.  
  500.    extern void globe_grid();
  501.    long oldpen, h1;
  502.    double lam, lat, temp;
  503.    static double lat_interval = 20.;
  504.    static double lam_interval = 30.;
  505.  
  506.    oldpen = rp->FgPen;                 /* save original pen color */
  507.    if (val!=FLAT && val!=MERCATOR)     /* draw grid for globe     */
  508.       globe_grid (val, lat0, lam0);
  509.    else  {                             /* otherwise grid for flat */
  510.                                        /*   or Mercator map       */
  511.       SetAPen (rp, BLACK);             /* set grid color to black */
  512.       for (lam=-180.; lam<=+180.; lam += lam_interval)  {
  513.          h1 = lam*HFACTOR + CENTERX;
  514.          Move (rp, h1, 0);
  515.          Draw (rp, h1, WHEIGHT-1);
  516.       }
  517.       for (lat=80.; lat>=-80.; lat -= lat_interval)  {
  518.          if (val==FLAT)
  519.             h1 = -lat*VFACTOR;
  520.          else
  521.             h1 = -log (tan((lat/2.+45.)*rad)) * M_VFACTOR;
  522.          h1 += CENTERY;
  523.          Move (rp, 0, h1);
  524.          Draw (rp, WWIDTH-1, h1);
  525.       }
  526.       SetAPen (rp, WHITE);             /* draw coordinate axes in white */
  527.       Move (rp, CENTERX, 0);
  528.       Draw (rp, CENTERX, WHEIGHT-1);
  529.       Move (rp, 0, CENTERY);
  530.       Draw (rp, WWIDTH-1, CENTERY);
  531.    }
  532.    SetAPen (rp, oldpen);               /* restore pen color */
  533. }
  534.  
  535. /* ============================================================= */
  536.  
  537. void globe_grid (val, lat0, lam0)      /* controls drawing globe grid */
  538.  
  539. double lat0, lam0;
  540. long val;
  541.  
  542. {
  543.  
  544.    extern void globe_grid_plot();
  545.    extern int limit_lam(), limit_lat();
  546.    double lat, lam, c0, s0, c1, s1, c2, s2, rlam, rlat;
  547.    double hp, scale, fac3, dlat, dlam, ddelt;
  548.    double lamp[2], latp[2];
  549.    int k;
  550.    long oldpen;
  551.    char first;
  552.    static double lat_interval = 20.;
  553.    static double lam_interval = 30.;
  554.    static double delta = 5.;
  555.    
  556.    lat0 *= rad;
  557.    c0 = cos (lat0);
  558.    s0 = sin (lat0);
  559.    dlat = lat_interval;
  560.    dlam = lam_interval;
  561.    ddelt = delta;
  562.    if (val==GLOBE)  {                  /* ordinary globe view */
  563.       hp = 0.;
  564.       scale = 1.;
  565.       fac3 = 1.;
  566.    }
  567.    else  {                             /* orbital view */
  568.       hp = etap;
  569.       scale = sqrt (1.-etap*etap);
  570.       fac3 = (facp/(facp-hp)) * scale;
  571.       if (view_height<=1200.)  {
  572.          dlat /= 4.;
  573.          dlam /= 4.;
  574.          ddelt /= 5.;
  575.       }
  576.    }
  577.    SetAPen (rp, BLACK);                /* grid lines in black */
  578.    for (lat=80.; lat>=-80.; lat-=dlat)  {  /* lines of equal latitude */
  579.       rlat = lat*rad;
  580.       if ((k=limit_lam (&lamp[0], rlat, lat0, hp))!=OK)
  581.          continue;                     /* skip if entirely out of view */
  582.       lamp[0] += lam0;
  583.       lamp[1] += lam0;
  584.       k = lamp[0]/ddelt - 1.;          /* express limits as multiple */
  585.       lamp[0] = k*ddelt;               /*   of ddelt                 */
  586.       k = lamp[1]/ddelt + 1.;
  587.       lamp[1] = k*ddelt;
  588.       c1 = cos (rlat);
  589.       s1 = sin (rlat);
  590.       first = TRUE;
  591.       if (lat==0.)
  592.          SetAPen (rp, WHITE);          /* draw equator in white */
  593.       for (lam=lamp[0]; lam<=lamp[1]; lam+=ddelt)  {
  594.          rlam = (lam-lam0)*rad;
  595.          if (rlam<-pi)
  596.             rlam += twopi;
  597.          if (rlam>+pi)
  598.             rlam -= twopi;
  599.          c2 = cos (rlam);
  600.          s2 = sin (rlam);
  601.          globe_grid_plot (rlat, rlam, c0, s0, c1, s1, c2, s2,
  602.                           val, hp, scale, fac3, &first);
  603.       }
  604.       if (lat==0.)
  605.          SetAPen (rp, BLACK);          /* reset pen to black */
  606.    }
  607.    for (lam=-180.; lam<+180.; lam+=dlam)  {  /* meridian circles */
  608.       rlam = (lam-lam0)*rad;
  609.       if (rlam<-pi)
  610.          rlam += twopi;
  611.       if (rlam>+pi)
  612.          rlam -= twopi;
  613.       if ((k=limit_lat (&latp[0], lat0, rlam, hp))!=OK)
  614.          continue;                     /* skip if entirely out of view */
  615.       k = latp[0]/ddelt + 1.;          /* express limits as multiple */
  616.       latp[0] = k*ddelt;               /*   of ddelt                 */
  617.       k = latp[1]/ddelt - 1;
  618.       latp[1] = k*ddelt;
  619.       if (latp[0]>=90.)                /* exclude North polar point */
  620.          latp[0] = 90. - ddelt;
  621.       if (latp[1]<=-90.)               /* exclude South polar point */
  622.          latp[1] = -90. + ddelt;
  623.       c2 = cos (rlam);
  624.       s2 = sin (rlam);
  625.       first = TRUE;
  626.       if (lam==0. || lam==-180.)       /* prime meridian in white */
  627.          SetAPen (rp, WHITE);
  628.       for (lat=latp[0]; lat>=latp[1]; lat-=ddelt)  {
  629.          rlat = lat*rad;
  630.          c1 = cos (rlat);
  631.          s1 = sin (rlat);
  632.          globe_grid_plot (rlat, rlam, c0, s0, c1, s1, c2, s2,
  633.                           val, hp, scale, fac3, &first);
  634.       }
  635.       if (lam==0. || lam==-180.)
  636.          SetAPen (rp, BLACK);          /* reset pen to black */
  637.    }
  638. }
  639.  
  640. /* ============================================================= */
  641.  
  642. void globe_grid_plot (lat, lam, c0, s0, c1, s1, c2, s2,
  643.                       val, hp, scale, fac3, first)
  644.  
  645. double lat, lam, c0, s0, c1, s1, c2, s2;   /* draws globe grids */
  646. long val;
  647. double hp, scale, fac3;
  648. char *first;
  649.  
  650. {
  651.  
  652.    extern void plotpoint(), globe_rim_fix();
  653.    char in_view;
  654.    short h1, h1c;                      /* x-dist. (pix) from center */
  655.    short h2, h2c;                      /* y-dist. (pix) from center */
  656.    short xs, ys;
  657.    long x, y;
  658.    double lamc, dlam;                  /* longitude */
  659.    double latc, dlat;                  /* latitude  */
  660.    double zp, facz, h1d, h2d, fac2;
  661.    static char prev_in_view;
  662.    static short h1prev, h2prev;
  663.    static double latprev, lamprev, zpprev;
  664.    
  665.    in_view = FALSE;                    /* get status of current point */
  666.    zp = s1*s0 + c1*c0*c2;
  667.    if (zp>=hp)  {                      /* zp > hp => in view */
  668.       in_view = TRUE;
  669.       h1d = HRADIUS * c1 * s2;
  670.       h2d = -VRADIUS * (s1*c0 - c1*s0*c2);
  671.       if (val!=GLOBE)  {
  672.          fac2 = (facp/(facp-zp))*scale;
  673.          h1d *= fac2;
  674.          h2d *= fac2;
  675.       }
  676.       h1 = h1d;
  677.       h2 = h2d;
  678.    }
  679.    if (*first==TRUE)  {                /* get status of first point */
  680.       *first = FALSE;
  681.       if (in_view==TRUE)  {            /* if first point is in view, */
  682.          x = h1 + CENTERX;             /*   move pen to first point  */
  683.          y = h2 + CENTERY;             /*   and plot it              */
  684.          Move (rp, x, y);
  685.          xs = x;
  686.          ys = y;
  687.          plotpoint (xs, ys);
  688.          h1prev = h1;
  689.          h2prev = h2;
  690.       }
  691.       prev_in_view = in_view;          /* save status of first point */
  692.       latprev = lat;
  693.       lamprev = lam;
  694.       zpprev = zp;
  695.       return;
  696.    }
  697.    if (in_view==TRUE)  {               /* if current point is in view, */
  698.       if (prev_in_view==FALSE)  {      /*   but previous point was not */
  699.          facz = zp / (zpprev-zp);      /*   in view, get rim point by  */
  700.          latc = lat - (latprev-lat)*facz; /* linear interpolation      */
  701.          lamc = lam - (lamprev-lam)*facz;
  702.          dlat = fabs (lat-latprev);
  703.          dlam = fabs (lam-lamprev);
  704.          if ( fabs (latc-latprev)> dlat || /* if rim point not between */
  705.               fabs (latc-lat)    > dlat)   /*   current and previous   */
  706.             latc = (lat+latprev)/2.;       /*   points, use midpoint   */
  707.          if ( fabs (lamc-lamprev)> dlam ||
  708.               fabs (lamc-lam)    > dlam )
  709.             lamc = (lam+lamprev)/2.;
  710.          c1 = cos (latc);
  711.          h1d = HRADIUS * c1 * sin (lamc);
  712.          h2d = -VRADIUS * (sin (latc)*c0 - c1*s0*cos (lamc));
  713.          if (val!=GLOBE)  {
  714.             h1d *= fac3;
  715.             h2d *= fac3;
  716.          }
  717.          h1c = h1d;
  718.          h2c = h2d;
  719.          globe_rim_fix (&h1c, &h2c);
  720.          x = h1c + CENTERX;            /* move to rim point & plot it */
  721.          y = h2c + CENTERY;
  722.          Move (rp, x, y);
  723.          xs = x;
  724.          ys = y;
  725.          plotpoint (xs, ys);
  726.          h1prev = h1c;
  727.          h2prev = h2c;
  728.       }
  729.       if (h1!=h1prev || h2!=h2prev)  {
  730.          x = h1 + CENTERX;             /* draw to current point */
  731.          y = h2 + CENTERY;
  732.          Draw (rp, x, y);
  733.          Move (rp, x, y);
  734.       }
  735.       h1prev = h1;
  736.       h2prev = h2;
  737.    }
  738.    else  {                             /* current point is out of view   */
  739.       if (prev_in_view==TRUE)  {       /* if previous point was in view, */
  740.          facz = zp / (zpprev-zp);      /*   get rim point by linear      */
  741.          latc = lat - (latprev-lat)*facz; /* interpolation               */
  742.          lamc = lam - (lamprev-lam)*facz;
  743.          dlat = fabs (lat-latprev);
  744.          dlam = fabs (lam-lamprev);
  745.          if ( fabs (latc-latprev)> dlat || /* if rim point not between */
  746.               fabs (latc-lat)    > dlat)   /*   current and previous   */
  747.             latc = (lat+latprev)/2.;       /*   points, use midpoint   */
  748.          if ( fabs (lamc-lamprev)> dlam ||
  749.               fabs (lamc-lam)    > dlam )
  750.             lamc = (lam+lamprev)/2.;
  751.          c1 = cos (latc);         
  752.          h1d = HRADIUS * c1 * sin (lamc);
  753.          h2d = -VRADIUS * (c0*sin (latc) - c1*s0 * cos (lamc));
  754.          if (val!=GLOBE)  {
  755.             h1d *= fac3;
  756.             h2d *= fac3;
  757.          }
  758.          h1c = h1d;
  759.          h2c = h2d;
  760.          globe_rim_fix (&h1c, &h2c);
  761.          x = h1c + CENTERX;            /* draw to rim point */
  762.          y = h2c + CENTERY;
  763.          Draw (rp, x, y);
  764.          Move (rp, x, y);
  765.       }
  766.    }
  767.    prev_in_view = in_view;             /* save status of current point */
  768.    latprev = lat;
  769.    lamprev = lam;
  770.    zpprev = zp;
  771. }
  772.  
  773. /* ============================================================= */
  774.  
  775. int limit_lam (lam, lat, lat0, etap)   /* computes limits on longitude */
  776.                                        /*   for constant latitude      */
  777. double *lam, lat, lat0, etap;
  778.  
  779. {
  780.    double alpha;
  781.    
  782.    alpha = (etap-sin(lat)*sin(lat0)) / (cos(lat)*cos(lat0));
  783.    if (alpha<=-1.)  {                  /* negative => lamda covers */
  784.       lam[0] = -180.;                  /*   entire hemisphere      */
  785.       lam[1] = +180.;
  786.       return (OK);
  787.    }
  788.    else if (alpha>=+1.)                /* positive => nothing in view */
  789.       return (NOT_OK);
  790.    else  {                             /* otherwise, compute limits */
  791.       lam[0] = -acos (alpha)/rad;
  792.       lam[1] = -lam[0];
  793.       return (OK);
  794.    }
  795. }
  796.  
  797. /* ============================================================= */
  798.  
  799. int limit_lat (lat, lat0, lam, etap)   /* computes limits on latitude */
  800.                                        /*   for constant longitude    */
  801. double *lat, lat0, lam, etap;
  802.  
  803. {
  804.    double radical, a, b, sum, fac1;
  805.    
  806.    a = sin (lat0);
  807.    b = cos (lat0) * cos (lam);
  808.    sum = a*a + b*b;
  809.    radical = sum - etap*etap;
  810.    if (radical<=0.)                    /* no real solutions */
  811.       return (NOT_OK);
  812.    else  {                             /* two real solutions      */
  813.       radical = sqrt (radical);        /* solve quadrtic equation */
  814.       fac1 = (a*etap + b*radical)/sum;
  815.       lat[0] = asin (fac1)/rad;
  816.       fac1 = (a*etap - b*radical)/sum;
  817.       lat[1] = asin (fac1)/rad;
  818.       if (lat[0]<lat[1])  {            /* put in correct order */
  819.          b = lat[0];
  820.          lat[0] = lat[1];
  821.          lat[1] = b;
  822.       }
  823.       if (a>etap)                      /* check North pole */
  824.          lat[0] = 90.;
  825.       if (a<-etap)                     /* check South pole */
  826.          lat[1] = -90.;
  827.       return (OK);
  828.    }
  829. }
  830.  
  831. /* ============================================================= */
  832.  
  833. void getcoord (x, y, val, lat, lam)    /* converts screen coordinates   */
  834.                                        /*   into latitude and longitude */
  835. int x, y;
  836. long val;
  837. double *lat, *lam;
  838.  
  839. {
  840.  
  841.    (*lam) = (x - CENTERX) / HFACTOR;
  842.    if (val==FLAT)
  843.       (*lat) = (CENTERY - y) / VFACTOR;
  844.    else
  845.       (*lat) = -90. + 2.*atan(exp((CENTERY-y)/M_VFACTOR))/rad;
  846. }
  847.  
  848. /* ============================================================= */
  849.  
  850. void globe (ws, ws_trig, lat0, lam0, val, fill) /* draws globe projections */
  851.  
  852. short *ws, *ws_trig;
  853. double lat0, lam0;
  854. long val;
  855. int fill;
  856.  
  857. {
  858.  
  859.    extern void plotpoint(), globe_fill(), globe_rim_fix();
  860.    char first, prev_in_view, in_view, latzero;
  861.    short h1, h1c, h1prev;              /* x-dist. (pix) from center */
  862.    short h2, h2c, h2prev;              /* y-dist. (pix) from center */
  863.    short np, narea, nrim_in, nrim_out;
  864.    short xs, ys;
  865.    int i, j;
  866.    long x, y, pen, first_color;
  867.    double lam, lamc, lamprev;          /* longitude */
  868.    double lat, latc, latprev;          /* latitude  */
  869.    double c0, s0, c1, s1, c2, zp, zpprev;
  870.    double hp, fac2, fac3, facz, scale;
  871.    double h1d, h2d, dlat, dlam, lam0p;
  872.    static double fac = 1./32767.0;
  873.  
  874.    latzero = FALSE;
  875.    if (val==GLOBE)  {                  /* ordinary globe view */
  876.       hp = 0.;
  877.       scale = 1.;
  878.       fac3 = 1.;
  879.    }
  880.    else  {                             /* orbital globe view */
  881.       hp = etap;
  882.       scale = sqrt (1.-etap*etap);
  883.       fac3 = (facp/(facp-hp)) * scale;
  884.    }
  885.    if (lat0==0.)  {
  886.       c0 = 1.;
  887.       s0 = 0.;
  888.       if (val==GLOBE)
  889.          latzero = TRUE;               /* equatorial, ordinary globe view */
  890.    }
  891.    else  {
  892.       lat0 *= rad;
  893.       c0 = cos (lat0);
  894.       s0 = sin (lat0);
  895.    }
  896.    first = TRUE;
  897.    pen = ORANGE;
  898.    SetAPen (rp, pen);
  899.    first_color = BLACK;
  900.    j = 0;
  901.    np = 0;
  902.    narea = 0;
  903.    nrim_in = 0;
  904.    nrim_out = 0;
  905.    lam0p = lam0 * rad;                 /* convert lam0 to radians */
  906.    for (i=0; i<MAXVAL; i+=2)  {
  907.       if (ws[i]==0 && ws[i+1]==0)  {   /* if end of area, then     */
  908.          ++narea;                      /*   color-fill and skip to */
  909.          first = TRUE;                 /*    next area             */
  910.          if (np>0 && fill==FILL)
  911.             globe_fill (np, narea, nrim_in, nrim_out, lat0, lam0p,
  912.                         first_color);
  913.          j = 0;
  914.          np = 0;
  915.          nrim_in = 0;
  916.          nrim_out = 0;
  917.          first_color = BLACK;
  918.          continue;
  919.       }
  920.       lat = ws[i];                     /* latitude */
  921.       lat = (lat/100.) * rad;
  922.       lam = ws[i+1];                   /* longitude */
  923.       lam = (lam/100. - lam0) * rad;
  924.       if (lam<-pi)
  925.          lam += twopi;
  926.       if (lam>pi)
  927.          lam -= twopi;
  928.       c1 = ws_trig[i];                 /* cosine of latitude */
  929.       c1 *= fac;
  930.       s1 = ws_trig[i+1];               /* sine of latitude */
  931.       s1 *= fac;
  932.       in_view = FALSE;                 /* get status of current point */
  933.       if (latzero==TRUE)  {            /* equatorial globe view */
  934.          zp = c1*cos (lam);
  935.          if (lam>=-pi2 && lam<=+pi2)  {
  936.             in_view = TRUE;
  937.             h1 = HRADIUS * c1 * sin (lam);
  938.             h2 = -VRADIUS * s1;
  939.          }
  940.       }
  941.       else  {                          /* oblique earth view */
  942.          c2 = cos (lam);
  943.          zp = s1*s0 + c1*c0*c2;
  944.          if (zp>=hp)  {                /* zp > hp => in view */
  945.             in_view = TRUE;
  946.             h1d = HRADIUS * c1 * sin (lam);
  947.             h2d = -VRADIUS * (s1*c0 -c1*s0*c2);
  948.             if (val!=GLOBE)  {
  949.                fac2 = (facp/(facp-zp)) * scale;
  950.                h1d *= fac2;
  951.                h2d *= fac2;
  952.             }
  953.             h1 = h1d;
  954.             h2 = h2d;
  955.          }
  956.       }
  957.       if (first==TRUE)  {              /* get status of first point */
  958.          first = FALSE;
  959.          if (in_view==TRUE)  {         /* if first point is in view, */
  960.             x = h1 + CENTERX;          /*   move pen to first point  */
  961.             y = h2 + CENTERY;          /*   and plot it              */
  962.             Move (rp, x, y);
  963.             first_color = ReadPixel (rp, x, y);
  964.             xs = x;
  965.             ys = y;
  966.             plotpoint (xs, ys);
  967.             h1prev = h1;
  968.             h2prev = h2;
  969.             area[0] = h1;              /* save first point for later use */
  970.             area[1] = h2;
  971.             j = 2;
  972.             np = 1;
  973.          }
  974.          prev_in_view = in_view;       /* save status of first point */
  975.          latprev = lat;
  976.          lamprev = lam;
  977.          zpprev = zp;
  978.          continue;
  979.       }
  980.       if (in_view==TRUE)  {            /* if current point is in view, */
  981.          if (prev_in_view==FALSE)  {   /*   but previous point was not */
  982.             facz = zp / (zpprev-zp);   /*   in view, get rim point by  */
  983.             latc = lat - (latprev-lat)*facz; /* linear interpolation   */
  984.             lamc = lam - (lamprev-lam)*facz;
  985.             dlat = fabs (lat-latprev);
  986.             dlam = fabs (lam-lamprev);
  987.             if ( fabs (latc-latprev)> dlat || /* if rim point not between */
  988.                  fabs (latc-lat)    > dlat)   /*   current and previous   */
  989.                latc = (lat+latprev)/2.;       /*   point, use midpoint    */
  990.             if ( fabs (lamc-lamprev)> dlam ||
  991.                  fabs (lamc-lam)    > dlam )
  992.                lamc = (lam+lamprev)/2.;
  993.             if (latzero==TRUE)  {
  994.                h1c = HRADIUS * cos (latc) * sin (lamc);
  995.                h2c = -VRADIUS * sin (latc);
  996.             }
  997.             else  {
  998.                c1 = cos (latc);
  999.                h1d = HRADIUS * c1 * sin (lamc);
  1000.                h2d = -VRADIUS * (sin (latc)*c0 - c1*s0*cos (lamc));
  1001.                if (val!=GLOBE)  {
  1002.                   h1d *= fac3;
  1003.                   h2d *= fac3;
  1004.                }
  1005.                h1c = h1d;
  1006.                h2c = h2d;
  1007.             }
  1008.             globe_rim_fix (&h1c, &h2c); /* correct the computed rim point */
  1009.             x = h1c + CENTERX;         /* move to rim point & plot it */
  1010.             y = h2c + CENTERY;
  1011.             Move (rp, x, y);
  1012.             xs = x;
  1013.             ys = y;
  1014.             plotpoint (xs, ys);
  1015.             h1prev = h1c;
  1016.             h2prev = h2c;
  1017.             area[j] = h1;              /* save coords of rim point */
  1018.             area[j+1] = h2;
  1019.             j += 2;
  1020.             if (nrim_in<MAXRIM)  {
  1021.                rim_in[nrim_in] = np;
  1022.                ++nrim_in;
  1023.             }
  1024.             ++np;
  1025.          }
  1026.          if (h1!=h1prev || h2!=h2prev)  {
  1027.             x = h1 + CENTERX;          /* draw to current point */
  1028.             y = h2 + CENTERY;
  1029.             Draw (rp, x, y);
  1030.             Move (rp, x, y);
  1031.             area[j] = h1;              /* save coords of current point */
  1032.             area[j+1] = h2;
  1033.             j += 2;
  1034.             ++np;
  1035.          }
  1036.          h1prev = h1;
  1037.          h2prev = h2;
  1038.       }
  1039.       else  {                          /* current point is out of view   */
  1040.          if (prev_in_view==TRUE)  {    /* if previous point was in view, */
  1041.             facz = zp / (zpprev-zp);   /*   get rim point by linear      */
  1042.             latc = lat - (latprev-lat)*facz; /* interpolation            */
  1043.             lamc = lam - (lamprev-lam)*facz;
  1044.             dlat = fabs (lat-latprev);
  1045.             dlam = fabs (lam-lamprev);
  1046.             if ( fabs (latc-latprev)> dlat || /* if rim point not between */
  1047.                  fabs (latc-lat)    > dlat)   /*   current and previous   */
  1048.                latc = (lat+latprev)/2.;       /*   point, use midpoint    */
  1049.             if ( fabs (lamc-lamprev)> dlam ||
  1050.                  fabs (lamc-lam)    > dlam )
  1051.                lamc = (lam+lamprev)/2.;
  1052.             if (latzero==TRUE)  {
  1053.                h1c = HRADIUS * cos (latc) * sin (lamc);
  1054.                h2c = -VRADIUS * sin (latc);
  1055.             }
  1056.             else  {
  1057.                c1 = cos (latc);
  1058.                h1d = HRADIUS * c1 * sin (lamc);
  1059.                h2d = -VRADIUS * (c0*sin(latc) - c1*s0 * cos(lamc));
  1060.                if (val!=GLOBE)  {
  1061.                   h1d *= fac3;
  1062.                   h2d *= fac3;
  1063.                }
  1064.                h1c = h1d;
  1065.                h2c = h2d;
  1066.             }
  1067.             globe_rim_fix (&h1c, &h2c); /* correct the computed rim point */ 
  1068.             x = h1c + CENTERX;         /* draw to rim point */
  1069.             y = h2c + CENTERY;
  1070.             Draw (rp, x, y);
  1071.             Move (rp, x, y);
  1072.             area[j] = h1;              /* save coords of rim point */
  1073.             area[j+1] = h2;
  1074.             j += 2;
  1075.             if (nrim_out<MAXRIM)  {
  1076.                rim_out[nrim_out] = np;
  1077.                ++nrim_out;
  1078.             }
  1079.             ++np;
  1080.          }
  1081.       }
  1082.       prev_in_view = in_view;          /* save status of current point */
  1083.       latprev = lat;
  1084.       lamprev = lam;
  1085.       zpprev = zp;
  1086.    }
  1087. }
  1088.  
  1089. /* ============================================================= */
  1090.  
  1091. void globe_fill (np, narea, nrim_in, nrim_out, lat0, lam0, first_color)
  1092.                                        /* fills area on globe */
  1093.  
  1094. short np, narea, nrim_in, nrim_out;
  1095. double lat0, lam0;
  1096. long first_color;
  1097.  
  1098. {
  1099.  
  1100.    extern void plotpoint ();
  1101.    int k;
  1102.    long x, y;
  1103.    short xs, ys;
  1104.    char is_lake;
  1105.    static long pen = ORANGE;
  1106.    static char area1 = FALSE;
  1107.  
  1108.    if (nrim_in!=nrim_out)
  1109.       return;
  1110.    if (nrim_in==0)  {                  /* fill areas with no rim points */
  1111.       is_lake = FALSE;
  1112.       for (k=0; k<num_lakes; ++k)      /* check for lake */
  1113.          if (narea==lakes[k])  {
  1114.             if (first_color==ORANGE)  {
  1115.                pen = BLUE;             /* make lake blue if drawn */
  1116.                is_lake = TRUE;         /*   into orange           */
  1117.                break;
  1118.             }
  1119.             else                       /* don't color-fill if drawn */
  1120.                return;                 /*   into blue               */
  1121.          }
  1122.       xs = CENTERX + area[0];
  1123.       ys = CENTERY + area[1];
  1124.       if (np==1)
  1125.          plotpoint (xs, ys);
  1126.       else if (np>=GLOBE_FILL_MIN)
  1127.          afill (np, pen);
  1128.       if (narea==1)                    /* check south polar region */
  1129.          area1 = TRUE;
  1130.       else if (narea==2)  {
  1131.          if (area1==TRUE)  {
  1132.             area1 = FALSE;
  1133.             if (lat0<0.)  {
  1134.                x = CENTERX;
  1135.                y = CENTERY + VRADIUS * cos (lat0);
  1136.                Flood (rp, 1, x, y);
  1137.             }
  1138.          }
  1139.       }
  1140.       if (is_lake==TRUE)
  1141.          pen = ORANGE;
  1142.    }
  1143. }
  1144.  
  1145. /* ============================================================= */
  1146.  
  1147. void globe_rim_fix (h1, h2)            /* find nearest actual rim point */
  1148.  
  1149. short *h1, *h2;
  1150.  
  1151. {
  1152.  
  1153.    short inc1, inc2;
  1154.    int i;
  1155.    long x, y, color;
  1156.    static int itmax = 20;
  1157.    
  1158.    if ((*h1)>=0)                       /* right half */
  1159.       inc1 = +1;
  1160.    else                                /* left half */
  1161.       inc1 = -1;
  1162.    if ((*h2)>=0)                       /* bottom half */
  1163.       inc2 = +1;
  1164.    else                                /* top half */
  1165.       inc2 = -1;
  1166.    x = (*h1) + CENTERX + 5*inc1;       /* coords of test pixel */
  1167.    y = (*h2) + CENTERY + 5*inc2;
  1168.    for (i=0; i<itmax; ++i)  {          /* look for nearest black pixel */
  1169.       if ((color = ReadPixel (rp, x, y))==BLACK)
  1170.          break;
  1171.       x += inc1;
  1172.       y += inc2;
  1173.    }
  1174.    x -= inc1;
  1175.    y -= inc2;
  1176.    for (i=0; i<itmax; ++i)  {          /* look back for previous blue one */
  1177.       if ((color = ReadPixel (rp, x, y))==BLUE || color==ORANGE)
  1178.          break;
  1179.       x -= inc1;
  1180.       y -= inc2;
  1181.    }
  1182.    (*h1) = x - CENTERX;
  1183.    (*h2) = y - CENTERY;
  1184. }
  1185.  
  1186. /* ============================================================= */
  1187.  
  1188. void stars ()                          /* draws stars into black background */
  1189.  
  1190. {
  1191.  
  1192.    extern double ran();
  1193.    static int nmax = 150;
  1194.    double t;
  1195.    short x, y;
  1196.    int i, nstars;
  1197.    long oldpen, color;
  1198.   
  1199.    oldpen = rp->FgPen;
  1200.    SetAPen (rp, WHITE);
  1201.    nstars = 0;
  1202.    for (i=0; ; ++i)  {
  1203.       t = ran();
  1204.       x = t * WWIDTH + 0.5;
  1205.       t = ran();
  1206.       y = t * WHEIGHT + 0.5;
  1207.       if (x<0)
  1208.          x = 0;
  1209.       if (x>WWIDTH-1)
  1210.          x = WWIDTH-1;
  1211.       if (y<0)
  1212.          y = 0;
  1213.       if (y>WHEIGHT-1)
  1214.          y = WHEIGHT-1;
  1215.       if ((color = ReadPixel (rp, x, y)) == BLACK)  {
  1216.          WritePixel (rp, x, y);
  1217.          ++nstars;
  1218.          if (nstars>=nmax)
  1219.             break;
  1220.       }
  1221.    }
  1222.    SetAPen (rp, oldpen);
  1223. }
  1224.  
  1225. /* ============================================================= */
  1226.  
  1227. void floodfill ()                      /* flood fills an area based */
  1228.                                        /*   on mouse position       */
  1229. {
  1230.    struct IntuiMessage *msgf;
  1231.    short x, y;
  1232.    long color;
  1233.    
  1234.    WaitPort (w->UserPort);             /* wait for message */
  1235.    while (1)  {
  1236.       msgf = (struct IntuiMessage *) GetMsg (w->UserPort);
  1237.       if (msgf==NULL)
  1238.          WaitPort (w->UserPort);
  1239.       else
  1240.          if (msgf->Code==SELECTDOWN)
  1241.             break;
  1242.    }
  1243.    x = msgf->MouseX;                   /*  get mouse coordinates  */
  1244.    y = msgf->MouseY;                   /*    and flood fill       */
  1245.    if ((color = ReadPixel (rp, x, y))==BLUE)
  1246.       Flood (rp, 1, x, y);
  1247. }
  1248.  
  1249. /* ============================================================= */
  1250.  
  1251. int getbox (x0, y0, x, y)              /* selects a region to draw to */
  1252.                                        /*   larger scale              */
  1253. int *x0, *y0, *x, *y;
  1254.  
  1255. {
  1256.  
  1257.    extern void drawbox();
  1258.    struct IntuiMessage *msgf;
  1259.    int xd, yd, x1, x2, y1, y2;
  1260.    char selectbutton;
  1261.    
  1262.    selectbutton = FALSE;
  1263.    SetDrMd (rp, JAM2 | COMPLEMENT);    /* turn on complement mode and */
  1264.    ModifyIDCMP (w, IDCMPFLAGS | MOUSEMOVE); /* enable mouse events    */
  1265.    WaitPort (w->UserPort);
  1266.    while (1)  {
  1267.       msgf = (struct IntuiMessage *) GetMsg (w->UserPort);
  1268.       if (msgf==NULL)
  1269.          WaitPort (w->UserPort);
  1270.       else if (msgf->Code==SELECTDOWN)  { /* if user pressed left button, */
  1271.          x1 = msgf->MouseX;               /*   get initial mouse position */
  1272.          y1 = msgf->MouseY;
  1273.          xd = x1;
  1274.          yd = y1;
  1275.          selectbutton = TRUE;
  1276.       }
  1277.       else if (selectbutton==TRUE && msgf->Class==MOUSEMOVE)  {
  1278.          x2 = msgf->MouseX;
  1279.          y2 = msgf->MouseY;
  1280.          drawbox (x1, y1, xd, yd);     /* erase old box */
  1281.          xd = x2;
  1282.          yd = y2;
  1283.          drawbox (x1, y1, xd, yd);     /* draw new box */
  1284.       }
  1285.       else if (selectbutton==TRUE && msgf->Code==SELECTUP)
  1286.          break;
  1287.    }
  1288.    x2 = msgf->MouseX;                  /* erase current box */
  1289.    y2 = msgf->MouseY;
  1290.    drawbox (x1, y1, xd, yd);
  1291.    SetDrMd (rp, JAM2);                 /* restore original drawing mode */
  1292.    ModifyIDCMP (w, IDCMPFLAGS);        /*   and disable mouse events    */
  1293.    *x0 = x1;
  1294.    *y0 = y1;
  1295.    *x = x2;
  1296.    *y = y2;
  1297.    if (x1==x2 && y1==y2)               /* error if box is of zero size */
  1298.       return (NOT_OK);
  1299.    if (x1>x2 && y1>y2)  {              /* ensure that the vertices of */
  1300.       *x0 = x2;                        /*   the box are in the proper */
  1301.       *y0 = y2;                        /*   order                     */
  1302.       *x = x1;
  1303.       *y = y1;
  1304.    }
  1305.    else if (x1<x2 && y1>y2)  {
  1306.       *x0 = x1;
  1307.       *y0 = y2;
  1308.       *x = x2;
  1309.       *y = y1;
  1310.    }
  1311.    else if (x1>x2 && y1<y2)  {
  1312.       *x0 = x2;
  1313.       *y0 = y1;
  1314.       *x = x1;
  1315.       *y = y2;
  1316.    }
  1317.    return (OK);
  1318. }
  1319.  
  1320. /* ============================================================= */
  1321.  
  1322. void drawbox (x1, y1, x2, y2)          /* draws a box */
  1323.  
  1324. int x1, y1, x2, y2;
  1325.  
  1326. {
  1327.  
  1328.    Move (rp, x1, y1);
  1329.    Draw (rp, x2, y1);
  1330.    Move (rp, x2, y1);
  1331.    Draw (rp, x2, y2);
  1332.    Move (rp, x2, y2);
  1333.    Draw (rp, x1, y2);
  1334.    Move (rp, x1, y2);
  1335.    Draw (rp, x1, y1);
  1336.    Move (rp, x1, y1);
  1337. }
  1338.  
  1339. /* ============================================================= */
  1340.  
  1341. void box (ws, latp, lamp, fill)        /* draws areas contained within */
  1342.                                        /*   a rectangular region       */
  1343. short *ws;
  1344. double *latp, *lamp;
  1345. int fill;
  1346.  
  1347. {
  1348.  
  1349.    extern void plotpoint(), drawbox();
  1350.    char first, prev_in_view, in_view, is_lake;
  1351.    short h1, h1c, h1prev;              /* x-dist. (pix) from center */
  1352.    short h2, h2c, h2prev;              /* y-dist. (pix) from center */
  1353.    short np, xs, ys, nrim_in, nrim_out;
  1354.    int i, j, k, narea;
  1355.    long x, y, pen, first_color;
  1356.    double lam, lamc, lamprev;          /* longitude */
  1357.    double lat, latc, latprev;          /* latitude  */
  1358.    double bwidth, bheight, bcx, bcy, xscale, yscale;
  1359.    double lat1, lat2, lam1, lam2;
  1360.  
  1361.    lat1 = latp[0];                     /* store values for box corners */
  1362.    lat2 = latp[1];                     /*   locally                    */
  1363.    lam1 = lamp[0];
  1364.    lam2 = lamp[1];
  1365.    bwidth = lam2 - lam1;               /* box width (degrees)         */
  1366.    bheight = lat1 - lat2;              /* box height (degrees)        */
  1367.    bcx = (lam1 + lam2)/2.;             /* x-coord of box center (deg) */
  1368.    bcy = (lat1 + lat2)/2.;             /* y-coord of box center (deg) */
  1369.    xscale = WWIDTH / bwidth;           /* horizontal scale (pix/deg)  */
  1370.    yscale = (WHEIGHT-10) / bheight;    /* vertical scale (pix/deg)    */
  1371.    j = 0;
  1372.    np= 0;
  1373.    narea = 0;
  1374.    first = TRUE;
  1375.    pen = ORANGE;
  1376.    SetAPen (rp, pen);
  1377.    first_color = BLACK;
  1378.  
  1379.    drawbox (0, 0, WWIDTH-1, WHEIGHT-10); /* outline the box */
  1380.  
  1381.    for (i=0; i<MAXVAL; i+=2)  {
  1382.       if (ws[i]==0 && ws[i+1]==0)  {   /* if end of area, then */
  1383.          ++narea;
  1384.          if (np>0 && fill==FILL)  {    /* color-fill if requested */
  1385.             if (nrim_in!=nrim_out)
  1386.                ;
  1387.             else if (nrim_in==0)  {
  1388.                is_lake = FALSE;
  1389.                for (k=0; k<num_lakes; ++k)
  1390.                   if (narea==lakes[k])  {
  1391.                      is_lake = TRUE;
  1392.                      if (first_color==ORANGE)
  1393.                         pen = BLUE;
  1394.                      break;
  1395.                   }
  1396.                if (is_lake==FALSE || (is_lake==TRUE && pen==BLUE))  {
  1397.                   xs = CENTERX + area[0];
  1398.                   ys = CENTERY + area[1];
  1399.                   if (np==1)
  1400.                      plotpoint (xs, ys);
  1401.                   else
  1402.                      afill (np, pen);
  1403.                }
  1404.             }
  1405.          }
  1406.          j = 0;
  1407.          np = 0;
  1408.          nrim_in = 0;
  1409.          nrim_out = 0;
  1410.          pen = ORANGE;
  1411.          first_color = BLACK;
  1412.          first = TRUE;
  1413.          continue;                     /* skip to next point */
  1414.       }
  1415.       lat = ws[i];                     /* latitude */
  1416.       lat = lat/100.;
  1417.       lam = ws[i+1];                   /* longitude */
  1418.       lam = lam/100.;
  1419.       in_view = FALSE;                 /* get status of current point */
  1420.       if ( (lat<=lat1 && lat>=lat2)
  1421.         && (lam>=lam1 && lam<=lam2))  {
  1422.          in_view = TRUE;
  1423.          h1 = (lam-bcx) * xscale;
  1424.          h2 = - (lat-bcy) * yscale;
  1425.       }
  1426.       if (first==TRUE)  {              /* check first point */
  1427.          first = FALSE;
  1428.          if (in_view==TRUE)  {         /* if first point is in view, */
  1429.             x = h1 + CENTERX;          /*   move pen to first point  */
  1430.             y = h2 + CENTERY;          /*   and plot it              */
  1431.             Move (rp, x, y);
  1432.             first_color = ReadPixel (rp, x, y);
  1433.             xs = x;
  1434.             ys = y;
  1435.             plotpoint (xs, ys);
  1436.             h1prev = h1;
  1437.             h2prev = h2;
  1438.             area[0] = h1;
  1439.             area[1] = h2;
  1440.             j = 2;
  1441.             np = 1;
  1442.          }
  1443.          prev_in_view = in_view;       /* save status of first point */
  1444.          latprev = lat;
  1445.          lamprev = lam;
  1446.          continue;
  1447.       }
  1448.       if (in_view==TRUE)  {
  1449.          if (prev_in_view==FALSE)  {   /* if prev. point was not in view, */
  1450.             if (lamprev<=lam1)  {      /*   find rim point                */
  1451.                lamc = lam1;
  1452.                if (latprev<=lat2)      /* lower left */
  1453.                   latc = lat2;
  1454.                else if (latprev>=lat1) /* upper left */
  1455.                   latc = lat1;
  1456.                else                    /* left center */
  1457.                   latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
  1458.             }
  1459.             else if (lamprev>=lam2)  {
  1460.                lamc = lam2;
  1461.                if (latprev>=lat1)      /* upper right */
  1462.                   latc = lat1;
  1463.                else if (latprev<=lat2) /* lower right */
  1464.                   latc = lat2;
  1465.                else                    /* right center */
  1466.                   latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
  1467.             }
  1468.             else  {
  1469.                if (latprev>=lat1)      /* top center */
  1470.                   latc = lat1;
  1471.                else                    /* bottom center */
  1472.                   latc = lat2;
  1473.                lamc = lam + (lamprev-lam)*(latc-lat)/(latprev-lat);
  1474.             }
  1475.             h1c = (lamc-bcx) * xscale;
  1476.             h2c = - (latc-bcy) * yscale;
  1477.             x = h1c + CENTERX;         /* move to rim point & plot it */
  1478.             y = h2c + CENTERY;
  1479.             Move (rp, x, y);
  1480.             xs = x;
  1481.             ys = y;
  1482.             plotpoint (xs, ys);
  1483.             h1prev = h1c;
  1484.             h2prev = h2c;
  1485.             area[j] = h1;
  1486.             area[j+1] = h2;
  1487.             j += 2;
  1488.             if (nrim_in<MAXRIM)  {
  1489.                rim_in[nrim_in] = np;
  1490.                ++nrim_in;
  1491.             }
  1492.             ++np;
  1493.          }
  1494.          if (h1!=h1prev || h2!=h2prev)  {
  1495.             x = h1 + CENTERX;          /* draw to current point */
  1496.             y = h2 + CENTERY;
  1497.             Draw (rp, x, y);
  1498.             Move (rp, x, y);
  1499.             area[j] = h1;
  1500.             area[j+1] = h2;
  1501.             j += 2;
  1502.             ++np;
  1503.          }
  1504.          h1prev = h1;
  1505.          h2prev = h2;
  1506.       }
  1507.       else  {                          /* else out of view */
  1508.          if (prev_in_view==TRUE)  {    /* if previous point was in view, */
  1509.             if (lam<=lam1)  {          /*   find rim point                */
  1510.                lamc = lam1;
  1511.                if (lat<=lat2)          /* lower left */
  1512.                   latc = lat2;
  1513.                else if (lat>=lat1)     /* upper left */
  1514.                   latc = lat1;
  1515.                else                    /* left center */
  1516.                   latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
  1517.             }
  1518.             else if (lam>=lam2)  {
  1519.                lamc = lam2;
  1520.                if (lat>=lat1)          /* upper right */
  1521.                   latc = lat1;
  1522.                else if (lat<=lat2)     /* lower right */
  1523.                   latc = lat2;
  1524.                else                    /* right center */
  1525.                   latc = lat + (latprev-lat)*(lamc-lam)/(lamprev-lam);
  1526.             }
  1527.             else  {
  1528.                if (lat>=lat1)          /* top center */
  1529.                   latc = lat1;
  1530.                else                    /* bottom center */
  1531.                   latc = lat2;
  1532.                lamc = lam + (lamprev-lam)*(latc-lat)/(latprev-lat);
  1533.             }
  1534.             h1c = (lamc-bcx) * xscale;
  1535.             h2c = - (latc-bcy) * yscale;
  1536.             x = h1c + CENTERX;         /* draw to rim point */
  1537.             y = h2c + CENTERY;
  1538.             Draw (rp, x, y);
  1539.             Move (rp, x, y);
  1540.             area[j] = h1;
  1541.             area[j+1] = h2;
  1542.             j += 2;
  1543.             if (nrim_out<MAXRIM)  {
  1544.                rim_out[nrim_out] = np;
  1545.                ++nrim_out;
  1546.             }
  1547.             ++np;
  1548.          }
  1549.       }
  1550.       prev_in_view = in_view;          /* save status of current point */
  1551.       latprev = lat;
  1552.       lamprev = lam;
  1553.    }
  1554. }
  1555.  
  1556. /* ============================================================= */
  1557.  
  1558. void afill (pairs, pen)                /* draws and fills a region  */
  1559.                                        /*   using AreaDraw function */
  1560. short pairs;                            /* how many pairs of words */
  1561. long pen;                              /* pen number to use       */
  1562.  
  1563. {
  1564.  
  1565.    short i, j;
  1566.    long x, y, oldpen;
  1567.  
  1568.    oldpen = rp->FgPen;
  1569.    SetAPen (rp, pen);
  1570.    x = area[0] + CENTERX;
  1571.    y = area[1] + CENTERY;
  1572.    if (x<0)
  1573.       x = 0;
  1574.    if (x>WWIDTH-1)
  1575.       x = WWIDTH-1;
  1576.    if (y<0)
  1577.       y = 0;
  1578.    if (y>WHEIGHT-1)
  1579.       y = WHEIGHT-1;
  1580.    AreaMove (rp, x, y);
  1581.    j = 2;
  1582.    for (i=1; i<pairs; i++)  {
  1583.       if (area[j]==0 && area[j+1]==0)
  1584.          break;
  1585.       x = area[j] + CENTERX;
  1586.       y = area[j+1] + CENTERY;
  1587.       if (x<0)
  1588.          x = 0;
  1589.       if (x>WWIDTH-1)
  1590.          x = WWIDTH-1;
  1591.       if (y<0)
  1592.          y = 0;
  1593.       if (y>WHEIGHT-1)
  1594.          y = WHEIGHT-1;
  1595.       AreaDraw (rp, x, y);
  1596.       j += 2;
  1597.    }   
  1598.    AreaEnd (rp);
  1599.    SetAPen (rp, oldpen);
  1600. }
  1601.  
  1602. /* ============================================================= */
  1603.  
  1604. void adraw (pairs)                     /* draws area outlines */
  1605.  
  1606. short pairs;
  1607.  
  1608. {
  1609.  
  1610.    short i, j;
  1611.    long x, y, oldpen;
  1612.    
  1613.    oldpen = rp->FgPen;
  1614.    SetAPen (rp, ORANGE);
  1615.    x = area[0] + CENTERX;
  1616.    y = area[1] + CENTERY;
  1617.    if (x<0)
  1618.       x = 0;
  1619.    if (x>WWIDTH-1)
  1620.       x = WWIDTH-1;
  1621.    if (y<0)
  1622.       y = 0;
  1623.    if (y>WHEIGHT-1)
  1624.       y = WHEIGHT-1;
  1625.    Move (rp, x, y);
  1626.    j = 2;
  1627.    for (i=1; i<pairs; ++i)  {
  1628.       if (area[j]==0 && area[j+1]==0)
  1629.          break;
  1630.       x = area[j] + CENTERX;
  1631.       y = area[j+1] + CENTERY;
  1632.       if (x<0)
  1633.          x = 0;
  1634.       if (x>WWIDTH-1)
  1635.          x = WWIDTH-1;
  1636.       if (y<0)
  1637.          y = 0;
  1638.       if (y>WHEIGHT-1)
  1639.          y = WHEIGHT-1;
  1640.       Draw (rp, x, y);
  1641.       Move (rp, x, y);
  1642.       j += 2;
  1643.    }
  1644.    SetAPen (rp, oldpen);
  1645. }
  1646.  
  1647. /* ============================================================= */
  1648.  
  1649. void plotpoint (x, y)                  /* plots a single point */
  1650.  
  1651. short x, y;
  1652.  
  1653. {
  1654.  
  1655.    if (x<0)
  1656.       x = 0;
  1657.    if (x>WWIDTH-1)
  1658.       x = WWIDTH-1;
  1659.    if (y<0)
  1660.       y = 0;
  1661.    if (y>WHEIGHT-1)
  1662.       y = WHEIGHT-1;
  1663.    WritePixel (rp, x, y);
  1664. }
  1665.  
  1666. /* ============================================================= */
  1667.  
  1668. short *svector (nl,nh)                 /* Allocates a short vector */
  1669.                                        /*   with range [nl..nh]    */
  1670. int nl, nh;
  1671.  
  1672. {
  1673.    extern char *calloc();
  1674.    extern void errors();
  1675.    short *v;
  1676.  
  1677.    v = (short*) calloc (1, (unsigned) (nh-nl+1)*sizeof(short));
  1678.    if (v==NULL)
  1679.       return (NULL);
  1680.    else
  1681.       return v-nl;
  1682. }
  1683.  
  1684. /* ============================================================= */
  1685.  
  1686. void free_svector (v,nl,nh)            /* Frees a short vector   */
  1687.                                        /*   allocated by svector */
  1688. short *v;
  1689. int nl, nh;
  1690.  
  1691. {
  1692.  
  1693.    free ((char*) (v+nl));
  1694. }
  1695.