/* * geography.c -- * * This file provides general purpose geography functions. * See the user documentation for more information. * * Copyright (c) 2004 Gordon D. Carrie. All rights reserved. * * Licensed under the Open Software License version 2.1 * * Please address questions and feedback to user0@tkgeomap.org * * @(#) $Id: geography.c,v 1.9 2007/06/26 21:49:30 tkgeomap Exp $ * ******************************************** * */ #include #include #include #include #include #include "geography.h" /* * The following constants enable conversion between floating point degree * and radian measurements and the Angle data type. * * Angles are stored as integer microdegrees. This helps reduce the effects * of round off. A microdegree is about 0.11 meters on the spherical Earth, * so points closer than this are assumed to be at the same location. */ #define UNITPERDEG 1.0e+6 #define DEGPERUNIT 1.0e-6 #define RADPERUNIT 0.017453292519943294892e-6 #define UNITPERRAD 57.29577951308232088e+6 #define MAXDEG (INT_MAX * DEGPERUNIT) #define MINDEG (INT_MIN * DEGPERUNIT) #define MAXRAD (INT_MAX * RADPERUNIT) #define MINRAD (INT_MIN * RADPERUNIT) /* * The following constants give Angle representations of some frequently used * angles. */ #define D360 360000000 #define D270 270000000 #define D180 180000000 #define D90 90000000 static double rearth = 6366707.01949; /* *---------------------------------------------------------------------- * * REarth -- * * This procedure retrieves the radius of the Earth. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ double REarth(void) { return rearth; } /* *---------------------------------------------------------------------- * * SetREarth -- * * This procedure sets the assumed radius of the Earth. * * Results: * None. * * Side effects: * A static variable is set. * *---------------------------------------------------------------------- */ void SetREarth(r) double r; /* New Earth radius */ { rearth = r; } /* *---------------------------------------------------------------------- * * BadAngle -- * * This procedure returns a bogus angle. It is used to indicate * error conditions. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle BadAngle(void) { return -INT_MAX; } /* *---------------------------------------------------------------------- * * AngleIsOK -- * * This procedure determines whether an angle value is valid. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int AngleIsOK(a) Angle a; /* Angle to evaluate */ { return a != BadAngle(); } /* *---------------------------------------------------------------------- * * AngleIsBad -- * * This procedure determines whether an angle value is invalid. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int AngleIsBad(a) Angle a; /* Angle to evaluate */ { return a == BadAngle(); } /* *---------------------------------------------------------------------- * * AngleFmDeg -- * * This procedure converts a degree value to an Angle value. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle AngleFmDeg(deg) double deg; /* Degree value to convert */ { return (deg > MAXDEG || deg < MINDEG) ? BadAngle() : UNITPERDEG * deg + (deg > 0.0 ? 0.5 : -0.5); } /* *---------------------------------------------------------------------- * * AngleToDeg -- * * This procedure converts an Angle value to degrees. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ double AngleToDeg(a) Angle a; /* Angle value to convert */ { return a * DEGPERUNIT; } /* *---------------------------------------------------------------------- * * AngleFmRad -- * * This procedure converts a radian value to an Angle value. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle AngleFmRad(rad) double rad; /* Radian value to convert */ { return (rad > MAXRAD || rad < MINRAD) ? BadAngle() : UNITPERRAD * rad + (rad > 0.0 ? 0.5 : -0.5); } /* *---------------------------------------------------------------------- * * AngleToRad -- * * This procedure converts an Angle value to radians. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ double AngleToRad(a) Angle a; /* Angle value to convert */ { return a * RADPERUNIT; } /* *---------------------------------------------------------------------- * * ISin -- * * This procedure computes the sine of an Angle value. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ double ISin(a) Angle a; /* Angle value */ { return sin(a * RADPERUNIT); } /* *---------------------------------------------------------------------- * * ICos -- * * This procedure computes the cosine of an Angle value. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ double ICos(a) Angle a; /* Angle value */ { return cos(a * RADPERUNIT); } /* *---------------------------------------------------------------------- * * GeoPtFmDeg -- * * This procedure sets a geopoint given latitude and longitude values * in degrees. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt GeoPtFmDeg(dLat, dLon) double dLat; /* Latitude, degrees */ double dLon; /* Longitude, degrees */ { GeoPt geoPt; /* Return value */ geoPt.lat = AngleFmDeg(dLat); geoPt.lon = AngleFmDeg(dLon); return (AngleIsBad(geoPt.lat) || AngleIsBad(geoPt.lon)) ? GeoPtNowhere() : geoPt; } /* *---------------------------------------------------------------------- * * GeoPtFmRad -- * * This procedure sets a geopoint given latitude and longitude values * in radians. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt GeoPtFmRad(dLat, dLon) double dLat; /* Latitude, radians */ double dLon; /* Longitude, radians */ { GeoPt geoPt; geoPt.lat = AngleFmRad(dLat); geoPt.lon = AngleFmRad(dLon); return (AngleIsBad(geoPt.lat) || AngleIsBad(geoPt.lon)) ? GeoPtNowhere() : geoPt; } /* *---------------------------------------------------------------------- * * GeoPtGetDeg -- * * This procedure retrieves latitude and longitude values from a * geopoint in degrees. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ void GeoPtGetDeg(geoPt, dLatPtr, dLonPtr) GeoPt geoPt; /* Geographic point */ double *dLatPtr; /* Recipient of latitude */ double *dLonPtr; /* Recipient of longitude */ { *dLatPtr = AngleToDeg(geoPt.lat); *dLonPtr = AngleToDeg(geoPt.lon); } /* *---------------------------------------------------------------------- * * GeoPtGetRad -- * * This procedure retrieves latitude and longitude values from a * geopoint in degrees. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ void GeoPtGetRad(geoPt, dLatPtr, dLonPtr) GeoPt geoPt; /* Geographic point */ double *dLatPtr; /* Recipient of latitude */ double *dLonPtr; /* Recipient of longitude */ { *dLatPtr = AngleToRad(geoPt.lat); *dLonPtr = AngleToRad(geoPt.lon); } /* *---------------------------------------------------------------------- * * GeoPtIsSomewhere -- * * This procedure determines whether a geopoint is valid. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int GeoPtIsSomewhere(geoPt) GeoPt geoPt; /* Geographic point to evaluate */ { return AngleIsOK(geoPt.lat) && AngleIsOK(geoPt.lon); } /* *---------------------------------------------------------------------- * * GeoPtIsNowhere -- * * This procedure determines whether a geopoint is invalid. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int GeoPtIsNowhere(geoPt) GeoPt geoPt; /* Geographic point to evaluate */ { return AngleIsBad(geoPt.lat) || AngleIsBad(geoPt.lon); } /* *---------------------------------------------------------------------- * * GeoPtNowhere -- * * This procedure returns an invalid geopoint. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt GeoPtNowhere(void) { GeoPt nowhere; /* Return value */ nowhere.lat = nowhere.lon = BadAngle(); return nowhere; } /* *---------------------------------------------------------------------- * * MapPtIsSomewhere -- * * This procedure determines whether a map-point is valid. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int MapPtIsSomewhere(mapPt) MapPt mapPt; /* Map point to evaluate */ { return mapPt.abs != FLT_MAX && mapPt.ord != FLT_MAX; } /* *---------------------------------------------------------------------- * * MapPtIsNowhere -- * * This procedure determines whether a map-point is invalid. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int MapPtIsNowhere(mapPt) MapPt mapPt; /* Map point to evaluate */ { return mapPt.abs == FLT_MAX || mapPt.ord == FLT_MAX; } /* *---------------------------------------------------------------------- * * MapPtNowhere -- * * This procedure returns an invalid map-point. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ MapPt MapPtNowhere(void) { MapPt nowhere = {FLT_MAX, FLT_MAX}; return nowhere; } /* *---------------------------------------------------------------------- * * LatLonToCart -- * * This procedure convert a geopoint to 3D Geocentric Cartesian * coordinates where Earth has unit radius * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ CartPt LatLonToCart(geoPt) GeoPt geoPt; /* Geographic point value to convert */ { double lat, lon; /* Latitude, longitude in radians */ double coslat; CartPt cpt; /* Return value */ GeoPtGetRad(geoPt, &lat, &lon); coslat = cos(lat); cpt.x = coslat * cos(lon); cpt.y = coslat * sin(lon); cpt.z = sin(lat); return cpt; } /* *---------------------------------------------------------------------- * * ScaleMapPt -- * * This procedure scales a map-point. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ MapPt ScaleMapPt(mapPt, scale) MapPt mapPt; /* Map point to scale */ double scale; /* Amount to scale it by */ { if (MapPtIsNowhere(mapPt) || scale <= 0.0) { return MapPtNowhere(); } mapPt.abs *= scale; mapPt.ord *= scale; return mapPt; } /* *---------------------------------------------------------------------- * * GeoDistance -- * * This procedure computes the distance between two geopoints. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle GeoDistance(p1, p2) GeoPt p1; /* First point */ GeoPt p2; /* Second point */ { double lat1, lon1, lat2, lon2; double sin_dlon_2, sin_dlat_2; double a; /* * Reference: * * R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope, * vol. 68, no. 2, 1984, p. 159 * * cited in: * http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1 */ GeoPtGetRad(p1, &lat1, &lon1); GeoPtGetRad(p2, &lat2, &lon2); sin_dlon_2 = sin(0.5 * (lon2 - lon1)); sin_dlat_2 = sin(0.5 * (lat2 - lat1)); a = sqrt(sin_dlat_2 * sin_dlat_2 + cos(lat1) * cos(lat2) * sin_dlon_2 * sin_dlon_2); return AngleFmRad(a > 1.0 ? M_PI : 2.0 * asin(a)); } /* *---------------------------------------------------------------------- * * GeoQuickDistance -- * * This procedure computes a measure of the distance between two geopoints. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ double GeoQuickDistance(p1, p2) GeoPt p1; /* First point */ GeoPt p2; /* Second point */ { double lat1, lon1, lat2, lon2; double sin_dlon_2, sin_dlat_2; /* * Reference: * * R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope, * vol. 68, no. 2, 1984, p. 159 * * cited in: * http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1 */ GeoPtGetRad(p1, &lat1, &lon1); GeoPtGetRad(p2, &lat2, &lon2); sin_dlon_2 = sin(0.5 * (lon2 - lon1)); sin_dlat_2 = sin(0.5 * (lat2 - lat1)); return sin_dlat_2 * sin_dlat_2 + cos(lat1) * cos(lat2) * sin_dlon_2 * sin_dlon_2; } /* *---------------------------------------------------------------------- * * Azimuth -- * * This procedure computes the azimuth (bearing) from one point to another. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle Azimuth(p1, p2) GeoPt p1; /* First point */ GeoPt p2; /* Second point */ { double lat1, lon1, lat2, lon2; double cosDLon, sinDLon, sinDLat, sinMLat; GeoPtGetRad(p1, &lat1, &lon1); GeoPtGetRad(p2, &lat2, &lon2); cosDLon = cos(lon2 - lon1); sinDLon = sin(lon2 - lon1); sinDLat = sin(lat1 - lat2); sinMLat = sin(lat2 + lat1); return AngleFmRad(atan2(cos(lat2) * sinDLon, 0.5 * (sinMLat - sinDLat - (sinMLat + sinDLat) * cosDLon))); } /* *---------------------------------------------------------------------- * * GeoStep -- * * This procedure finds the point a given distance and direction from * a another point. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt GeoStep(geoPt, iDir, iDist) GeoPt geoPt; /* Starting point */ Angle iDir; /* Which way to step */ Angle iDist; /* How far to step */ { double dist, dir; double lat, lon; double cos_dist, sin_dist, cos_dir, sin_dir; double cos_lat, cos_lon, sin_lon, sin_lat; double x, y, z, h_1, h_2, dh; dist = AngleToRad(iDist); cos_dist = cos(dist); sin_dist = sin(dist); dir = AngleToRad(iDir); cos_dir = cos(dir); sin_dir = sin(dir); GeoPtGetRad(geoPt, &lat, &lon); cos_lat = cos(lat); cos_lon = cos(lon); sin_lon = sin(lon); sin_lat = sin(lat); x = cos_dist * cos_lon * cos_lat - sin_dir * sin_dist * sin_lon - cos_lon * cos_dir * sin_dist * sin_lat; y = sin_dir * cos_lon * sin_dist + cos_dist * cos_lat * sin_lon - cos_dir * sin_dist * sin_lon * sin_lat; z = cos_lat * cos_dir * sin_dist + cos_dist * sin_lat; /* * Compute x^2 + y^2 */ h_1 = cos_dist * cos_lat - cos_dir * sin_dist * sin_lat; h_2 = (sin_dir) * (sin_dist); dh = h_1 * h_1 + h_2 * h_2; lat = (dh == 0.0) ? (z > 0 ? M_PI_2 : -M_PI_2) : atan(z / sqrt(dh)); lon = atan2(y, x); return GeoPtFmRad(lat, lon); } /* *---------------------------------------------------------------------- * * GCircleX -- * * This procedure computes the intersection of two great circles. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt GCircleX(ln1pt1, ln1pt2, ln2pt1, ln2pt2) GeoPt ln1pt1; /* First point on first great circle */ GeoPt ln1pt2; /* Second point on first great circle */ GeoPt ln2pt1; /* First point on second great circle */ GeoPt ln2pt2; /* Second point on second great circle */ { CartPt ln1pt1c, ln1pt2c, ln2pt1c, ln2pt2c; /* Points on the great circles in * Cartesian Geocentric coordinates */ CartPt p1; /* Normal to plane containing ln1pt1 ln1pt2 */ CartPt p2; /* Normal to plane containing ln2pt2 ln2pt2 */ CartPt c; /* Vector perpendicular to p1 and p2 */ double cs; /* Normalizing factor */ CartPt m; /* Mean of the four input points */ double dp, dm; /* Distance from mean to +-(cx, cy, cz) */ double lat, lon; /* Lat and lon of return value */ /* * Points in 3D Geocentric Cartesian Coordinates */ ln1pt1c = LatLonToCart(ln1pt1); ln1pt2c = LatLonToCart(ln1pt2); ln2pt1c = LatLonToCart(ln2pt1); ln2pt2c = LatLonToCart(ln2pt2); /* * Normal to plane containing ln1pt1 and ln1pt2 (given by cross product) */ p1.x = ln1pt1c.y * ln1pt2c.z - ln1pt1c.z * ln1pt2c.y; p1.y = ln1pt1c.z * ln1pt2c.x - ln1pt1c.x * ln1pt2c.z; p1.z = ln1pt1c.x * ln1pt2c.y - ln1pt1c.y * ln1pt2c.x; /* * Normal to plane containing ln2pt1 and ln2pt2 (given by cross product) */ p2.x = ln2pt1c.y * ln2pt2c.z - ln2pt1c.z * ln2pt2c.y; p2.y = ln2pt1c.z * ln2pt2c.x - ln2pt1c.x * ln2pt2c.z; p2.z = ln2pt1c.x * ln2pt2c.y - ln2pt1c.y * ln2pt2c.x; /* * Perpendicular to p1 and p2 (another cross product) */ c.x = p1.y * p2.z - p1.z * p2.y; c.y = p1.z * p2.x - p1.x * p2.z; c.z = p1.x * p2.y - p1.y * p2.x; /* * If all pts on one line, undefined intersection */ if (c.x == 0.0 && c.y == 0.0 && c.z == 0.0) { return GeoPtNowhere(); } /* * Normalize */ cs = 1.0 / sqrt(c.x * c.x + c.y * c.y + c.z * c.z); c.x *= cs; c.y *= cs; c.z *= cs; /* * Intersection is at a pole */ if (c.x == 0.0 && c.y == 0.0) { return GeoPtFmDeg((c.z > 0) ? 90.0 : -90.0, 0.0); } /* * Find which of the two intersections of the two great circles is * closer to the mean of the four input points */ m.x = 0.25 * (ln1pt1c.x + ln1pt2c.x + ln2pt2c.x + ln2pt2c.x); m.y = 0.25 * (ln1pt1c.y + ln1pt2c.y + ln2pt2c.y + ln2pt2c.y); m.z = 0.25 * (ln1pt1c.z + ln1pt2c.z + ln2pt2c.z + ln2pt2c.z); dp = (m.x-c.x) * (m.x-c.x) + (m.y-c.y) * (m.y-c.y) + (m.z-c.z) * (m.z-c.z); dm = (m.x+c.x) * (m.x+c.x) + (m.y+c.y) * (m.y+c.y) + (m.z+c.z) * (m.z+c.z); if (dm < dp) { /* * Mean point is closer to (-c.x, -c.y, -c.z) */ c.x = -c.x; c.y = -c.y; c.z = -c.z; } lat = atan(c.z / sqrt(c.x * c.x + c.y * c.y)); lon = atan2(c.y, c.x); return GeoPtFmRad(lat, lon); } /* *---------------------------------------------------------------------- * * DomainLat -- * * This procedure puts an angle into the range -90.0 <= angle <= 90.0 * without changing its sine. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle DomainLat(lat) Angle lat; /* Latitude value */ { if (lat > D360) { lat = lat - D360 * (lat / D360); } else if (lat < 0) { lat = lat + D360 * (-lat / D360 + 1); } if (lat > D90 && lat < D270) { lat = D180 - lat; } else if (lat >= D270) { lat = lat - D360; } return lat; } /* *---------------------------------------------------------------------- * * DomainLon -- * * This procedure puts an angle within 180 degrees of a reference angle. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle DomainLon(lon, refLon) Angle lon; /* Longitude value */ Angle refLon; /* Reference longitude */ { if (lon == refLon) { return lon; } if (lon > refLon + D360) { lon -= D360 * ((lon - refLon) / D360); } else if (lon < refLon - D360) { lon += D360 * ((refLon - lon) / D360); } if (lon > refLon + D180) { return lon - D360; } else if (lon < refLon - D180) { return lon + D360; } else { return lon; } } /* *---------------------------------------------------------------------- * * GwchLon -- * * This procedure is equivalent to DomainLon(lon, 0.0) * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ Angle GwchLon(lon) Angle lon; /* Longitude value */ { if (lon == 0) { return lon; } if (lon > D360) { lon -= D360 * ((lon ) / D360); } else if (lon < -D360) { lon += D360 * ((-lon) / D360); } if (lon > D180) { return lon - D360; } else if (lon < -D180) { return lon + D360; } else { return lon; } } /* *---------------------------------------------------------------------- * * DomainLonPt -- * * This procedure returns a geographic point whose longitude is within * 180 degrees of a reference longitude. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt DomainLonPt(geoPt, refLon) GeoPt geoPt; /* Geographic point */ Angle refLon; /* Reference longitude */ { geoPt.lat = DomainLat(geoPt.lat); geoPt.lon = DomainLon(geoPt.lon, refLon); return geoPt; } /* *---------------------------------------------------------------------- * * GwchLonPt -- * * This procedure is equivalent to DomainLonPt(geoPt, 0.0). * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ GeoPt GwchLonPt(geoPt) GeoPt geoPt; /* Geographic point */ { return DomainLonPt(geoPt, 0); } /* *---------------------------------------------------------------------- * * LonCmp -- * * This procedure compares two longitudes. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ enum LonSgn LonCmp(lon0, lon1) Angle lon0; /* First longitude */ Angle lon1; /* Second longitude */ { lon0 = DomainLon(lon0, lon1); return (lon0 < lon1) ? West : (lon0 > lon1) ? East : PrMd; } /* *---------------------------------------------------------------------- * * LonBtwn -- * * This procedure determines whether a meridian lies between two others. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int LonBtwn(lon, lon0, lon1) Angle lon; /* Longitude to evaluate */ Angle lon0; /* Longitude at one end of interval */ Angle lon1; /* Longitude at other end of interval */ { lon0 = DomainLon(lon0, lon); lon1 = DomainLon(lon1, lon); return ((lon0 > lon1) ? lon0 - lon1 : lon1 - lon0) < D180 && ((lon0 < lon && lon < lon1) || (lon1 < lon && lon < lon0)); } /* *---------------------------------------------------------------------- * * LonBtwn1 -- * * This procedure determines whether a meridian lies between or coincides * with two others, * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int LonBtwn1(lon, lon0, lon1) Angle lon; /* Longitude to evaluate */ Angle lon0; /* Longitude at one end of interval */ Angle lon1; /* Longitude at other end of interval */ { lon0 = DomainLon(lon0, lon); lon1 = DomainLon(lon1, lon); return ((lon0 > lon1) ? lon0 - lon1 : lon1 - lon0) < D180 && ((lon0 < lon && lon <= lon1) || (lon1 < lon && lon <= lon0)); } /* *---------------------------------------------------------------------- * * LatCmp -- * * This procedure compares two latitudes. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ enum LatSgn LatCmp(lat0, lat1) Angle lat0; /* First latitude */ Angle lat1; /* Second latitude */ { return (lat0 < lat1) ? South : (lat0 > lat1) ? North : Eq; } /* *---------------------------------------------------------------------- * * AngleCmp -- * * This procedure compares to angles. * * Results: * See the user documentation. * * Side effects: * None. * *---------------------------------------------------------------------- */ int AngleCmp(a0, a1) Angle a0; /* First angle to compare */ Angle a1; /* Second angle to compare */ { return (a0 < a1) ? -1 : (a0 > a1) ? 1 : 0; }