Garbo/Garbo.cdr/pc/source/star.lzh/star.18

Text File | 1990-04-06 | 31.9 KB | 1,103 lines

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 18 (of 32)."
# Contents:  starchart/ssup.c.ab
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'starchart/ssup.c.ab' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'starchart/ssup.c.ab'\"
else
echo shar: Extracting \"'starchart/ssup.c.ab'\" \(30231 characters\)
sed "s/^X//" >'starchart/ssup.c.ab' <<'END_OF_FILE'
X  double xroot1, xroot2;
X  double yr1, yr2, r1, r2;
X  int n;
X  int x1, y1;
X  double lat_1, lon_1;
X  int in;
X
X  if (fabs(x_2 - x_1) < 1e-5) {        /* Line has infinite slope */
X    *x = x_1;
X    *y = a*x_1*x_1 + b;
X    *int_1 = TRUE;
X  } else {                /* Line slope may be calculated */
X    c = (y_2 - y_1)/(x_2 - x_1);
X    d = y_1 - c * x_1;
X    if (a < 1e-5) {            /* virtually a straight line y = b */
X      if (fabs(c) < 1e-5) {        /* Constant y */
X    n = 0;
X      } else {
X    xroot1 = (b - d)/c;
X    n = 1;
X      };
X    } else {
X      quadrat(a, -c, b-d, &xroot1, &xroot2, &n);
X    };
X    if (n == 0) {            /* No intersection */
X      *int_1 = FALSE;
X    } else if (n == 1) {        /* One intersection */
X      *x = xroot1;
X      *y = a*xroot1*xroot1 + b;
X      *int_1 = TRUE;
X    } else {                /* Two intersections */
X      yr1 = c*xroot1 + d;
X      yr2 = c*xroot2 + d;
X      r1 = (xroot1-x_1)*(xroot1-x_1) + (yr1-y_1)*(yr1-y_1)
X        + (xroot1-x_2)*(xroot1-x_2) + (yr1-y_2)*(yr1-y_2);
X      r2 = (xroot2-x_1)*(xroot2-x_1) + (yr2-y_1)*(yr2-y_1)
X        + (xroot2-x_2)*(xroot2-x_2) + (yr2-y_2)*(yr2-y_2);
X      if (r1 > r2) {
X    *x = xroot2;
X    *y = yr2;
X    *int_1 = TRUE;
X      } else {
X    *x = xroot1;
X    *y = yr1;
X    *int_1 = TRUE;
X      };
X    }
X  }
X
X  if (*int_1)
X    if (((((y_1-Fuz) <= *y) && (*y <= (y_2+Fuz)))
X     || (((y_2-Fuz) <= *y) && (*y <= (y_1+Fuz))))
X    && ((((x_1-Fuz) <= *x) && (*x <= (x_2+Fuz)))
X        || (((x_2-Fuz) <= *x) && (*x <= (x_1+Fuz)))))
X      *int_1 = TRUE;
X    else
X      *int_1 = FALSE;
X
X  if (*int_1) {
X    inv_gt(*x, *y, &lat_1, &lon_1);
X    xform(lat_1, lon_1, &x1, &y1, &in);
X    if (!in) *int_1 = FALSE;
X  }
X}
X
X
Xellip_intersect(x_1, y_1, x_2, y_2, a, b, y0, x, y, int_1)
X     double x_1, y_1, x_2, y_2, *x, *y;
X     double a, b, y0;        /* x*x/a + (y-y0)*(y-y0)/b - 1 = 0 */
X     int *int_1;
X{
X  double c, d;
X  double xroot1, xroot2;
X  double yr1, yr2, r1, r2;
X  int n;
X  int x1, y1;
X  double lat_1, lon_1;
X  int in;
X
X  if (fabs(x_2 - x_1) < 1e-5) {        /* Line has infinite slope */
X    xroot1 = xroot2 = *x = x_1;
X    if (x_1*x_1/a > 1.0) {
X      n = 0;
X      *int_1 = FALSE;
X    } else if (x_1*x_1/a == 1.0) {
X      yr1 = y0;
X      n = 1;
X      *int_1 = TRUE;
X    } else {
X      yr1 = y0 + sqrt(b*(1 - x_1*x_1/a));
X      yr2 = y0 - sqrt(b*(1 - x_1*x_1/a));
X      n = 2;
X      *int_1 = TRUE;
X    }
X  } else {                /* Line slope may be calculated */
X    c = (y_2 - y_1)/(x_2 - x_1);
X    d = y_1 - c * x_1;
X    quadrat(b + a*c*c, 2*a*c*(d - y0), a*(d - y0)*(d - y0) - a*b,
X        &xroot1, &xroot2, &n);
X    if (n == 0) {        /* No intersection */
X      *int_1 = FALSE;
X    } else if (n == 1) {    /* One intersection */
X      *x = xroot1;
X      *y = c*xroot1 + d;
X      *int_1 = TRUE;
X    } else {            /* Two intersections */
X      yr1 = c*xroot1 + d;
X      yr2 = c*xroot2 + d;
X      *int_1 = TRUE;
X    };
X  };
X
X
X  if (n == 2) {
X    r1 = (xroot1-x_1)*(xroot1-x_1) + (yr1-y_1)*(yr1-y_1)
X          + (xroot1-x_2)*(xroot1-x_2) + (yr1-y_2)*(yr1-y_2);
X    r2 = (xroot2-x_1)*(xroot2-x_1) + (yr2-y_1)*(yr2-y_1)
X          + (xroot2-x_2)*(xroot2-x_2) + (yr2-y_2)*(yr2-y_2);
X    if (r1 > r2) {
X      *x = xroot2;
X      *y = yr2;
X      *int_1 = TRUE;
X    } else {
X      *x = xroot1;
X      *y = yr1;
X      *int_1 = TRUE;
X    }
X  }
X
X  if (*int_1)
X    if ((((y_1 <= *y) && (*y <= y_2)) || ((y_2 <= *y) && (*y <= y_1)))
X    && (((x_1 <= *x) && (*x <= x_2)) || ((x_2 <= *x) && (*x <= x_1))))
X      *int_1 = TRUE;
X    else
X      *int_1 = FALSE;
X
X  if (*int_1) {
X    inv_gt(*x, *y, &lat_1, &lon_1);
X    xform(lat_1, lon_1, &x1, &y1, &in);
X    if (!in) *int_1 = FALSE;
X  }
X}
X
X
Xcirc_intersect(x_1, y_1, x_2, y_2, r, x1, y1, int_1, x2, y2, int_2)
X     double x_1, y_1, x_2, y_2, *x1, *y1, *x2, *y2;
X     double r;
X     int *int_1, *int_2;
X{
X  double c, d;
X  double xroot1, xroot2;
X  double yr1, yr2, r1, r2;
X  int n;
X  int xt1, yt1;
X  double lat_1, lon_1;
X  int in;
X
X  if (fabs(x_2 - x_1) < 1e-5) {        /* Line has infinite slope */
X    xroot1 = xroot2 = *x1 = *x2 = x_1;
X    if (fabs(r) > fabs(x_1)) {
X      yr1 = sqrt(r*r - x_1*x_1);
X      yr2 = -yr1;
X      n = 2;
X      *int_1 = *int_2 = TRUE;
X    } else if (fabs(r) == fabs(x_1)) {
X      yr1 = 0;
X      n = 1;
X      *int_1 = *int_2 = TRUE;
X    } else {
X      n = 0;
X      *int_1 = *int_2 = FALSE;
X    };
X  } else {                /* Line slope may be calculated */
X    c = (y_2 - y_1)/(x_2 - x_1);
X    d = y_1 - c * x_1;
X    quadrat(1 + c*c, 2*c*d, d*d - r*r,
X        &xroot1, &xroot2, &n);
X    if (n == 0) {        /* No intersection */
X      *int_1 = *int_2 = FALSE;
X    } else if (n == 1) {    /* One intersection */
X      *x1 = xroot1;
X      *y1 = c*xroot1 + d;
X      *int_1 = TRUE;
X      *int_2 = FALSE;
X    } else {            /* Two intersections */
X      yr1 = c*xroot1 + d;
X      yr2 = c*xroot2 + d;
X      *int_1 = *int_2 = TRUE;
X    };
X  };
X
X  if (n == 2) {
X    r1 = (xroot1-x_1)*(xroot1-x_1) + (yr1-y_1)*(yr1-y_1)
X          + (xroot1-x_2)*(xroot1-x_2) + (yr1-y_2)*(yr1-y_2);
X    r2 = (xroot2-x_1)*(xroot2-x_1) + (yr2-y_1)*(yr2-y_1)
X          + (xroot2-x_2)*(xroot2-x_2) + (yr2-y_2)*(yr2-y_2);
X    if (r1 > r2) {
X      *x1 = xroot2;
X      *y1 = yr2;
X      *x2 = xroot1;
X      *y2 = yr1;
X      *int_1 = *int_2 = TRUE;
X    } else {
X      *x1 = xroot1;
X      *y1 = yr1;
X      *x2 = xroot2;
X      *y2 = yr2;
X      *int_1 = *int_2 = TRUE;
X    }
X  }
X
X  if (*int_1)
X    if ((((y_1 <= *y1) && (*y1 <= y_2)) || ((y_2 <= *y1) && (*y1 <= y_1)))
X    && (((x_1 <= *x1) && (*x1 <= x_2)) || ((x_2 <= *x1) && (*x1 <= x_1))))
X      *int_1 = TRUE;
X    else
X      *int_1 = FALSE;
X
X  if (*int_1) {
X    inv_gt(*x1, *y1, &lat_1, &lon_1);
X    xform(lat_1, lon_1, &xt1, &yt1, &in);
X    if (!in) *int_1 = FALSE;
X  }
X
X  if (*int_2)
X    if ((((y_1 <= *y2) && (*y2 <= y_2)) || ((y_2 <= *y2) && (*y2 <= y_1)))
X    && (((x_1 <= *x2) && (*x2 <= x_2)) || ((x_2 <= *x2) && (*x2 <= x_1))))
X      *int_2 = TRUE;
X    else
X      *int_2 = FALSE;
X
X  if (*int_2) {
X    inv_gt(*x2, *y2, &lat_1, &lon_1);
X    xform(lat_1, lon_1, &xt1, &yt1, &in);
X    if (!in) *int_2 = FALSE;
X  }
X}
X
X
Xquadrat(a, b, c, x_1, x_2, n)
X     double a, b, c, *x_1, *x_2;
X     int *n;
X{
X  double t;
X
X  if (a == 0) {
X    *n = 0;
X  } else {
X    t = b * b - 4 * a * c;
X    if (t < 0) {
X      *n = 0;
X    } else if (t == 0) {
X      *x_1 = -b/(2*a);
X      *n = 1;
X    } else {
X      *x_1 = (-b + sqrt(t))/(2*a);
X      *x_2 = (-b - sqrt(t))/(2*a);
X      *n = 2;
X    };
X  };
X}
X
X/* Draw a curved line between points 1 and 2.
Xclipxform has been called, xloc1, yloc1, xloc2, yloc2 are in bounds */
X
Xdrawcurveline(lat1, lon1, lat2, lon2,
X          xloc1, yloc1, xloc2, yloc2, line_style, great_circle, clevel)
X     int xloc1, yloc1, xloc2, yloc2;
X     double lat1, lon1, lat2, lon2;
X     int line_style;
X     int great_circle;     /* draw as great circle */
X     int clevel;         /* safety valve */
X{
X  double mlat, mlon;    /* midpoint lat and long */
X  double tlat1, tlat2;    /* temporary */
X  int txloc1, tyloc1, in;
X  int txloc2, tyloc2;
X  double slope1;
X  int mxloc, myloc;    /* transformed */
X  int mpx, mpy;        /* from given x,y */
X  int inregion;
X
X/* ra difference should be less than 180 degrees: take shortest path */
X  if ((xfs_proj_mode == STEREOGR) || (xfs_proj_mode == GNOMONIC)
X      || (xfs_proj_mode == ORTHOGR))
X    if ((lon1 - lon2) > 180.0) lon1 -= 360.0;
X    else if ((lon2 - lon1) > 180.0) lon2 -= 360.0;
X
X/* Adjust so lon1 and lon2 are continuous across region: see xform */
X/* needed so midpoint is correct */
X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
X    if ((xf_west < 0.0) && (lon1>(xf_west+360.0))) lon1 -= 360.0;
X    if ((xf_east > 360.0) && (lon1<(xf_east-360.0))) lon1 += 360.0;
X    if ((xf_west < 0.0) && (lon2>(xf_west+360.0))) lon2 -= 360.0;
X    if ((xf_east > 360.0) && (lon2<(xf_east-360.0))) lon2 += 360.0;
X
X    /* path crosses boundary of map */
X    /* if the distance from point1 to left (East) + point2 to right
X       or point2 to right and point2 to left
X       is less than distance from point1 to point2,
X       then we have a problem. */
X    if (xfs_wide_warn)
X      if ((fabs(lon1-xf_east) + fabs(xf_west-lon2)) < (fabs(lon1-lon2))) {
X    slope1 = (lat2 - lat1)/(fabs(lon1-xf_east) + fabs(xf_west-lon2));
X
X    tlat1 = lat1 + slope1 * fabs(lon1-(xf_east-xf_c_scale));
X    xform(tlat1, xf_east-xf_c_scale, &txloc1, &tyloc1, &in);
X    drawcurveline(lat1, lon1, tlat1, xf_east-xf_c_scale, xloc1, yloc1,
X              txloc1, tyloc1, line_style, great_circle, 0);
X
X    tlat2 = lat2 - slope1 * fabs(lon2-(xf_west+xf_c_scale));
X    xform(tlat2, xf_west+xf_c_scale, &txloc2, &tyloc2, &in);
X    drawcurveline(tlat2, xf_west+xf_c_scale, lat2, lon2, txloc2, tyloc2,
X              xloc2, yloc2, line_style, great_circle, 0);
X    return;
X      } else if ((fabs(lon2-xf_east)+fabs(xf_west-lon1)) < (fabs(lon1-lon2))) {
X    slope1 = (lat1 - lat2)/(fabs(lon2-xf_east) + fabs(xf_west-lon1));
X
X    tlat1 = lat2 + slope1 * fabs(lon2-(xf_east-xf_c_scale));
X    xform(tlat1, xf_east-xf_c_scale, &txloc1, &tyloc1, &in);
X    drawcurveline(lat2, lon2, tlat1, xf_east-xf_c_scale, xloc2, yloc2,
X              txloc1, tyloc1, line_style, great_circle, 0);
X
X    tlat2 = lat1 - slope1 * fabs(lon1-(xf_west+xf_c_scale));
X    xform(tlat2, xf_west+xf_c_scale, &txloc2, &tyloc2, &in);
X    drawcurveline(tlat2, xf_west+xf_c_scale, lat1, lon1, txloc2, tyloc2,
X              xloc1, yloc1, line_style, great_circle, 0);
X    return;
X      };
X  }
X
X
X  if (great_circle) {
X    gcmidpoint(lat1, lon1, lat2, lon2, &mlat, &mlon);
X  } else {
X    mlat = (lat1 + lat2)/2;
X    mlon = (lon1 + lon2)/2;
X  };
X
X
X
X  xform(mlat, mlon, &mxloc, &myloc, &inregion);
X
X  mpx = (xloc1 + xloc2)/2;
X  mpy = (yloc1 + yloc2)/2;
X
X  if (((abs(mpx - mxloc) + abs(mpy - myloc)) > 0) && (clevel < 100)) {
X    /* center is not where it should be */
X    drawcurveline(lat1, lon1, mlat, mlon, xloc1, yloc1,
X          mxloc, myloc, line_style, great_circle, ++clevel);
X    drawcurveline(mlat, mlon, lat2, lon2, mxloc, myloc,
X          xloc2, yloc2, line_style, great_circle, ++clevel);
X  } else {
X    D_movedraw(xloc1, yloc1, xloc2, yloc2, line_style);
X  }
X}
X
Xgcmidpoint(lat1, lon1, lat2, lon2, pmlat, pmlon)
X     double lat1, lon1, lat2, lon2, *pmlat, *pmlon;
X{
X  double l1, m1, n1;
X  double l2, m2, n2;
X  double l3, m3, n3;
X
X  /* transform from (ra, dec) to l m n direction cosines */
X  l1 = DCOS(lat1)*DCOS(lon1);
X  m1 = DCOS(lat1)*DSIN(lon1);
X  n1 = DSIN(lat1);
X
X  l2 = DCOS(lat2)*DCOS(lon2);
X  m2 = DCOS(lat2)*DSIN(lon2);
X  n2 = DSIN(lat2);
X
X  l3 = l1 + l2;
X  m3 = m1 + m2;
X  n3 = n1 + n2;
X  n3 /= sqrt(l3*l3 + m3*m3 + n3*n3);
X
X  *pmlon = atan2(m3, l3) / 0.0174532925199;
X  if ((*pmlon < 0) && (lon1 > 0) && (lon2 > 0)) *pmlon += 360.0;
X  *pmlat = asin(n3) / 0.0174532925199;
X}
X
X
X
X/* Area handling */
X/* Works often, but is far from rigorous */
Xstatic struct t_area_struct {
X  double lat, lon;
X  int great_circle;
X  int clipped_at;
X} area_raw[100], area_clip[200], area_1[200], area_2[200];
Xstatic int nsegments = 0, nseg_clip = 0, nseg_1 = 0, nseg_2 = 0;
X
Xareastart(lat, lon, great_circle)
X     double lat, lon;
X     int great_circle;
X{
X  if (nsegments != 0)
X    areafinish();
X
X  nsegments = 0;
X  area_raw[nsegments].lat = lat;
X  area_raw[nsegments].lon = lon;
X  area_raw[nsegments].great_circle = great_circle;
X  nsegments++;
X}
X
X
Xareaadd(lat, lon, great_circle)
X     double lat, lon;
X     int great_circle;
X{
X  area_raw[nsegments].lat = lat;
X  area_raw[nsegments].lon = lon;
X  area_raw[nsegments].great_circle = great_circle;
X  nsegments++;
X}
X
Xint wrap_area();
X
Xareafinish()
X{
X  int i;
X  int xloc, yloc, xloc2, yloc2;
X  int inregion;
X  double tlat1, tlon1, tlat2, tlon2, tlat3, tlon3, tlat4, tlon4;
X  double slope1;
X  int done, old_n;
X  int curr_area;
X  int gt12;
X  int wraps;
X  int done_clip;
X  done_clip = FALSE;
X  nseg_clip = nseg_1 = nseg_2 = 0;
X  wraps = FALSE;
X
X
X  /* skip if no point is within region */
X  done = TRUE;
X  for (i = 0; i < nsegments; i++) {
X    xform(area_raw[i].lat, area_raw[i].lon, &xloc, &yloc, &inregion);
X    if (inregion) done = FALSE;
X  };
X  if (done) return;
X
X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
X    if (xfs_wide_warn)
X      for (i = 0; i < (nsegments - 1); i++)
X    if (wrap_area(area_raw[i].lat, area_raw[i].lon,
X              area_raw[i+1].lat, area_raw[i+1].lon,
X              &tlat1, &tlon1, &tlat2, &tlon2))
X      wraps = TRUE;
X
X    if (!xfs_wide_warn) {
X      /* Use algorithm to clip area_raw to rectangular region */
X      /* Use area_1 and area_2 as scratch */
X
X      /* clip against west from area_raw to area_1 */
X      tlat1 = area_raw[nsegments-1].lat;
X      tlon1 = area_raw[nsegments-1].lon;
X      for (i = 0, nseg_1 = 0; i < nsegments; i++) {
X    tlat4 = area_raw[i].lat;
X    tlon4 = area_raw[i].lon;
X
X    if (eastof(tlon4, xf_west)) {
X      if (eastof(tlon1, xf_west)) {
X        /* last point and new point (4) are in, add new point */
X      } else {
X        /* new point in, but old out.
X         Calculate intersection point (2), add it and new point */
X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
X
X        tlat2 = tlat1 + slope1 * (xf_west-tlon1);
X        tlon2 = xf_west;
X        area_1[nseg_1].lat = tlat2;
X        area_1[nseg_1].lon = tlon2;
X        area_1[nseg_1].great_circle = area_raw[i].great_circle;
X        area_1[nseg_1].clipped_at = WEST_CLIP;
X        nseg_1++;
X      }
X      area_1[nseg_1].lat = tlat4;
X      area_1[nseg_1].lon = tlon4;
X      area_1[nseg_1].great_circle = area_raw[i].great_circle;
X      area_1[nseg_1].clipped_at = 0;
X      nseg_1++;
X    } else {
X      if (eastof(tlon1, xf_west)) {
X        /* new point out, but old point in.
X           Calculate intersection and add it */
X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
X
X        tlat2 = tlat1 + slope1 * (xf_west-tlon1);
X        tlon2 = xf_west;
X        area_1[nseg_1].lat = tlat2;
X        area_1[nseg_1].lon = tlon2;
X        area_1[nseg_1].great_circle = area_raw[i].great_circle;
X        area_1[nseg_1].clipped_at = WEST_CLIP;
X        nseg_1++;
X      }
X    }
X    tlat1 = tlat4;
X    tlon1 = tlon4;
X      }
X
X
X      /* clip against south from area_1 to area_2 */
X      tlat1 = area_1[nseg_1-1].lat;
X      tlon1 = area_1[nseg_1-1].lon;
X      for (i = 0, nseg_2 = 0; i < nseg_1; i++) {
X    tlat4 = area_1[i].lat;
X    tlon4 = area_1[i].lon;
X
X    if (tlat4 > xf_south) {
X      if (tlat1 > xf_south) {
X        /* last point and new point (4) are in, add new point */
X      } else {
X        /* new point in, but old out.
X         Calculate intersection point (2), add it and new point */
X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
X
X        tlat2 = xf_south;
X        tlon2 = tlon1 + slope1 * (xf_south-tlat1);
X        area_2[nseg_2].lat = tlat2;
X        area_2[nseg_2].lon = tlon2;
X        area_2[nseg_2].great_circle = area_1[i].great_circle;
X        area_2[nseg_2].clipped_at = SOUTH_CLIP;
X        nseg_2++;
X      }
X      area_2[nseg_2].lat = tlat4;
X      area_2[nseg_2].lon = tlon4;
X      area_2[nseg_2].great_circle = area_1[i].great_circle;
X      area_2[nseg_2].clipped_at = 0;
X      nseg_2++;
X    } else {
X      if (tlat1 > xf_south) {
X        /* new point out, but old point in.
X           Calculate intersection and add it */
X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
X
X        tlat2 = xf_south;
X        tlon2 = tlon1 + slope1 * (xf_south-tlat1);
X        area_2[nseg_2].lat = tlat2;
X        area_2[nseg_2].lon = tlon2;
X        area_2[nseg_2].great_circle = area_1[i].great_circle;
X        area_2[nseg_2].clipped_at = SOUTH_CLIP;
X        nseg_2++;
X      }
X    }
X    tlat1 = tlat4;
X    tlon1 = tlon4;
X      }
X
X      /* clip against east from area_2 to area_1 */
X      tlat1 = area_2[nseg_2-1].lat;
X      tlon1 = area_2[nseg_2-1].lon;
X      for (i = 0, nseg_1 = 0; i < nseg_2; i++) {
X    tlat4 = area_2[i].lat;
X    tlon4 = area_2[i].lon;
X
X    if (westof(tlon4, xf_east)) {
X      if (westof(tlon1, xf_east)) {
X        /* last point and new point (4) are in, add new point */
X      } else {
X        /* new point in, but old out.
X         Calculate intersection point (2), add it and new point */
X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
X
X        tlat2 = tlat1 + slope1 * (xf_east-tlon1);
X        tlon2 = xf_east;
X        area_1[nseg_1].lat = tlat2;
X        area_1[nseg_1].lon = tlon2;
X        area_1[nseg_1].great_circle = area_2[i].great_circle;
X        area_1[nseg_1].clipped_at = EAST_CLIP;
X        nseg_1++;
X      }
X      area_1[nseg_1].lat = tlat4;
X      area_1[nseg_1].lon = tlon4;
X      area_1[nseg_1].great_circle = area_2[i].great_circle;
X      area_1[nseg_1].clipped_at = 0;
X      nseg_1++;
X    } else {
X      if (westof(tlon1, xf_east)) {
X        /* new point out, but old point in.
X           Calculate intersection and add it */
X        slope1 = (tlat4 - tlat1) / (tlon4 - tlon1);
X
X        tlat2 = tlat1 + slope1 * (xf_east-tlon1);
X        tlon2 = xf_east;
X        area_1[nseg_1].lat = tlat2;
X        area_1[nseg_1].lon = tlon2;
X        area_1[nseg_1].great_circle = area_2[i].great_circle;
X        area_1[nseg_1].clipped_at = EAST_CLIP;
X        nseg_1++;
X      }
X    }
X    tlat1 = tlat4;
X    tlon1 = tlon4;
X      }
X
X      /* clip against north from area_1 to area_clip */
X      tlat1 = area_1[nseg_1-1].lat;
X      tlon1 = area_1[nseg_1-1].lon;
X      for (i = 0, nseg_clip = 0; i < nseg_1; i++) {
X    tlat4 = area_1[i].lat;
X    tlon4 = area_1[i].lon;
X
X    if (tlat4 < xf_north) {
X      if (tlat1 < xf_north) {
X        /* last point and new point (4) are in, add new point */
X      } else {
X        /* new point in, but old out.
X         Calculate intersection point (2), add it and new point */
X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
X
X        tlat2 = xf_north;
X        tlon2 = tlon1 + slope1 * (xf_north-tlat1);
X        area_clip[nseg_clip].lat = tlat2;
X        area_clip[nseg_clip].lon = tlon2;
X        area_clip[nseg_clip].great_circle = area_1[i].great_circle;
X        area_clip[nseg_clip].clipped_at = NORTH_CLIP;
X        nseg_clip++;
X      }
X      area_clip[nseg_clip].lat = tlat4;
X      area_clip[nseg_clip].lon = tlon4;
X      area_clip[nseg_clip].great_circle = area_1[i].great_circle;
X      area_clip[nseg_clip].clipped_at = 0;
X      nseg_clip++;
X    } else {
X      if (tlat1 < xf_north) {
X        /* new point out, but old point in.
X           Calculate intersection and add it */
X        slope1 = (tlon4 - tlon1) / (tlat4 - tlat1);
X
X        tlat2 = xf_north;
X        tlon2 = tlon1 + slope1 * (xf_north-tlat1);
X        area_clip[nseg_clip].lat = tlat2;
X        area_clip[nseg_clip].lon = tlon2;
X        area_clip[nseg_clip].great_circle = area_1[i].great_circle;
X        area_clip[nseg_clip].clipped_at = NORTH_CLIP;
X        nseg_clip++;
X      }
X    }
X    tlat1 = tlat4;
X    tlon1 = tlon4;
X      }
X
X      if ((area_clip[0].lat != area_clip[nseg_clip-1].lat)
X      || (area_clip[0].lon != area_clip[nseg_clip-1].lon)) {
X    area_clip[nseg_clip].lat = area_clip[0].lat;
X    area_clip[nseg_clip].lon = area_clip[0].lon;
X    area_clip[nseg_clip].great_circle = area_raw[i].great_circle;
X    area_clip[nseg_clip].clipped_at = WEST_CLIP;
X    nseg_clip++;
X      }
X
X      if (nseg_clip > 0) doarea(area_clip, nseg_clip);
X
X      done_clip = TRUE;
X    };
X  };
X
X  if (!done_clip) {
X    /* Scan list of segments, constructing clipped polygon */
X    for (i = 0; i < (nsegments - 1); i++) {
X      if (clipr_xform(area_raw[i].lat, area_raw[i].lon,
X              area_raw[i+1].lat, area_raw[i+1].lon,
X              &xloc, &yloc, &xloc2, &yloc2, area_raw[i+1].great_circle,
X              &tlat1, &tlon1, &tlat2, &tlon2)) {
X    addsegment(tlat1, tlon1, tlat2, tlon2,
X           (nseg_clip > 0) ? area_clip[nseg_clip-1].clipped_at : 0,
X           clip_at1, clip_at2,
X           area_raw[i].great_circle, area_raw[i+1].great_circle);
X      };
X    };
X
X    /* ensure that area is closed. */
X    /* redraw first segment */
X    if ((nseg_clip > 0) && ((area_clip[nseg_clip-1].lat != area_clip[0].lat)
X               || (area_clip[nseg_clip-1].lon != area_clip[0].lon))) {
X      i = 0;
X      done = FALSE;
X      old_n = nseg_clip;
X      do {
X    if (clipr_xform(area_raw[i].lat, area_raw[i].lon,
X            area_raw[i+1].lat, area_raw[i+1].lon,
X            &xloc, &yloc, &xloc2, &yloc2,
X            area_raw[i+1].great_circle,
X            &tlat1, &tlon1, &tlat2, &tlon2)) {
X      addsegment(tlat1, tlon1, tlat2, tlon2,
X             (nseg_clip > 0) ? area_clip[nseg_clip-1].clipped_at : 0,
X             clip_at1, clip_at2,
X             area_raw[i].great_circle, area_raw[i+1].great_circle);
X      done = TRUE;
X    }
X    i++;
X      } while ((i < (nsegments - 1)) && (!done));
X      if ((old_n+1) < nseg_clip) nseg_clip--; /* only needed one added */
X    }
X
X
X    /* Now fix problem of areas crossing from left to right across boundary */
X    /* if the distance from point1 to left (East) + point2 to right
X       or point2 to right and point2 to left
X       is less than distance from point1 to point2,
X       then we have a problem. */
X    curr_area = 1;
X    for (i = 0; i < (nseg_clip - 1); i++) {
X      if (wraps && wrap_area(area_clip[i].lat, area_clip[i].lon,
X                 area_clip[i+1].lat, area_clip[i+1].lon,
X                 &tlat1, &tlon1, &tlat2, &tlon2)) {
X    /* Crosses boundary */
X    switch (curr_area) {
X    case 1:
X                /* finish segment in part 1,
X                   begin part 2 */
X      area_1[nseg_1].lat = area_clip[i].lat;
X      area_1[nseg_1].lon = area_clip[i].lon;
X      area_1[nseg_1].great_circle = area_clip[i].great_circle;
X      area_1[nseg_1].clipped_at = area_clip[i].clipped_at;
X      nseg_1++;
X
X      area_1[nseg_1].lat = tlat1;
X      area_1[nseg_1].lon = tlon1;
X      area_1[nseg_1].great_circle = area_clip[i].great_circle;
X      area_1[nseg_1].clipped_at = area_clip[i].clipped_at;
X      nseg_1++;
X
X      curr_area = 2;
X
X      area_2[nseg_2].lat = tlat2;
X      area_2[nseg_2].lon = tlon2;
X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
X      nseg_2++;
X
X      gt12 = (tlon1 > tlon2); /* area 1 greater than area 2 */
X      break;
X    case 2:
X                /* finish part 2, draw it,
X                   close gap in part 1 */
X      if (((gt12) && (tlon1 < tlon2))
X          || ((!gt12) && (tlon1 > tlon2))){
X        tlat3 = tlat1;
X        tlon3 = tlon1;
X        tlat1 = tlat2;
X        tlon1 = tlon2;
X        tlat2 = tlat3;
X        tlon2 = tlon3;
X      };
X      area_2[nseg_2].lat = area_clip[i].lat;
X      area_2[nseg_2].lon = area_clip[i].lon;
X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
X      nseg_2++;
X
X      area_2[nseg_2].lat = tlat2;
X      area_2[nseg_2].lon = tlon2;
X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
X      nseg_2++;
X
X      area_2[nseg_2].lat = area_2[0].lat;
X      area_2[nseg_2].lon = area_2[0].lon;
X      area_2[nseg_2].great_circle = TRUE;
X      area_2[nseg_2].clipped_at = area_2[0].clipped_at;
X      nseg_2++;
X
X      doarea(area_2, nseg_2);
X
X      curr_area = 1;
X
X      area_1[nseg_1].lat = tlat1;
X      area_1[nseg_1].lon = tlon1;
X      area_1[nseg_1].great_circle = TRUE;
X      area_1[nseg_1].clipped_at = area_1[nseg_1-1].clipped_at;
X      nseg_1++;
X
X      break;
X    default:
X      break;
X    };
X      } else {
X    switch (curr_area) {
X    case 1:
X      area_1[nseg_1].lat = area_clip[i].lat;
X      area_1[nseg_1].lon = area_clip[i].lon;
X      area_1[nseg_1].great_circle = area_clip[i].great_circle;
X      area_1[nseg_1].clipped_at = area_clip[i].clipped_at;
X      nseg_1++;
X      break;
X    case 2:
X      area_2[nseg_2].lat = area_clip[i].lat;
X      area_2[nseg_2].lon = area_clip[i].lon;
X      area_2[nseg_2].great_circle = area_clip[i].great_circle;
X      area_2[nseg_2].clipped_at = area_clip[i].clipped_at;
X      nseg_2++;
X      break;
X    default:
X      break;
X    }
X      }
X    };
X    if (nseg_clip > 0) {
X      area_1[nseg_1].lat = area_clip[nseg_clip-1].lat;
X      area_1[nseg_1].lon = area_clip[nseg_clip-1].lon;
X      area_1[nseg_1].great_circle = area_clip[nseg_clip-1].great_circle;
X      area_1[nseg_1].clipped_at = area_clip[nseg_clip-1].clipped_at;
X      nseg_1++;
X    };
X
X    if (nseg_1 > 0) doarea(area_1, nseg_1);
X  };
X
X  /* Mark it done */
X  nsegments = 0;
X}
X
Xdoarea(area, nseg)
X     struct t_area_struct area[];
X     int nseg;
X{
X  int i;
X  int xloc, yloc, xloc2, yloc2;
X  double tlat1, tlon1, tlat2, tlon2;
X
X  for (i = 0; i < (nseg - 1); i++) {
X    if (clipr_xform(area[i].lat, area[i].lon,
X            area[i+1].lat, area[i+1].lon,
X            &xloc, &yloc, &xloc2, &yloc2, area[i+1].great_circle,
X            &tlat1, &tlon1, &tlat2, &tlon2)) {
X
X      if (i == 0) D_areamove(xloc, yloc);
X      /* Here we can add tests to add segments along the boundary of
X     projection modes other than SANSONS and RECTANGULAR */
X
X      addareacurved(tlat1, tlon1, tlat2, tlon2,
X            xloc, yloc, xloc2, yloc2,
X            area[i+1].great_circle, 0);
X      if (i == (nseg - 2)) D_areafill(xloc2, yloc2);
X    }
X  }
X}
X
Xint wrap_area(lat1, lon1, lat2, lon2, plat1, plon1, plat2, plon2)
X     double lat1, lon1, lat2, lon2;
X     double *plat1, *plon1, *plat2, *plon2;
X{
X  double slope1, tlat1, tlat2;
X
X/* Adjust so lon1 and lon2 are continuous across region: see xform */
X/* needed so midpoint is correct */
X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
X    if ((xf_west < 0.0) && (lon1>(xf_west+360.0))) lon1 -= 360.0;
X    if ((xf_east > 360.0) && (lon1<(xf_east-360.0))) lon1 += 360.0;
X    if ((xf_west < 0.0) && (lon2>(xf_west+360.0))) lon2 -= 360.0;
X    if ((xf_east > 360.0) && (lon2<(xf_east-360.0))) lon2 += 360.0;
X
X    /* path crosses boundary of map */
X    /* if the distance from point1 to left (East) + point2 to right
X       or point2 to right and point2 to left
X       is less than distance from point1 to point2,
X       then we have a problem. */
X    if (xfs_wide_warn)
X      if ((fabs(lon1-xf_east) + fabs(xf_west-lon2)) < (fabs(lon1-lon2))) {
X    /* point 1 close to east boundary */
X    slope1 = (lat2 - lat1)/(fabs(lon1-xf_east) + fabs(xf_west-lon2));
X
X    tlat1 = lat1 + slope1 * fabs(lon1-(xf_east-xf_c_scale));
X
X    /* segment from (lat1, lon1) to (tlat1, xf_east-xf_c_scale)
X       is on one side of boundary */
X    *plat1 = tlat1;
X    *plon1 = xf_east-xf_c_scale;
X    if (*plon1 > 360) *plon1 -= 360;
X    if (*plon1 < 0) *plon1 += 360;
X
X    tlat2 = lat2 - slope1 * fabs(lon2-(xf_west+xf_c_scale));
X    /* segment from (tlat2, xf_west+xf_c_scale) to (lat2, lon2)
X       is on one side of boundary */
X    *plat2 = tlat2;
X    *plon2 = xf_west+xf_c_scale;
X    if (*plon2 > 360) *plon2 -= 360;
X    if (*plon2 < 0) *plon2 += 360;
X    return TRUE;
X      } else if ((fabs(lon2-xf_east)+fabs(xf_west-lon1)) < (fabs(lon1-lon2))) {
X    /* point 1 close to west boundary */
X    slope1 = (lat1 - lat2)/(fabs(lon2-xf_east) + fabs(xf_west-lon1));
X
X    tlat1 = lat2 + slope1 * fabs(lon2-(xf_east-xf_c_scale));
X
X    /* segment from (lat2, lon2) to (tlat1, xf_east-xf_c_scale)
X       is on one side of boundary */
X    *plat2 = tlat1;
X    *plon2 = xf_east-xf_c_scale;
X    if (*plon2 > 360) *plon2 -= 360;
X    if (*plon2 < 0) *plon2 += 360;
X
X    tlat2 = lat1 - slope1 * fabs(lon1-(xf_west+xf_c_scale));
X    /* segment from (tlat2, xf_west+xf_c_scale) to (lat1, lon1)
X       is on one side of boundary */
X    *plat1 = tlat2;
X    *plon1 = xf_west+xf_c_scale;
X    if (*plon1 > 360) *plon1 -= 360;
X    if (*plon1 < 0) *plon1 += 360;
X    return TRUE;
X      };
X  }
X  return FALSE;
X}
X
X
Xaddsegment(tlat1, tlon1, tlat2, tlon2, last_clip, clip_1, clip_2,
X       gc_1, gc_2)
X     double tlat1, tlon1, tlat2, tlon2;
X     int last_clip, clip_1, clip_2;
X     int gc_1, gc_2;
X{
X
X  /* just add a segment */
X  area_clip[nseg_clip].lat = tlat1;
X  area_clip[nseg_clip].lon = tlon1;
X  area_clip[nseg_clip].great_circle = gc_1;
X  area_clip[nseg_clip].clipped_at = clip_1;
X  nseg_clip++;
X  area_clip[nseg_clip].lat = tlat2;
X  area_clip[nseg_clip].lon = tlon2;
X  area_clip[nseg_clip].great_circle = gc_2;
X  area_clip[nseg_clip].clipped_at = clip_2;
X  nseg_clip++;
X}
X
X
X
X/* addareacurved: add a set of boundary segments between points 1 and 2
X   clipxform has been called, xloc1, yloc1, xloc2, yloc2 are in bounds */
Xaddareacurved(lat1, lon1, lat2, lon2,
X          xloc1, yloc1, xloc2, yloc2, great_circle, clevel)
X     int xloc1, yloc1, xloc2, yloc2;
X     double lat1, lon1, lat2, lon2;
X     int great_circle;     /* draw as great circle */
X     int clevel;         /* safety valve */
X{
X  double mlat, mlon;    /* midpoint lat and long */
X  int mxloc, myloc;    /* transformed */
X  int mpx, mpy;        /* from given x,y */
X  int inregion;
X
X/* ra difference should be less than 180 degrees: take shortest path */
X  if ((xfs_proj_mode == STEREOGR) || (xfs_proj_mode == GNOMONIC)
X      || (xfs_proj_mode == ORTHOGR))
X    if ((lon1 - lon2) > 180.0) lon1 -= 360.0;
X    else if ((lon2 - lon1) > 180.0) lon2 -= 360.0;
X
X/* Adjust so lon1 and lon2 are continuous across region: see xform */
X/* needed so midpoint is correct */
X  if ((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR)) {
X    if ((xf_west < 0.0) && (lon1>(xf_west+360.0))) lon1 -= 360.0;
X    if ((xf_east > 360.0) && (lon1<(xf_east-360.0))) lon1 += 360.0;
X    if ((xf_west < 0.0) && (lon2>(xf_west+360.0))) lon2 -= 360.0;
X    if ((xf_east > 360.0) && (lon2<(xf_east-360.0))) lon2 += 360.0;
X  }
X
X
X  if (great_circle) {
X    gcmidpoint(lat1, lon1, lat2, lon2, &mlat, &mlon);
X  } else {
X    mlat = (lat1 + lat2)/2;
X    mlon = (lon1 + lon2)/2;
X  };
X
X  xform(mlat, mlon, &mxloc, &myloc, &inregion);
X
X  mpx = (xloc1 + xloc2)/2;
X  mpy = (yloc1 + yloc2)/2;
X
X  if (((abs(mpx - mxloc) + abs(mpy - myloc)) > 0) && (clevel < 100)) {
X    /* center is not where it should be */
X    addareacurved(lat1, lon1, mlat, mlon,
X          xloc1, yloc1, mxloc, myloc, great_circle, ++clevel);
X    addareacurved(mlat, mlon, lat2, lon2,
X          mxloc, myloc, xloc2, yloc2, great_circle, ++clevel);
X  } else {
X    D_areaadd(xloc2, yloc2);
X  }
X}
X
X
X
X
X/* Translate two digit encoded size string */
X/* Maximum parsed: 89000 secs of arc = 24.666 degrees */
X/* Warning, should use 32 bit integers to allow this value */
Xlong size_obj(size_str)
X     char *size_str;
X{
X  char chr1, chr2;
X
X  chr1 = islower(size_str[0]) ? toupper(size_str[0]) : size_str[0];
X  chr2 = islower(size_str[1]) ? toupper(size_str[1]) : size_str[1];
X
X  if ((chr1 == ' ') && (chr2 == ' ')) return -1;
X  if (chr1 == ' ') return chr2 - '0';
X  if (isdigit(chr1)) return ((chr1-'0')*10L + (chr2-'0'));
X  if (chr1 < 'J') return ((chr1-'A'+1)*100L + (chr2-'0')*10L);
X  if (chr1 < 'S') return ((chr1-'I')*1000L + (chr2-'0')*100L);
X  if (chr1 <= 'Z') return ((chr1-'R')*10000L + (chr2-'0')*1000L);
X  return -1;
X}
X
X/*
XExamples:
X
X"  "        -1
X" 6"        6
X"09"        9
X"73"        73
X"A0"        100
X"C3"        330
X"D5"        450
X"I6"        960
X"J2"        1200
X"R3"        9300
X"S6"        16000
X"Z0"        80000
X"Z9"        89000
X*/
X
X/* Precession functions */
X/* precession based on formulae from
X   Astronomical Almanac, 1988, p B19
X*/
Xstatic double M, N;
Xinitprecess(eout)
X     double eout;
X{
X  double T;
X
X  T = (eout - 2000.0) / 100.0;
X  M = (1.2812323 + (0.0003879 + 0.0000101 * T) * T) * T;
X  N = (0.5567530 - (0.0001185 - 0.0000116 * T) * T) * T;
X}
X
Xdo_precess(lat, lon)
X     double *lat, *lon;
X{
X  double am, dm, a2, d2;
X
X  /* am = /alpha_m, dm = /delta_m */
X
X  am = *lon + (M + N * DSIN(*lon) * DTAN(*lat))/2.0;
X  dm = *lat + N*DCOS(am)/2.0;
X  a2 = *lon + M + N*DSIN(am)*DTAN(dm);
X  d2 = *lat + N * DCOS(am);
X
X  if (a2 >= 360.0) a2 -= 360.0;
X  if (a2 < 0.0) a2 += 360.0;
X
X  *lon = a2;
X  *lat = d2;
X}
X
END_OF_FILE
if test 30231 -ne `wc -c <'starchart/ssup.c.ab'`; then
    echo shar: \"'starchart/ssup.c.ab'\" unpacked with wrong size!
fi
# end of 'starchart/ssup.c.ab'
fi
echo shar: End of archive 18 \(of 32\).
cp /dev/null ark18isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 32 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0