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

Text File | 1990-04-06 | 32.8 KB | 1,141 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 19 (of 32)."
# Contents:  starchart/ssup.c.aa
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'starchart/ssup.c.aa' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'starchart/ssup.c.aa'\"
else
echo shar: Extracting \"'starchart/ssup.c.aa'\" \(31121 characters\)
sed "s/^X//" >'starchart/ssup.c.aa' <<'END_OF_FILE'
X/*
X * starsupp.c, extracted from starchart.c (now starmain.c)
X * starchart.c -- version 2, September 1987
X *        revision 2.1 December, 1987
X *        revision 3.0beta January, 1989
X *
X * Portions Copyright (c) 1987 by Alan Paeth (awpaeth@watcgl)
X *
X * Copyright (c) 1990 by Craig Counterman. All rights reserved.
X *
X * This software may be redistributed freely, not sold.
X * This copyright notice and disclaimer of warranty must remain
X *    unchanged.
X *
X * No representation is made about the suitability of this
X * software for any purpose.  It is provided "as is" without express or
X * implied warranty, to the extent permitted by applicable law.
X *
X */
X
X
X/*
X ! Version 2 modification authors:
X !
X !   [a] Petri Launiainen, [e-mail] (with Jyrki Yli-Nokari, [e-mail])
X !   [b] Bob Tidd, [e-mail]
X !   [c] Alan Paeth, awpaeth@watcgl
X !
X */
X
Xstatic char rcsid[]="$Header: starsupp.c,v 2.13 90/03/10 15:31:51 ccount Exp $";
X
X#include <stdio.h>
X#include <math.h>
X#ifndef SYSV
X#include <strings.h>
X#else
X#include <string.h>
X#endif
X#include <ctype.h>
X
X#include "star3.h"
X
X#ifndef READMODE
X#define READMODE "r"
X#endif
X#define OPENFAIL 0
X#define LINELEN 82
X
X
X/* PI / 180 = .0174532925199 */
X#define DCOS(x) (cos((x)*.0174532925199))
X#define DSIN(x) (sin((x)*.0174532925199))
X#define DTAN(x) (tan((x)*.0174532925199))
X#define DASIN(x) (asin(x)/.0174532925199)
X#define DATAN2(x,y) (atan2(x,y)/.0174532925199)
X#define MAX(a,b) ((a)>(b)?(a):(b))
X#define MIN(a,b) ((a)<(b)?(a):(b))
X
X/* externals from starmain.c */
Xextern mapwindow *mapwin[];
Xextern int numwins;
X
X
X/* exportable */
Xdouble xf_west, xf_east, xf_north, xf_south, xf_bottom;
Xint xf_xcen, xf_ycen, xf_ybot;
Xint xf_w_left, xf_w_right, xf_w_top, xf_w_bot;
X
Xdouble xf_c_scale;
X
X/* local storage */
Xstatic int xfs_proj_mode;
X/* sin_dlcen = sin (phi0) = sin (declination of center), cos_dlcen similar */
Xstatic double xfs_ra_cen, xfs_dl_cen, sin_dlcen, cos_dlcen, chart_scale;
Xstatic double xfs_scale;        /* Formerly yscale */
Xstatic double xfs_inv;
Xstatic int xfs_wide_warn;
X
X
X/* TRANSFORMATION FUNCTIONS */
X
X/**
X    stereographic projection:
X        Imagine object projected on a sphere of radius 1.
X        Center of chart is against a piece of paper.
X        You are looking at the center of the chart through the
X                center of the sphere.
X        The objects are seen projected on the piece of paper.
X
X    to do it:
X        1) get angle from center to object, call this theta.
X        2) get direction from center to object, call this phi.
X        3) projection of object will be 2*tan(theta/2) from the center
X                of the paper, in the direction phi.
X
X    to do steps 1 & 2, use alt-azimuth formula, with theta = 90 - alt.
X
X    scale: when theta = scale, projected object = min(ww,wh)/2 from center.
X        i.e. 2*tan(scale/2)*xfs_scale = min(ww,wh)/2.
X        or xfs_scale = (min(ww,wh)/2)/(2*tan(scale/2))
X                  = min(ww,wh)/(4*tan(scale/2))
X        put factor of 2 in xfs_scale to save a little computation
X           xfs_scale = min(ww,wh)/(2*tan(scale/2))
X        R = tan(theta/2)*xfs_scale.
X*/
X/**
X    alt-azimuth:
X      sin(alt) = sin(obs_lat)sin(obj_dl) + cos(obs_lat)cos(obj_dl)cos(hour)
X
X      cos(azi) = (cos(obs_lat)sin(obj_dl) - sin(obs_lat)cos(obj_dl)cos(hour))
X            / cos(alt)
X      sin(azi) = (-cos(obj_dl)sin(hour)) / cos(alt)
X      tan(azi) = -cos((obj_dl)sin(hour)
X              / (cos(obs_lat)sin(obj_dl) - sin(obs_lat)cos(obj_dl)cos(hour))
X
X
X    with alt = 90 - theta, azi = phi, hour = lon - racen, obj_dl =lat,
X                and obs_lat = dlcen, this becomes:
X
X    cos(theta) = sin(dlcen)sin(lat) + cos(dlcen)cos(lat)cos(lon - racen)
X
X    tan(phi) = -cos(lat)sin(lon - racen)
X              / (cos(dlcen)sin(lat) - sin(dlcen)cos(lat)cos(lon - racen))
X
X
X        racen and dlcen are constants for a chart, so introduce variables
X            racen, sin_dlcen, cos_dlcen.
X        also add chart_scale which is chart->scale in radians.
X        initxform sets these static variables.
X
X
X
X    Other projections from book.
X**/
X
X
Xinitxform(win)
X     mapwindow *win;
X{
X  double adj, xscale, tscale;
X
X  xfs_proj_mode = win->proj_mode;
X
X  if (win->scale <= 0.0) return (FALSE);
X  if (win->height == 0) return (FALSE);
X
X  xfs_ra_cen = win->racen;
X  xfs_dl_cen = win->dlcen;
X  xf_xcen = win->x_offset + (win->width)/2;
X  xf_ycen = win->y_offset + (win->height)/2;
X  xf_ybot = win->y_offset;
X
X  xf_north = (win->dlcen + win->scale / 2);
X  xf_south = (win->dlcen - win->scale / 2);
X
X  if (xf_north > 90.0) xf_north = 90.0;
X  if (xf_south < -90.0) xf_south = -90.0;
X
X  if (win->invert) {
X    xfs_inv = -1.0;
X    xf_bottom = xf_north;
X  } else {
X    xfs_inv = 1.0;
X    xf_bottom = xf_south;
X  }
X
X  if (xfs_proj_mode == SANSONS) {
X/*
X * calculate xf_east and xf_west by calculating the widest range of lon
X * which will be within the chart.
X * xscale is other than win->scale in order to widen the horizontal viewing
X * area, which otherwise shrinks near the poles under Sanson's projection
X * this happens in polar maps which do not span the celestial equator
X */
X    adj = 1.0;
X    if (xf_north * xf_south > 0.0)
X      adj = MAX(DCOS(xf_north), DCOS(xf_south));
X    xscale = win->scale/adj;
X
X    tscale = xscale*win->width/win->height / 2.0;
X    if (tscale > 180.0) tscale = 180.0;
X  } else if (xfs_proj_mode == RECTANGULAR) {
X    tscale = win->scale * win->width / (2*win->height);
X    if (tscale > 180.0) tscale = 180.0;
X  };
X
X  xf_east  = win->racen + tscale;
X  xf_west  = win->racen - tscale;
X
X  /* set warning, may have problems in SANSONS or RECTANGULAR
X     with lines which should wrap around */
X  if (((xfs_proj_mode == SANSONS) || (xfs_proj_mode == RECTANGULAR))
X      && (tscale > 90.0)) xfs_wide_warn = TRUE;
X  else xfs_wide_warn = FALSE;
X
X
X  xf_w_left = win->x_offset;
X  xf_w_right = win->x_offset + win->width;
X  xf_w_bot = win->y_offset;
X  xf_w_top = win->y_offset + win->height;
X
X  switch (xfs_proj_mode) {
X  case GNOMONIC:
X  case ORTHOGR:
X  case STEREOGR:
X    sin_dlcen = DSIN(win->dlcen);
X    cos_dlcen = DCOS(win->dlcen);
X    chart_scale = win->scale * .0174532925199; /* Radians */
X    break;
X  case SANSONS:
X  case RECTANGULAR:
X  default:
X    break;
X  }
X
X/* xf_c_scale is the size in degrees which one pixel occupies on the map */
X/* xfs_scale is the conversion factor for size of the picture
X     (= R in some formulas for stereographic,
X     gnomonic and orthographic projections) */
X  if (xfs_proj_mode == STEREOGR) {
X    xfs_scale = MIN(win->height, win->width) / (4.0*DTAN(win->scale/2.0));
X    xf_c_scale = win->c_scale = 1.0/(2.0 * DTAN(0.5) * xfs_scale);
X  } else if (xfs_proj_mode == SANSONS) {
X    xfs_scale = win->height / win->scale;
X    xf_c_scale = win->c_scale = win->scale / win->height;
X  } else if (xfs_proj_mode == GNOMONIC) {
X    xfs_scale = MIN(win->height, win->width) / (2.0 * DTAN(win->scale/2.0));
X    xf_c_scale = win->c_scale = 1.0/(DTAN(1.0) * xfs_scale);
X  } else if (xfs_proj_mode == ORTHOGR) {
X    xfs_scale = MIN(win->height, win->width) / (2.0 * DSIN(win->scale/2.0));
X    xf_c_scale = win->c_scale = 1.0/(DSIN(1.0) * xfs_scale);
X  } else if (xfs_proj_mode == RECTANGULAR) {
X    xfs_scale = win->height / win->scale;
X    xf_c_scale = win->c_scale = 1.0/ xfs_scale;
X  }
X
X  /* initialize gnomonic transform function */
X  init_gt(win);
X
X  return (TRUE);
X}
X
X
X
Xxform(lat, lon, xloc, yloc, inregion)
X     double lat, lon;
X     int *xloc, *yloc, *inregion;
X{
X  double theta, rac_l;
X  double denom;
X  double Dcoslat, Dsinlat, Dcosrac_l, Dsinrac_l;
X  /* Dcoslat, Dsinlat: of object latitude in degrees = phi
X     Dcosrac_l, Dsinrac_l: of object ra - longditude of center = d(lambda) */
X
X  switch (xfs_proj_mode) {
X  case SANSONS:
X /*
X  * This is Sanson's Sinusoidal projection. Its properties:
X  *   (1) area preserving
X  *   (2) preserves linearity along y axis (declination/azimuth)
X  */
X/* because of the (xfs_ra_cen-lon) below, lon must be continuous across the
X   plotted region.
X   xf_west is xfs_ra_cen - (scale factor),
X   xf_east is xfs_ra_cen + (scale factor),
X   so xf_west may be negative, and xf_east may be > 360.0.
X   lon should be 0 <= lon < 360.0.  we must bring lon into the range.
X   if xf_west < 0 and (xf_west + 360) < lon then lon is in the right half
X       of the chart and needs to be adjusted
X   if xf_east > 360 and (xf_east - 360) > lon then lon is in the left half
X       of the chart and needs to be adjusted
X*/
X
X
X    if ((xf_west < 0.0) && (lon>(xf_west+360.0))) lon -= 360.0;
X    if ((xf_east > 360.0) && (lon<(xf_east-360.0))) lon += 360.0;
X
X    *xloc = xf_xcen + (xfs_ra_cen-lon)*xfs_scale*DCOS(lat) + 0.5;
X    *yloc = xf_ybot + (int)xfs_inv*((lat - xf_bottom) * xfs_scale) + 0.5;
X    *inregion = ((lon >= xf_west) && (lon <= xf_east) &&
X         (lat >= xf_south) && (lat <= xf_north));
X    break;
X  case STEREOGR:
X/* stereographic projection */
X    rac_l = lon - xfs_ra_cen;
X    Dsinlat = DSIN(lat);
X    Dcoslat = DCOS(lat);
X    Dcosrac_l = DCOS(rac_l);
X    Dsinrac_l = DSIN(rac_l);
X    theta = acos(sin_dlcen*Dsinlat + cos_dlcen*Dcoslat*Dcosrac_l);
X
X    *inregion = (theta <= chart_scale);
X    if (*inregion) {
X      denom = (1+sin_dlcen*Dsinlat + cos_dlcen*Dcoslat*Dcosrac_l) / xfs_scale;
X      *xloc = xf_xcen - 2*Dcoslat*Dsinrac_l / denom + 0.5;
X      *yloc = xf_ycen + 2*xfs_inv*(cos_dlcen*Dsinlat
X                 - sin_dlcen*Dcoslat*Dcosrac_l) / denom + 0.5;
X    }
X    break;
X  case GNOMONIC:
X/* gnomonic projection */
X    rac_l = lon - xfs_ra_cen;
X    Dsinlat = DSIN(lat);
X    Dcoslat = DCOS(lat);
X    Dcosrac_l = DCOS(rac_l);
X    Dsinrac_l = DSIN(rac_l);
X
X    theta = acos(sin_dlcen*Dsinlat + cos_dlcen*Dcoslat*Dcosrac_l);
X
X    if (theta <= 1.57) { /* avoid wrapping */
X      denom = (sin_dlcen*Dsinlat + cos_dlcen*Dcoslat*Dcosrac_l)    / xfs_scale;
X      *yloc = xf_ycen +
X    (int)xfs_inv*
X     (cos_dlcen * Dsinlat - sin_dlcen * Dcoslat * Dcosrac_l) / denom + 0.5;
X      *xloc = xf_xcen - Dcoslat * Dsinrac_l / denom + 0.5;
X      *inregion = ((*xloc >= xf_w_left) && (*xloc <= xf_w_right)
X           && (*yloc <= xf_w_top) && (*yloc >= xf_w_bot));
X    } else *inregion = FALSE;
X    break;
X  case ORTHOGR:
X    rac_l = lon - xfs_ra_cen;
X    Dsinlat = DSIN(lat);
X    Dcoslat = DCOS(lat);
X    Dcosrac_l = DCOS(rac_l);
X    Dsinrac_l = DSIN(rac_l);
X
X    theta = acos(sin_dlcen*Dsinlat + cos_dlcen*Dcoslat*Dcosrac_l);
X
X    if (theta <= 1.57) { /* avoid wrapping */
X      *yloc = xf_ycen +
X    (int)xfs_inv * xfs_scale
X      * (cos_dlcen * Dsinlat - sin_dlcen * Dcoslat * Dcosrac_l);
X      *xloc = xf_xcen - xfs_scale * Dcoslat * Dsinrac_l;
X      *inregion = ((*xloc >= xf_w_left) && (*xloc <= xf_w_right)
X           && (*yloc <= xf_w_top) && (*yloc >= xf_w_bot));
X    } else *inregion = FALSE;
X    break;
X  case RECTANGULAR:
X    if ((xf_west < 0.0) && (lon>(xf_west+360.0))) lon -= 360.0;
X    if ((xf_east > 360.0) && (lon<(xf_east-360.0))) lon += 360.0;
X
X    *yloc = xf_ycen + (int)xfs_inv * xfs_scale * (lat - xfs_dl_cen);
X    *xloc = xf_xcen + xfs_scale * (xfs_ra_cen - lon);
X    *inregion = ((*xloc >= xf_w_left) && (*xloc <= xf_w_right)
X         && (*yloc <= xf_w_top) && (*yloc >= xf_w_bot));
X    break;
X  default:
X    break;
X  }
X}
X
X/* Given x and y of a point on the display,
X   return the latitude and longitude */
Xint invxform(x, y, latp, lonp)
X     int x, y;
X     double *latp, *lonp;
X{
X  int i;
X  int winnum;
X  mapwindow *win;
X
X  /* temporaries to hold values set by initxform */
X  double t_xf_west, t_xf_east, t_xf_north, t_xf_south, t_xf_bottom;
X  int t_xf_w_left, t_xf_w_right, t_xf_w_top, t_xf_w_bot;
X  int t_xf_xcen, t_xf_ycen, t_xf_ybot;
X  double t_xf_c_scale;
X
X  int t_xfs_proj_mode;
X  double t_xfs_ra_cen, t_sin_dlcen, t_cos_dlcen, t_chart_scale;
X  double t_xfs_scale;
X  double t_xfs_inv;
X
X
X  double rho;
X  double R, theta;
X  double l, m, n;
X  double l_, m_, n_;
X
X
X
X  *latp = 0.0;
X  *lonp = 0.0;
X
X  /* First, find which mapwindow the point is in */
X  for (i = 0; i < numwins; i++) {
X    if ((x >= mapwin[i]->x_offset) && (y >= mapwin[i]->y_offset)
X    && (x <= (mapwin[i]->x_offset+mapwin[i]->width))
X    && (y <= (mapwin[i]->y_offset+mapwin[i]->height)))
X      /* point is in window i */
X      break;
X  }
X  if (i == numwins) return -1; /* outside all windows */
X
X  winnum = i;
X  win = mapwin[winnum];
X
X  /* Now, initialize inverse transformation for window winnum */
X  t_xf_west = xf_west;
X  t_xf_east = xf_east;
X  t_xf_north = xf_north;
X  t_xf_south = xf_south;
X  t_xf_bottom = xf_bottom;
X
X  t_xf_xcen = xf_xcen;
X  t_xf_ycen = xf_ycen;
X  t_xf_ybot = xf_ybot;
X
X  t_xf_w_left = xf_w_left;
X  t_xf_w_right = xf_w_right;
X  t_xf_w_top = xf_w_top;
X  t_xf_w_bot = xf_w_bot;
X
X  t_xf_c_scale = xf_c_scale;
X
X  t_xfs_proj_mode = xfs_proj_mode;
X
X  t_xfs_ra_cen = xfs_ra_cen;
X
X  t_sin_dlcen = sin_dlcen;
X  t_cos_dlcen = cos_dlcen;
X  t_chart_scale = chart_scale;
X
X  t_xfs_scale = xfs_scale;
X
X  t_xfs_inv = xfs_inv;
X
X  initxform(win);
X
X  /* Calculate lat and lon */
X  switch (win->proj_mode) {
X  case SANSONS:
X    *latp = (y - xf_ybot)/xfs_scale*xfs_inv + xf_bottom;
X    *lonp = -((x - xf_xcen)/(xfs_scale*DCOS(*latp)) - xfs_ra_cen);
X    break;
X  case GNOMONIC:
X  case ORTHOGR:
X  case STEREOGR:
X    x -= xf_xcen;
X    y -= xf_ycen;
X    y *= xfs_inv;
X
X    R = sqrt((double) (x*x + y*y));
X    theta = atan2((double) y, (double) x);
X
X    /* rho is the angle from the center of the display to the object on the
X       unit sphere. */
X    switch (win-> proj_mode) {
X    case STEREOGR:
X      rho = 2.0*atan(R/(2.0*xfs_scale));
X      break;
X    case GNOMONIC:
X      rho = atan(R/xfs_scale);
X      break;
X    case ORTHOGR:
X      rho = asin(R/xfs_scale);
X      break;
X    }
X
X    /* transform from (rho, theta) to l m n direction cosines */
X    l = sin(rho)*cos(theta);    /* rho and theta are in radians */
X    m = sin(rho)*sin(theta);
X    n = cos(rho);
X
X    /* transform to new declination at center
X       new axes rotated about x axis (l) */
X    l_ = l;
X    m_ = m*sin_dlcen - n*cos_dlcen;
X    n_ = m*cos_dlcen + n*sin_dlcen;
X
X    /* calculate lon and lat */
X    *lonp = atan2(l_, m_) /  0.0174532925199 + xfs_ra_cen - 180.0;
X    *latp = 90 - acos(n_) /  0.0174532925199;
X
X    break;
X  case RECTANGULAR:
X    *latp = (y - xf_ycen) / (xfs_inv * xfs_scale) + xfs_dl_cen;
X    *lonp = (xf_xcen - x) / xfs_scale + xfs_ra_cen;
X    break;
X  default:            /* error */
X    winnum = -1;
X  }
X
X
X  /* restore initxform's variables */
X  xf_west = t_xf_west;
X  xf_east = t_xf_east;
X  xf_north = t_xf_north;
X  xf_south = t_xf_south;
X  xf_bottom = t_xf_bottom;
X
X  xf_xcen = t_xf_xcen;
X  xf_ycen = t_xf_ycen;
X  xf_ybot = t_xf_ybot;
X
X  xf_w_left = t_xf_w_left;
X  xf_w_right = t_xf_w_right;
X  xf_w_top = t_xf_w_top;
X  xf_w_bot = t_xf_w_bot;
X
X  xf_c_scale = t_xf_c_scale;
X
X  xfs_proj_mode = t_xfs_proj_mode;
X
X  xfs_ra_cen = t_xfs_ra_cen;
X
X  sin_dlcen = t_sin_dlcen;
X  cos_dlcen = t_cos_dlcen;
X  chart_scale = t_chart_scale;
X
X  xfs_scale = t_xfs_scale;
X
X  xfs_inv = t_xfs_inv;
X
X  if (*lonp >= 360.0) *lonp -= 360.0;
X  if (*lonp < 0.0) *lonp += 360.0;
X
X  return winnum;
X}
X
X
X
X/* Gnomonic transformation
X   Used to draw vectors as great circles
X   in Gnomonic transform a great circle is projected to a line.
X*/
Xstatic double gt_sin_dlcen, gt_cos_dlcen, gt_chart_scale;
Xstatic double gt_scale;
X
X/* endpoints of west and east boundaries in SANSONS and RECTANGULAR */
Xstatic double gt_wx1, gt_wy1, gt_wx2, gt_wy2;
Xstatic double gt_ex1, gt_ey1, gt_ex2, gt_ey2;
X
X/* midpoint, a and b for north and south boundaries.
X   y = a*x*x + y0 for parabola (b == 0.0)
X   1 = x*x/a + (y-y0)*(y-y0)/b for ellipse */
Xstatic double gt_ny0, gt_na, gt_nb;
Xstatic double gt_sy0, gt_sa, gt_sb;
X
X/* radius for STEREOGRAPHIC, GNOMONIC, and ORTHOGRAPHIC */
Xstatic double gt_r;
X
X/* Can we clip to boundaries analytically? */
Xstatic int gt_use_boundaries;
X
X
Xinit_gt(win)
X     mapwindow *win;
X{
X  double x_1, x_2, y_1, y_2;
X  double r_theta;
X  double adj;
X
X  gt_use_boundaries = TRUE;
X
X  gt_sin_dlcen = DSIN(win->dlcen);
X  gt_cos_dlcen = DCOS(win->dlcen);
X  gt_chart_scale = win->scale * .0174532925199; /* Radians */
X
X/* gt_scale is the conversion factor for size of the picture ( = R) */
X  if (xfs_proj_mode == STEREOGR)
X    gt_scale = MIN(win->height, win->width) / (2.0 * DTAN(win->scale));
X  else
X    gt_scale = MIN(win->height, win->width) / (2.0 * DTAN(win->scale/2.0));
X
X  adj = xf_c_scale*0.9;            /* use boundaries slightly
X                     more restricted than full plot */
X
X  /* calculate boundaries of region */
X  switch (xfs_proj_mode) {
X  case SANSONS:
X  case RECTANGULAR:
X    do_gt(xf_south+adj, xf_west+adj, >_wx1, >_wy1, &r_theta);
X    gt_use_boundaries &= (r_theta <= 1.57);
X    do_gt(xf_north-adj, xf_west+adj, >_wx2, >_wy2, &r_theta);
X    gt_use_boundaries &= (r_theta <= 1.57);
X    do_gt(xf_south+adj, xf_east-adj, >_ex1, >_ey1, &r_theta);
X    gt_use_boundaries &= (r_theta <= 1.57);
X    do_gt(xf_north-adj, xf_east-adj, >_ex2, >_ey2, &r_theta);
X    gt_use_boundaries &= (r_theta <= 1.57);
X
X    do_gt(xf_north-adj, xfs_ra_cen, &x_1, &y_1, &r_theta);
X    gt_use_boundaries &= (r_theta <= 1.57);
X    gt_ny0 = y_1;
X
X    if (fabs(xf_north-adj) > (90 - fabs(win->dlcen))) {
X      /* ellipse */
X      do_gt(xf_north-adj, xfs_ra_cen + 180, &x_2, &y_2, &r_theta);
X      gt_use_boundaries &= (r_theta <= 1.57);
X      gt_nb = (y_2 - y_1)/2;
X      gt_ny0 = y_1 + gt_nb ;
X      gt_nb = gt_nb*gt_nb;
X      do_gt(xf_north-adj, xfs_ra_cen + 90, &x_1, &y_1, &r_theta);
X      gt_use_boundaries &= (r_theta <= 1.57);
X
X      gt_na = x_1*x_1 / (1 - (y_1 - gt_ny0)*(y_1 - gt_ny0)/gt_nb);
X    } else {
X      /* parabola */
X      if (gt_use_boundaries) {
X    gt_nb = 0.0;
X    gt_na = (gt_ey2 - gt_ny0) / (gt_ex2 * gt_ex2);
X    /* error if gt_ex2 == 0.0, as when r_theta was > 1.57 */
X      }
X    }
X
X    do_gt(xf_south+adj, xfs_ra_cen, &x_1, &y_1, &r_theta);
X    gt_use_boundaries &= (r_theta <= 1.57);
X    gt_sy0 = y_1;
X
X    if (fabs(xf_south+adj) > (90 - fabs(win->dlcen))) {
X      /* ellipse */
X      do_gt(xf_south+adj, xfs_ra_cen + 180, &x_2, &y_2, &r_theta);
X      gt_use_boundaries &= (r_theta <= 1.57);
X      gt_sb = (y_2 - y_1)/2;
X      gt_sy0 = y_1 - gt_sb ;
X      gt_sb = gt_sb*gt_sb;
X
X      do_gt(xf_south+adj, xfs_ra_cen + 90, &x_1, &y_1, &r_theta);
X      gt_use_boundaries &= (r_theta <= 1.57);
X      do_gt(xf_south+adj, xfs_ra_cen + 270, &x_2, &y_2, &r_theta);
X      gt_use_boundaries &= (r_theta <= 1.57);
X      gt_sa = (x_2 - x_1)/2;
X      gt_sa = gt_sa*gt_sa;
X    } else {
X      /* parabola */
X      if (gt_use_boundaries) {
X    gt_sb = 0.0;
X    gt_sa = (gt_ey1 - gt_sy0) / (gt_ex1 * gt_ex1);
X    /* error if gt_ex2 == 0.0, as when r_theta was > 1.57 */
X      }
X    }
X    break;
X  case STEREOGR:
X    gt_r = MIN(win->height, win->width)/2.0 - 1;
X    break;
X  case ORTHOGR:
X    gt_use_boundaries = FALSE; /* can't handle this analytically */
X    break;
X  case GNOMONIC:
X    gt_wx1 = gt_wx2 = xf_w_left - xf_xcen + 1;
X    gt_ex1 = gt_ex2 = xf_w_right - xf_xcen - 1;
X    gt_ey1 = gt_wy1 = xf_w_bot - xf_ycen + 1;
X    gt_ey2 = gt_wy2 = xf_w_top - xf_ycen - 1;
X    break;
X  default:            /* error */
X    break;
X  }
X}
X
X
X/* Note, returns xloc and yloc as doubles */
Xdo_gt(lat, lon, xloc, yloc, r_theta)
X     double lat, lon, *r_theta;
X     double *xloc, *yloc;
X{
X  double theta, rac_l;
X  double denom;
X  double Dcoslat, Dsinlat, Dcosrac_l, Dsinrac_l;
X  /* Dcoslat, Dsinlat: of object latitude in degrees = phi
X     Dcosrac_l, Dsinrac_l: of object ra - longditude of center = d(lambda) */
X
X  rac_l = lon - xfs_ra_cen;
X  Dsinlat = DSIN(lat);
X  Dcoslat = DCOS(lat);
X  Dcosrac_l = DCOS(rac_l);
X  Dsinrac_l = DSIN(rac_l);
X
X  *r_theta =
X    theta = acos(gt_sin_dlcen*Dsinlat + gt_cos_dlcen*Dcoslat*Dcosrac_l);
X
X  if (theta <= 1.57) { /* avoid wrapping */
X    denom = (gt_sin_dlcen*Dsinlat + gt_cos_dlcen*Dcoslat*Dcosrac_l) / gt_scale;
X    *yloc = xfs_inv*
X    (gt_cos_dlcen * Dsinlat - gt_sin_dlcen * Dcoslat * Dcosrac_l) / denom;
X    *xloc = - Dcoslat * Dsinrac_l / denom;
X  };
X}
X
X/* Given x and y of a point on the display,
X   return the latitude and longitude */
Xinv_gt(x, y, latp, lonp)
X     double x, y;
X     double *latp, *lonp;
X{
X  double rho;
X  double R, theta;
X  double l, m, n;
X  double l_, m_, n_;
X
X  y *= xfs_inv;
X
X  *latp = 0.0;
X  *lonp = 0.0;
X
X  /* Calculate lat and lon */
X  R = sqrt((double) (x*x + y*y));
X  theta = atan2((double) y, (double) x);
X
X  /* rho is the angle from the center of the display to the object on the
X     unit sphere. */
X  rho = atan(R/gt_scale);
X
X  /* transform from (rho, theta) to l m n direction cosines */
X  l = sin(rho)*cos(theta);    /* rho and theta are in radians */
X  m = sin(rho)*sin(theta);
X  n = cos(rho);
X
X  /* transform to new declination at center
X     new axes rotated about x axis (l) */
X  l_ = l;
X  m_ = m*gt_sin_dlcen - n*gt_cos_dlcen;
X  n_ = m*gt_cos_dlcen + n*gt_sin_dlcen;
X
X  /* calculate lon and lat */
X  *lonp = atan2(l_, m_) / 0.0174532925199 + xfs_ra_cen - 180.0;
X  if (n_ > 1) n_ = 1;
X  if (n_ < -1) n_ = -1;
X  *latp = 90 - acos(n_) / 0.0174532925199;
X
X  if (*lonp >= 360.0) *lonp -= 360.0;
X  if (*lonp < 0.0) *lonp += 360.0;
X}
X
X
X/*
X * clipping extentions (ccount)
X */
X
X/* defines and clipped_at are for use in area drawing,
X   to indicate that a corner has fallen in the area */
X#define NO_CLIP 0
X#define WEST_CLIP 1
X#define EAST_CLIP 2
X#define NORTH_CLIP 4
X#define SOUTH_CLIP 8
X#define RADIUS_CLIP 16
X#define NE_CORNER 6
X#define SE_CORNER 10
X#define NW_CORNER 5
X#define SW_CORNER 9
X
Xstatic int clip_at1, clip_at2;
X
X/* return transformed values clipped so both endpoints are inregion */
X/* return lats and lons */
X/* The line is a great circle on the celestial sphere,
X   or a line in lat and long (e.g. the line 0h 0d to 2h 10d
X                              contains the point at 0.5h 5d).
X
X   Gnomonic transformation maps a great circle to a line.
X
X   There are three possibilities:
X   1) both endpoints are in the region.
X   2) one endpoint is in the region.
X   3) both endpoints are outside of the drawn region.
X
XIn case 1, nothing needs to be done.
XIn case 2, a second point at the intersection of the line and the boundary of
Xthe drawn region is calculated.
XIn case 3, it is possible that some segment of the line is in the drawn region,
Xand two points along the boundary are calculated.
X
XThe boundary of the drawn region, projected through gnomonic transformation,
Xare:
XSansons: lines on east and west, a line, a parabola, or an ellipse
X                to the north and south.
X        ellipse iff declination of boundary > (90 - decl. of center)
XStereographic: circle
XGnomonic: rectangle
XOrthographic: rectangle in orthographic is very comples, use a circle.
XSimple: same as sansons.
X*/
X
X
Xclipr_xform(lat1, lon1, lat2, lon2, xloc1, yloc1, xloc2, yloc2, great_circle,
X      plat1, plon1, plat2, plon2)
X     double  lon1, lat1, lon2, lat2;
X     int *xloc1, *yloc1, *xloc2, *yloc2;
X     int great_circle;     /* draw as great circle */
X     double  *plon1, *plat1, *plon2, *plat2;
X{
X  int Lisin, Risin;
X  double x_1, y_1, x_2, y_2;
X  double theta_1, theta_2;
X  int int_w, int_e, int_n, int_s, int_r1, int_r2;
X  double xw, xe, xn, xs, xr1, xr2, yw, ye, yn, ys, yr1, yr2;
X
X  double Llat, Llon, Rlat, Rlon, Mlat, Mlon;
X  int Lx, Ly, Rx, Ry, Mx, My;
X  int inL, inR, inM;
X
X
X  *plon1 = lon1;
X  *plon2 = lon2;
X  *plat1 = lat1;
X  *plat2 = lat2;
X  clip_at1 = clip_at2 = NO_CLIP;
X  xform(lat1, lon1, xloc1, yloc1, &Lisin);
X  xform(lat2, lon2, xloc2, yloc2, &Risin);
X  if (Lisin && Risin)        /* is already ok: case 1 */
X    return TRUE;
X
X  if (great_circle && gt_use_boundaries) {
X    /* Transform to gnomonic */
X    do_gt(lat1, lon1, &x_1, &y_1, &theta_1);
X    do_gt(lat2, lon2, &x_2, &y_2, &theta_2);
X
X    if ((theta_1 > 1.57) || (theta_2 > 1.57)) /* out of field, skip */
X      return FALSE;
X
X    /* Find intersections with boundaries */
X    switch (xfs_proj_mode) {
X    case SANSONS:
X    case RECTANGULAR:
X      line_intersect(x_1, y_1, x_2, y_2, gt_wx1, gt_wy1, gt_wx2, gt_wy2,
X             &xw, &yw, &int_w);
X      line_intersect(x_1, y_1, x_2, y_2, gt_ex1, gt_ey1, gt_ex2, gt_ey2,
X             &xe, &ye, &int_e);
X
X      if (gt_nb == 0.0) {          /* parabola */
X    para_intersect(x_1, y_1, x_2, y_2, gt_na, gt_ny0,
X               &xn, &yn, &int_n);
X      } else {
X    ellip_intersect(x_1, y_1, x_2, y_2, gt_na, gt_nb, gt_ny0,
X            &xn, &yn, &int_n);
X      }
X
X      if (gt_sb == 0.0) {          /* parabola */
X    para_intersect(x_1, y_1, x_2, y_2, gt_sa, gt_sy0,
X               &xs, &ys, &int_s);
X      } else {
X    ellip_intersect(x_1, y_1, x_2, y_2, gt_sa, gt_sb, gt_sy0,
X            &xs, &ys, &int_s);
X      }
X      int_r1 = int_r2 = FALSE;
X      break;
X    case GNOMONIC:
X      line_intersect(x_1, y_1, x_2, y_2, gt_wx1, gt_wy1, gt_wx2, gt_wy2,
X             &xw, &yw, &int_w);
X      line_intersect(x_1, y_1, x_2, y_2, gt_ex1, gt_ey1, gt_ex2, gt_ey2,
X             &xe, &ye, &int_e);
X      line_intersect(x_1, y_1, x_2, y_2, gt_ex2, gt_ey2, gt_wx2, gt_wy2,
X             &xn, &yn, &int_n);
X      line_intersect(x_1, y_1, x_2, y_2, gt_wx1, gt_wy1, gt_ex1, gt_ey1,
X             &xs, &ys, &int_s);
X      int_r1 = int_r2 = FALSE;
X      break;
X    case STEREOGR:
X      circ_intersect(x_1, y_1, x_2, y_2, gt_r, &xr1, &yr1, &int_r1,
X              &xr2, &yr2, &int_r2);
X      int_w = int_n = int_s = int_e = FALSE;
X    case ORTHOGR:
X      break;
X    default:                  /* error */
X      break;
X    };
X
X
X   if (!(!Lisin && !Risin)) {          /* case 2 */
X      if (int_w) {
X    x_1 = xw; y_1 = yw;
X    if (Risin)
X      clip_at1 = WEST_CLIP;
X    else
X      clip_at2 = WEST_CLIP;
X      } else if (int_e) {
X    x_1 = xe; y_1 = ye;
X    if (Risin)
X      clip_at1 = EAST_CLIP;
X    else
X      clip_at2 = EAST_CLIP;
X      } else if (int_n) {
X    x_1 = xn; y_1 = yn;
X    if (Risin)
X      clip_at1 = NORTH_CLIP;
X    else
X      clip_at2 = NORTH_CLIP;
X      } else if (int_s) {
X    x_1 = xs; y_1 = ys;
X    if (Risin)
X      clip_at1 = SOUTH_CLIP;
X    else
X      clip_at2 = SOUTH_CLIP;
X      } else if (int_r1) {
X    x_1 = xr1; y_1 = yr1;
X    if (Risin)
X      clip_at1 = RADIUS_CLIP;
X    else
X      clip_at2 = RADIUS_CLIP;
X      } else {
X/*    fprintf(stderr, "Error drawing vector\n");
X    fprintf(stderr, "from (%.3f %.3f) to (%.3f %.3f)\n",
X        lat1, lon1, lat2, lon2);*/
X    return FALSE;
X      };
X      if (Lisin) {              /* Need to find new point 2 */
X    inv_gt(x_1, y_1, plat2, plon2);
X    xform(*plat2, *plon2, xloc2, yloc2, &inM);
X      } else {                  /* Need to find new point 1 */
X    inv_gt(x_1, y_1, plat1, plon1);
X    xform(*plat1, *plon1, xloc1, yloc1, &inM);
X      };
X    } else {                  /* case 3 */
X      if (int_w && int_e) {
X    x_1 = xw; y_1 = yw;
X    x_2 = xe; y_2 = ye;
X    clip_at1 = WEST_CLIP;
X    clip_at2 = EAST_CLIP;
X      } else if (int_w && int_n) {
X    x_1 = xw; y_1 = yw;
X    x_2 = xn; y_2 = yn;
X    clip_at1 = WEST_CLIP;
X    clip_at2 = NORTH_CLIP;
X      } else if (int_w && int_s) {
X    x_1 = xw; y_1 = yw;
X    x_2 = xs; y_2 = ys;
X    clip_at1 = WEST_CLIP;
X    clip_at2 = SOUTH_CLIP;
X      } else if (int_e && int_n) {
X    x_1 = xe; y_1 = ye;
X    x_2 = xn; y_2 = yn;
X    clip_at1 = EAST_CLIP;
X    clip_at2 = NORTH_CLIP;
X      } else if (int_e && int_s) {
X    x_1 = xe; y_1 = ye;
X    x_2 = xs; y_2 = ys;
X    clip_at1 = EAST_CLIP;
X    clip_at2 = SOUTH_CLIP;
X      } else if (int_n && int_s) {
X    x_1 = xn; y_1 = yn;
X    x_2 = xs; y_2 = ys;
X    clip_at1 = NORTH_CLIP;
X    clip_at2 = SOUTH_CLIP;
X      } else if (int_r1 && int_r2) {
X    x_1 = xr1; y_1 = yr1;
X    x_2 = xr2; y_2 = yr2;
X    clip_at1 = clip_at2 = RADIUS_CLIP;
X      } else return FALSE;
X
X      inv_gt(x_1, y_1, plat1, plon1);
X      inv_gt(x_2, y_2, plat2, plon2);
X      xform(*plat1, *plon1, xloc1, yloc1, &inM);
X      xform(*plat2, *plon2, xloc2, yloc2, &inM);
X    }
X    return TRUE;
X  } else {                  /* find boundaries by bisection */
X
X    if (!Lisin && !Risin)     /* is hopeless */
X      return FALSE;
X
X    /* Now, one side is in, and the other out.  Make sure we won't have
X       problems with crossing 0h */
X    /* If the difference between lon1 and lon2 is greater than
X       the difference if you subtract 360 from the larger,
X       then shift the larger by 360 degrees */
X
X    if (fabs(MAX(lon1,lon2) - MIN(lon1,lon2))
X    > fabs(MAX(lon1,lon2) - 360.0 - MIN(lon1,lon2)))
X      if (lon2 > 180.0) lon2 -= 360.0;
X      else lon1 -= 360.0;
X
X    Llat = lat1;
X    Llon = lon1;
X    Rlat = lat2;
X    Rlon = lon2;
X    xform(Llat, Llon, &Lx, &Ly, &inL);
X    xform(Rlat, Rlon, &Rx, &Ry, &inR);
X
X
X    /* One endpoint is in.
X       Now use bisection to find point at edge */
X    do {
X      if (great_circle) {
X    gcmidpoint(Llat, Llon, Rlat, Rlon, &Mlat, &Mlon);
X      } else {
X    Mlat = (Llat + Rlat) / 2.0;
X    Mlon = (Llon + Rlon) / 2.0;
X      };
X      xform(Mlat, Mlon, &Mx, &My, &inM);
X
X      if (inL)             /* L in R out */
X    if (inM) {         /* between M and R */
X      Llat = Mlat;
X      Llon = Mlon;
X      inL = inM;
X      Lx = Mx;
X      Ly = My;
X    } else {         /* between M and L */
X      Rlat = Mlat;
X      Rlon = Mlon;
X      inR = inM;
X      Rx = Mx;
X      Ry = My;
X    }
X      else             /* L out R in */
X    if (inM) {         /* between M and L */
X      Rlat = Mlat;
X      Rlon = Mlon;
X      inR = inM;
X      Rx = Mx;
X      Ry = My;
X    } else {         /* between M and R */
X      Llat = Mlat;
X      Llon = Mlon;
X      inL = inM;
X      Lx = Mx;
X      Ly = My;
X    };
X    } while ((fabs((Llat - Rlat)) > xf_c_scale) ||
X         (fabs((Llon - Rlon)) > xf_c_scale));
X
X    if (Lisin) {         /* Left point is in,
X                    bisection found right point */
X      *xloc2 = Lx;         /* Use Lx, Ly, since they're inside bounds */
X      *yloc2 = Ly;
X      *plon2 = Llon;
X      *plat2 = Llat;
X    } else {             /* Bisection found left point */
X      *xloc1 = Rx;         /* Use Rx, Ry, since they're inside bounds */
X      *yloc1 = Ry;
X      *plon1 = Rlon;
X      *plat1 = Rlat;
X    }
X
X    return TRUE;
X  }
X}
X
X/* For the macintosh using MPW C */
X#ifdef macintosh
X#pragma segment 3b
X#endif
X
X#define Fuz 0.1
X
X/* calculate and return the intersection point of two lines given
X   two points on each line */
Xline_intersect(x_1, y_1, x_2, y_2, x_3, y_3, x_4, y_4, x, y, int_1)
X     double x_1, y_1, x_2, y_2, x_3, y_3, x_4, y_4;
X     double *x, *y;
X     int *int_1;            /* FALSE if parallel */
X{
X  double a, b, c, d;
X  int x1, y1;
X  double lat_1, lon_1;
X  int in;
X
X  if (fabs(x_2 - x_1) > 1e-5) {    /* Slope may be calculated */
X    a = (y_2 - y_1)/(x_2 - x_1);
X    b = y_1 - a * x_1;
X    if ((fabs(x_4 - x_3) < 1e-5)) {    /* This slope is infinite */
X      /* calculate intersection */
X      *x = x_3;
X      *y = a*x_3 + b;
X      *int_1 = TRUE;
X    } else {                /* Both slopes may be calculated */
X      c = (y_4 - y_3)/(x_4 - x_3);
X      d = y_3 - c * x_3;
X      if (fabs(a - c) < 1e-5) {        /* Slopes the same, no intersection */
X    *int_1 = FALSE;
X      } else {                /* calculate intersection */
X    *x = (d - b)/(a - c);
X    *y = (a*d - b*c)/(a - c);
X    *int_1 = TRUE;
X      };
X    };
X  } else {                /* Slope is infinite */
X    if ((fabs(x_4 - x_3) < 1e-5)) {    /* this slope is also infinite */
X      *int_1 = FALSE;
X    } else {                /* There's an intersection */
X      c = (y_4 - y_3)/(x_4 - x_3);
X      d = y_3 - c * x_3;
X      *x = x_1;
X      *y = c*x_1 + d;
X      *int_1 = TRUE;
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
Xpara_intersect(x_1, y_1, x_2, y_2, a, b, x, y, int_1)
X     double x_1, y_1, x_2, y_2, *x, *y;
X     double a, b;        /* y = a*x*x + b */
X     int *int_1;
X{
X  double c, d;
END_OF_FILE
if test 31121 -ne `wc -c <'starchart/ssup.c.aa'`; then
    echo shar: \"'starchart/ssup.c.aa'\" unpacked with wrong size!
fi
# end of 'starchart/ssup.c.aa'
fi
echo shar: End of archive 19 \(of 32\).
cp /dev/null ark19isdone
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