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 */