[comp.graphics] Great Circle algorithms: summary

rdd@walt.cc.utexas.edu (Robert Dorsett) (10/10/89)

A few weeks ago, I wrote:


>I'm looking for references to algorithms that can take two
>coordinates on a sphere, and produce a set of points that
>describe the shortest path between those two points.  The closer
>these algorithms come to carto- graphic techniques of plotting
>great circle tracks in a two-dimensional space, the better.

I got a lot of replies; the differences in the ways people have suggested
the problem be solved are quite interesting.  Three people actually sent code;
the code's at the end.

Thanks to everyone who replied.

This message will be archived on rascal.ics.utexas.edu, where it will be avail-
able by anonymous ftp in misc/av/+navaids



Robert Dorsett                                   
Internet: rdd@rascal.ics.utexas.edu               
UUCP: ...cs.utexas.edu!rascal.ics.utexas.edu!rdd  



--------------



From: ben@lamm.mth.msu.edu (Ben Lotto) 

There is a very easy way of producing the great circle route: Let
(a,b,c) and (d,e,f) be the two points you are trying to join. 
The great circle route joining (a,b,c) and (d,e,f) is the line
segment joining the two points projected back to the sphere.

It is easy to get the 3-dimensional coordinates from the latitude and
longitude; if u and v are the latitude and longitude (I would orient
things so that north and west are positive and east and south are
negative; this gives a right handed system with the north pole at
(0,0,1) and the point of the equator at 0 degrees longitude at (1,0,0)),
then
	a = cos v * cos u
	b = sin v * cos u
	c = sin u

The equation for the line segment:
	x(t) = ta + (1-t)d
	y(t) = tb + (1-t)e
	z(t) = tc + (1-t)f
for t in the interval [0,1].

To get the points on the sphere, divide by the length:
	X(t) = x(t)/(x(t)^2 + y(t)^2 + z(t)^2)^{1/2}
	Y(t) = y(t)/(x(t)^2 + y(t)^2 + z(t)^2)^{1/2}
	Z(t) = z(t)/(x(t)^2 + y(t)^2 + z(t)^2)^{1/2}

To get latitude and longitude:
	Lat(t) = arcsin(Z(t))
	Long(t) = arcsin(Y(t)/cos(Lat(t))

Hope this is what you want.

-B. A. Lotto  (ben@nsf1.mth.msu.edu)
Department of Mathematics/Michigan State University/East Lansing, MI  48824


--------------


>	Long(t) = arcsin(Y(t)/cos(Lat(t))

Robert> program works fine for many cases, but falls apart when the
Robert> output longitude is more than pi/2 or less than -pi/2.  My

I guess the best thing to do is to break things up into cases:

Let A = arcsin(Y(t)/cos(Lat(t)).

Case		Long(t)
====		=======
X>0		A
X<0,Y>0		pi - A
X<0,Y<0		- pi - A	(note that A<0 in this case)

I am still using the convention here that Lat is + for east, - for west.

Let me know if this works better.

-B. A. Lotto  (ben@nsf1.mth.msu.edu)
Department of Mathematics/Michigan State University/East Lansing, MI  48824



--------------



From ted@NMSU.Edu Fri Sep 15 09:29:43 1989


how about an algorithm instead of a reference.

convert the points to unit sphere projections.

form the plane containing the origin and the two points.

intersect the sphere with the plane.


if you want to compute incremental points along the way, here is a
recursive bisection algorithm that computes the a point along the
great circle that is at most delta (an angle) away from the initial
point.

_next_point(a,b:point, delta: angle)
if dot(a,b) <= cos(delta) return b
else return next_point(a,(a+b)/2,delta)

next_point(xa,xb:spherical_coordinates, delta: angle)
return _next_point(convert(xa), convert(xb), delta)

if you want to walk along this bisection path, you can easily modify
this algorithm to output the points generated, or to call a procedure
with the generated points.



--------------



From: thomson@cs.utah.edu (Rich Thomson)

Well, if what you want to do is plot great circle tracks in a cartographic
projection (either sphere based or ellipsoid based), then you really should
check out this book:

        Map Projections -- A Working Manual
	    by John P. Snyder
	U.S. Geological Survey Professional Paper 1395
	(Supersedes USGS Bulletin 1532)
	U.S.G. Printing Office, Washington: 1987
	~385 pages, including projection collage poster

He describes all the nitty-gritty of moving from (lat,long) pairs to (x,y)
pairs, for both the sphere and the ellipsoid, for many commonly used map
projections (including those relevant to sattelite mapping).

I have a piece of postscript code that demonstrates this usage to turn
lattitude, longitude coordinates into (x,y) page positions.  I'm using this
as a starting point for making a PostScript prologue that I can wrap around
aerial photo data so it will plot out the coverage of the individual
frames.  When this is ready, I think I'll be posting it to the net
somewhere (probably comp.lang.postscript or comp.graphics).

						-- Rich
-- 
Rich Thomson	thomson@cs.utah.edu  {bellcore,hplabs,uunet}!utah-cs!thomson
"Tyranny, like hell, is not easily conquered; yet we have this consolation with
us, that the harder the conflict, the more glorious the triumph. What we obtain
too cheap, we esteem too lightly." Thomas Paine, _The Crisis_, Dec. 23rd, 1776



--------------



From: Geoff Peck <geoff@apple.com>

The HP-25 applications book, which came packaged with the calculator,
had several useful navigation algorithms.  If you can't find a copy of the
book (I know, it's an old calculator), I might be able to dig mine up and
type in a summary of the algorithms.



--------------


From: cdl%mpl@ucsd.edu (Carl Lowenstein)

You might try looking where the navigators look --
	_American Practical Navigator_, Bowditch, 
	U.S. Gov't Printing Office, various dates since 1802.

My copy (1962 edition) has on p. 232 section 822:
	Great-circle sailing by computation.

This should handle most of what you want, as soon as you figure out what
a haversine is.

		hav (x) = (1.0 - cos (x)) / 2.0

Some people worry that the earth is not a sphere.


-- 
	carl lowenstein		marine physical lab	u.c. san diego
	{decvax|ucbvax}	!ucsd!mplvax!cdl
	cdl@mpl.ucsd.edu



--------------



From: mdf@ziebmef.mef.org (Matthew Francey)

  Let the "start" point be p1 (a vector) and the end point be p2 (another
vector).
  Let these vectors be unit vectors as well.

  calculate W = (p1 cross p2).
  normalize W, ie W = W/|W|.
  calculate Q = (W cross p1).

  any point on the great circle between p1 and p2 can now be given as:

	p = p1*cos(x) + Q*sin(x),

  where x is on [0, acos(p1 dot p2)]

  A nice clean calculation that only needs a single square root and some
sin and cos evaluations, and gives 3d coordinates directly (which can be
fed into whatever-is-your-current-favourite-map-projection).

  If you want to stay away from  cos() and sin() calls, then another
technique is to just find the midpoint of the great circle between p1 and p2
by normalizing (p1 + p2)/2....  this can be applied recursively to find
mid points between midpointes etc.  There is a point where this technique
will loose out to the above...
-- 
Name: Matthew Francey			 Address: N43o34'13.5" W79o34'33.3" 86m
mdf@ziebmef.mef.org                uunet!{utgpu!moore,attcan!telly}!ziebmef!mdf



--------------



From: uunet.UU.NET!cbmvax!mitchell@cs.utexas.edu (Fred Mitchell - QA)

I have done exactly that in a graphics package that I'm writing. However,
I don't have the code for it handy, so I will attempt to give you a
generalization of what I did. If you need more help, then just E-Mail
me.

Consider your sphere to be a unit sphere, i.e, a radius of 1. You are
given the two positions of interest as unit vectors (it should be trivial
to convert whatever sperical coordinate system you're using to that).
Let's call them vectors U1 and U2 .

Find the unit vector V that is perpendicular to U1 & U2 by:
V = U1 X U2.

You may then trace out a great circle by rotating U1 about V until you
reach U2. The details of this eludes me at present :-( but
if you can get back to me later, I can give you a better description.

Unless someone beats me to it! :-)
	-Mitchell



--------------



From thomson@cs.utah.edu Fri Sep 15 19:22:53 1989


Glad to hear that this is what you were looking for.  Here's a little
example code that projects coordinates for the State of Utah using the
State Plane Coordinate System (a set of projection characterizing data for
each state in the US).  This is taken directly from the book I told you
about.  He goes into great detail about the situations where you'd want to
use the different projections and when its appropriate to use the sphere or
ellipsoid equations.  I know there's a dearth of comments in this code, but
if you have the book it should be easy to see where this came from.  BTW, I
obtained the book from the USGS office in San Francisco for $20.  Given the
professional paper number, you shouldn't have any trouble ordering it from
the USGS.  Good Luck.

							-- Rich

Rich Thomson	thomson@cs.utah.edu  {bellcore,hplabs,uunet}!utah-cs!thomson
"Tyranny, like hell, is not easily conquered; yet we have this consolation with
us, that the harder the conflict, the more glorious the triumph. What we obtain
too cheap, we esteem too lightly." Thomas Paine, _The Crisis_, Dec. 23rd, 1776

----
#include <stdio.h>
#include <math.h>

double dms(degrees, minutes, seconds)
    double degrees, minutes, seconds;
    {
    return degrees + (minutes + seconds/60.0)/60.0;
    }

double deg_2_rad(degrees)
    double degrees;
    {
    return (degrees*M_PI)/180.0;
    }

static double R,			/* radius of sphere */
	phi_0, lambda_0,		/* origin of rectangular coord. */
	phi_1_central, phi_2_central,	/* standard parallels */
	n,				/* cone constant */
	F,
	rho,
	rho_0,				/* radius of lat. circle */
	theta;				/* angle from central merid. */

init_project(scale)
    double scale;
    {
    R = 6370997.0/scale;		/* radius of scaled sphere, meters */

    phi_0 = deg_2_rad(dms(38.0, 20.0, 0.0)); /* origin, 38d20' N */
    lambda_0 = deg_2_rad(dms(111.0, 30.0, 0.0)); /* 111d30' W */

    				/* standard parallels for central Utah Zone */
    phi_1_central = deg_2_rad(dms(39.0, 1.0, 0.0)); /* 39d01' N */
    phi_2_central = deg_2_rad(dms(40.0, 39.0, 0.0)); /* 40d39' N */

    n = log(cos(phi_1_central) / cos(phi_2_central)) /
	log(tan(M_PI_4 + phi_2_central/2.0) /
	    tan(M_PI_4 + phi_1_central/2.0));

    F = cos(phi_1_central) * pow(tan(M_PI_4 + phi_1_central/2.0)/n, n);

    rho_0 = R*F/pow(tan(M_PI_4 + phi_1_central/2.0), n);
    }

project(lambda, phi, x, y)
    double lambda, phi,			/* coords of interest, degrees */
	*x, *y;				/* project coords, units of R */
    {
    double theta = n*(deg_2_rad(lambda) - lambda_0),
	rho = R*F/pow(tan(M_PI_4 + deg_2_rad(phi)/2.0), n);

    *x = rho*sin(theta);
    *y = rho_0 - rho*cos(theta);
    }

double map_scale(argc, argv)
    int argc;
    char *argv[];
    {
    if (argc == 1)
	return 1.0;
    else
	return atof(argv[argc-1]);
    }

main(argc, argv)
    int argc;
    char *argv[];
    {
    double lambda, phi, x, y;

    init_project(map_scale(argc, argv));

    for (lambda = 109.0; lambda <= 114.0; lambda += 0.5)
	{
	for (phi = 37.0; phi <= 42.0; phi += 0.5)
	    {
	    project(lambda, phi, &x, &y);
	    printf("%16.4f W,%16.4f N  %16.4f cm,%16.4f cm\n",
		   lambda, phi, 100.0*x, 100.0*y);
	    }
	}
    }



--------------



From: dmi@peregrine.COM (Dean Inada)

#include <stdio.h>
#include <math.h>
#define pi	3.141593635389793238462643
double lat1,lon1,lat2,lon2;
main(argc,argv)
int argc;
char *argv[];
{
	double cos1,sin1,cos2,sin2,lat,lon,cosd,sind;
	double r,dlat,clon;
	double d,z;
	fscanf(stdin," %lf %lf ",&lat2,&lon2);
	cos2 = cos((90.0-lat2)/180*pi);
	sin2 = sqrt(1-cos2*cos2);
	fprintf(stdout,"%lf %lf\n",lat2,lon2);
	while( !feof(stdin) ){
		lat1 = lat2; lon1 = lon2;  cos1 = cos2; sin1 = sin2;
		fscanf(stdin," %lf %lf ",&lat2,&lon2);
		cos2 = cos((90.0-lat2)/180*pi);
		sin2 = sqrt(1-cos2*cos2);
		cosd = cos((lon1-lon2)/180*pi);
		sind = sin((lon1-lon2)/180*pi);
		cd = cos1*cos2+sin1*sin2*cosd;
		/* d = acos(cd)/pi*180; */ /* distance between the points */
		r = (1+cd)/2;
		z = (cos2+cos1)/2;
		clat = z/sqrt(r);
		dlon = acos((sin1+sin2*cosd)/(2*sqrt(r-z*z)));
		if( sind < 0 ){ dlon = -dlon; };
		fprintf(stdout,"%lf %lf\n",90-(acos(clat)/pi*180),lon1-(dlon/pi*180));
		fprintf(stdout,"%lf %lf\n",lat2,lon2);
	}
}
/* you can save half the input and output conversions by keeping */
/* cos(lat) instead of lat */
/* also, if stdin is piped back from stdout, cd will be constant in the loop */



--------------

From: anonymous :-)



/* file: gclib.c */

#include	<math.h>
#include	"gdadef.h"
#include	"defuncs.h"



#define	pi	3.141593
#define	negpi	-pi
#define	twopi	6.283185
#define pi2	1.570796
#define pi4	0.7853982

#define DTR	0.01745329



#define	SUCCESS	-1
#define FAILURE	0


#define INTVAL	20
#define RANGE	1e-5
#define RNGX	1*DTR		/* used to qualify vertex points */



extern float	lng_sdist();
extern float	lng_width();
extern int	drawing;



	
/***************************************************************

routine gcdist

***************************************************************/
gcdist(slng,slat,elng,elat,vlng,vlat)
double	slng,slat;		/* start long/lat */
double	elng,elat;		/* end long/lat */
double	*vlng,*vlat;		/* vertex long/lat returned as result */

{
int	i;
double	tbl[4][3];		/* to record four sets of vertext points */

double	hav_rho, hav_d, hav_dlng, hav_dlat, dist;
double  dlng, dlat;
double 	coslat,coelat;

double	svdlo,evdlo;		/* temporary term */
double	svlat,evlat;		/* latitude of src/dst vertex */
double	svlng,evlng;		/* longitude of src/dst vertex */
double	sterm,eterm;		/* longitude to apply to src/dst to loc vtx */
double	scrse,ecrse;		/* source and destination courses */

/* compute shortest longitudinal difference between points specified */

dlng = lng_sdist(slng,elng);


/* compute latitudinal difference */

if (sgn(slat) == sgn(elat))
	dlat = abs(slat - elat);
else
	dlat = abs(slat) + abs(elat);
	
/*
 check if latitude difference is zero and either source or
 destination is zero: i.e. equatorial great circle.
*/

if (dlat == 0.0 && slat == 0.0)
	{
	*vlat = 0.0;		/* response caller expects */
	*vlng = 0.0;
	return(SUCCESS);
	}

/*
 compute great circle distance.  Algorithim comes
 from American Practical Navigator - great circle sailings.
*/

hav_dlng = haversin(dlng);
hav_dlat = haversin(dlat);

hav_rho = hav_dlng * cos(slat) * cos(elat);
hav_d = hav_rho + hav_dlat;

dist = acos(1.0 - (2.0 * hav_d));


/* compute latitude and delta-longitude for source point */

coelat = pi2 - elat;
scrse = asin(sin(dlng)*sin(coelat)/sin(dist));  /* source course */
svlat = acos(cos(slat) * sin(scrse));
svdlo = cos(scrse)/sin(svlat);
if ( almost(svdlo,1.0,RANGE) )
	sterm = pi2;
else
	sterm = asin(svdlo);

/* compute latitude and delta-longitude for destination point */

coslat = pi2 - slat;
ecrse = asin(sin(dlng)*sin(coslat)/sin(dist));	/* destination course */
evlat = acos(cos(elat) * sin(ecrse));
evdlo = cos(ecrse) / sin(evlat);
if ( almost(evdlo,1.0,RANGE) )
	eterm = pi2;
else
	eterm = asin(evdlo);

/*
 now the problem is one of actually determining where the vertex points
 are.  They are at:
	(slng + sterm, elng - eterm) or
	(slng - sterm, elng + eterm) or
	(slng + sterm, elng + eterm) or
	(slng - sterm, elng - eterm)
 Since vertex points are 180 apart, we can use this fact to determine
 which is actually the correct set.  Note that set may refer to the same
 vertex point whereapon it is deamed the correct point.
*/

tbl[0][0] = lng_sdist(		/* first set - add sterm; sub eterm */
		tbl[0][1] = slng + sterm,
		tbl[0][2] = elng - eterm);

tbl[1][0] = lng_sdist(		/* second set - sub sterm; add eterm */
		tbl[1][1] = slng - sterm,
		tbl[1][2] = elng + eterm);

tbl[2][0] = lng_sdist(		/* third set - add both terms */
		tbl[2][1] = slng + sterm,
		tbl[2][2] = elng + eterm);

tbl[3][0] = lng_sdist(		/* forth set - subtract both terms */
		tbl[3][1] = slng - sterm,
		tbl[3][2] = elng - eterm);

/*
 loop to locate correct point; determine latitude values for
 each vertex point; return lat/long of vertex point to caller
*/

for(i=0; i<4; i++)		/* check four sets of points */
	{
	dlng = tbl[i][0];		/* fetch longitude difference */
	svlng = tbl[i][1];		/* fetch src longitude vertex */
	evlng = tbl[i][2];		/* fetch dst longitude vertex */
	
	/* normalize to -pi to pi range */

	if (svlng > pi)
		svlng -= twopi;
	if (svlng < -pi)
		svlng += twopi;
	
	/*
	 check if points are 180 degrees opposed (i.e. longitude difference
	 will equal pi) or if points are one and the same (i.e. longitude
	 difference will equal 0).
	*/

	if (almost(dlng,pi,RNGX))
		{		

		/*
		 found correct set and points are distinct 
		 determine which point is closer to the 
		 the vertex point (sterm is distance from source point to 
		 the vertex; eterm is the distance from destination point to 
		 the vertex;
		*/

		if (sterm < eterm)
			{	/* vertex is closer to source point */
			svlat = sgn(slat) * svlat;
			evlat = -svlat;
			}
		else
			{	/* vertex is closer to dest point */
			evlat = sgn(elat) * evlat;
			svlat = -evlat;
			}

		/* return source vertex longitude and latitude */

		*vlat = svlat;
		*vlng = svlng;
		return(SUCCESS);
		}

	else if (almost(dlng,0,RNGX))
		{

		/*
		 vertex points are very close (i.e. difference
		 is almost zero). therefore svlng ~ evlng;
		 determine the sign of the latitude: sign takes
		 the sign of the source point or the dest point
		 depending on which is closer to the vertex.
		*/

		if (sterm < eterm)
			{	/* vertex is closer to src pt */
			svlat *= sgn(slat);
			}
		else
			{	/* vertex is closer to dst pt */
			svlat *= sgn(elat);
			}

		*vlat = svlat;
		*vlng = svlng;
		return(SUCCESS);
		}
	}

return(FAILURE);
}



/********************************************************************

routine gcplot

********************************************************************/
double	gcplot(vlat,lng_interval)
double	vlat;			/* lat of gc vertex */
double  lng_interval;		/* long angle from vertex */

{

/*
 algorithim from American Practical Navigator. compute and return to
 caller, the latitude of a point on the great circle (with specified
 vertex) at the specified angle from the vertex longitude.  Note
 that the latitude found is valid for the interval on either side
 of the vertex.
*/

return(atan(cos(lng_interval)*tan(vlat) ) );
}



/***************************************************************

routine get_disp_pnt

***************************************************************/
get_disp_pnt(tbl,pntcnt,scnwn)
float	tbl[][2];		/* table contain screen values */
int	pntcnt;			/* number of points to inspect */
float	scnwn[];		/* screen coordinates */

{
int	index, sign, i;

float	x;

sign = -1;
index = pntcnt >> 1;			/* start with mid point */
for (i = 1; i < pntcnt; i++)
	{
	x = tbl[index][X];		/* fetch screen value */
	if (scnwn[LEFT] <= x && x <= scnwn[RIGHT])
		return(index);
	else
		{
		index = index + (sign * i);	/* alternate about mid pnt */
		sign = -sign;		/* reverse sign for next iteration */
		}
	}
return(-1);
}




/***************************************************************

routine gcdraw

***************************************************************/
gcdraw(type,mapwn,route,src,dst,lnlen,gplen)
int	type;	/* projection type of maps described by mapwn */

/*
 unprojected map coordinates (long/lats in radians).
*/
float	mapwn[];
int	route;	/* long or short great circle route */

/*
 lat/long coordinates of the source (starting) point.
 lat/long corrdinates of destination point.
*/
float	src[];
float	dst[];

float	lnlen;	/* length of line segment in centimeters */
float	gplen;	/* length of gap segment in centimeters */

{
double 	latvtx, lngvtx;

float	wnd_descr[4];	/* window width and left/bottom boundary */

float	prown[4];	/* projected coordinate space */
float	vwpwn[4];	/* current viewport */
float	scnwn[4];	/* screen coordinate space */

float	pltpt[2];	/* unprojected great circle point to plot */
float	propt[2];	/* projected great circle point to plot */
float	scnpt[2];	/* gc point in scrn coord. ready to plot */

float	tbl[INTVAL+1][2];		/* table of scn coord for plotting */
float	earth_tbl[INTVAL+1][2];		/* table of earth coord */
float	d_lng;
float	delta;

int	i;

/*
 perform transformations on map coordinates to minimize work 
 done inside plotting loop.
*/
				
proj_wnd(type,mapwn,prown);	/* convert to projected values */
get_vwp_size(vwpwn);	/* fetch viewport size */
get_wnd_size(scnwn);	/* fetch current screen size (for later use) */

/*
 correct for aspect ratio distortion.  modify the dimensions of the 
 projected coordinate space (prown) so as to eliminate distortion in
 mapping screen space on to the current viewport (vwpwn).  note that
 prown is merely updated in place.
*/

adj_wnd(type,vwpwn,prown,prown);	/* correct aspect ratio distortion */

/*
 creat descriptor describing window which will be passed to map_pnt.
 this saves some computational overhead.
 Contains width (x direction) and starting X value. Also, height (y
 direction and starting Y value.
*/

wnd_descr[0] = lng_width(prown[LEFT],prown[RIGHT]);
wnd_descr[1] = prown[LEFT];
wnd_descr[2] = prown[TOP] - prown[BOTTOM];
wnd_descr[3] = prown[BOTTOM];

/*
 compute lognitudinal difference between great circle start and end points.
 note this calculation takes into account a counter clockwise
 rotation around the globe (longitude increasing from starting to 
 end points on the great circle).
*/

d_lng = lng_width(src[LNG],dst[LNG]);

/*
 determine if long or short great circle is desired by the caller
 and if appropiate, switch points (since plotting algorithim wants
 to generate points on the great circle by incrementing longitude values).
*/

if ( (d_lng > pi && route == GC_SHORT) || (d_lng < pi && route == GC_LONG) )

	{
	/* exchange points to insure proper plotting direction*/

	swap_pts(src,dst);

	/* recompute difference in longitude for new configuration */

	d_lng = lng_width(src[LNG],dst[LNG]);
	}

/*
 check if great circle follows the equator or is a meridian or
 if vertex points can't be computed. 
 if so, great circle is drawn as a straight line without
 performing any calculations for a great circle of general case.
*/

if	(
  (d_lng == 0.0)
	||
  (gcdist(src[LNG],src[LAT],dst[LNG],dst[LAT],&lngvtx,&latvtx) == FAILURE)
	)

	{
	
	/*
	 equatorial or meridian great circle - or vertex points can't
	 be found therefore put source and destination points
	 into tbl for subsequent conversion below. plot as a straight line.
	*/

	tbl[0][X] = src[LNG];
	tbl[0][Y] = src[LAT];
	tbl[1][X] = dst[LNG];
	tbl[1][Y] = dst[LAT];

	/*
	 loop to convert two points stored in table from 
	 earth coordinates into screen coordinates. return 
	 converted value back to tbl for plotting by plot_tbl;
	*/

	for (i=0; i<2; i++)
		{
		earth_tbl[i][LNG] = pltpt[X] = tbl[i][X];
		earth_tbl[i][LAT] = pltpt[Y] = tbl[i][Y];
		pltpt[Y] = tbl[i][Y];
		proj_pnt(type,pltpt,propt);
		map_pnt(0.0,wnd_descr,propt,scnwn,scnpt);
		tbl[i][X] = scnpt[X];
		tbl[i][Y] = scnpt[Y];
		}

	/*
	 send tbl to be plotted. note that only one segment is to
	 be drawn by plot_tbl.
	*/

	plot_tbl(tbl,1,lnlen,gplen);
	i = get_disp_pnt(tbl,2,scnwn);
	if (i < 0)
		{
		src[LNG] = mapwn[LEFT];
		src[LAT] = mapwn[BOTTOM];
		}
	else
		{
		src[LNG] = earth_tbl[i][LNG];
		src[LAT] = earth_tbl[i][LAT];
		}
	return;
	}	

if ( almost(d_lng,pi,RANGE) )
	{
	printf("gcdraw: meridan gc wraps over pole\n");
	return;
	}

/*
 compute longitude interval to apply at each iteration through
 the plotting loop.
*/

delta = d_lng/INTVAL;		/* delta longitude */

/*
 the following is the plotting loop which will draw the great
 circle.  Remember that pltpt is in unprojected coordinate space;
 propt is in projected coordinate space, and scnpt is in screen
 coordinate space. delta is added to plot point to compute a
 new longitude; gcplot is called to compute the latitude of
 the point on the great circle.
*/

for (i=0; i<= INTVAL; i++)	/* for INTVAL segments: */
	{

	/* determine next longitude value to plot - pltpt */

	if (i == 0)
		{	/* use initial values */
		pltpt[X] = src[LNG];
		pltpt[Y] = src[LAT];
		}
	else
		{	/* compute lng and lat value */
		pltpt[X] = (i==INTVAL)? dst[LNG] : pltpt[X] + delta;
		d_lng = lng_sdist(pltpt[X],lngvtx);
		pltpt[Y] = gcplot(latvtx,d_lng);
		}

	normalize_pnt(pltpt);	/* wrap longitude if necessary */

	/*
	 translate great circle point (pltpt) into project coordinate
	 space.
	*/

	proj_pnt(type,pltpt,propt);

	/*
	 map projected point from projected window space into screen
	 space.
	*/

	if (propt[X] > prown[RIGHT] && propt[X] <= prown[RIGHT] + delta)
		map_pnt(delta,wnd_descr,propt,scnwn,scnpt);
	else
		map_pnt(0.0,wnd_descr,propt,scnwn,scnpt);

	/*
	 save data value in table for subsequent plotting by 
	 plot_func.
	*/

	tbl[i][X] = scnpt[X];
	tbl[i][Y] = scnpt[Y];

	/*
	 save true longitude and latitude of point on the
	 great circle.  returne to caller as a possible connector
	 point.
	*/

	earth_tbl[i][LNG] = pltpt[X];
	earth_tbl[i][LAT] = pltpt[Y];
	}


/*
 plot the great circle by calling plot_tbl. table values were generated
 above. Note that line length and gap length are determined by 
 the caller of gcdraw.
*/

plot_tbl(tbl,INTVAL,lnlen,gplen);

/*
 locate a displayable point on the great circle near midway
 in the range of longitudes.  note that number of points (INTVAL+1)
 is passed to function.
*/

i = get_disp_pnt(tbl,INTVAL+1,scnwn);

/*
 check if a display able point was found in the table of points
 descriping the great circle.  if so, return the longitude/latitude
 values describing the point in src point descriptor; otherwise
 return the coordinates of the lower left corner of the map window.
*/

if (i < 0)
	{	/* no displayable point exists on the great circle */
	src[LNG] = mapwn[LEFT];
	src[LAT] = mapwn[BOTTOM];
	}
else
	{	/* displayable point exits with index i in table */
	src[LNG] = earth_tbl[i][LNG];
	src[LAT] = earth_tbl[i][LAT];
	}
return;
}


/* end of gclib.c */