earthdistance.c
上传用户:blenddy
上传日期:2007-01-07
资源大小:6495k
文件大小:2k
源码类别:

数据库系统

开发平台:

Unix_Linux

  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <string.h>
  4. #include <postgres.h>
  5. #include <utils/geo_decls.h> /* for Pt */
  6. #include <utils/palloc.h> /* for palloc */
  7. /* Earth's radius is in statute miles. */
  8. const EARTH_RADIUS = 3958.747716;
  9. const TWO_PI = 2.0 * M_PI;
  10. /******************************************************
  11.  *
  12.  * degtorad - convert degrees to radians
  13.  *
  14.  * arg: double, angle in degrees
  15.  *
  16.  * returns: double, same angle in radians
  17.  ******************************************************/
  18. static double
  19. degtorad(double degrees)
  20. {
  21. return (degrees / 360.0) * TWO_PI;
  22. }
  23. /******************************************************
  24.  *
  25.  * geo_distance - distance between points
  26.  *
  27.  * args:
  28.  *  a pair of points - for each point,
  29.  *    x-coordinate is longitude in degrees west of Greenwich
  30.  *    y-coordinate is latitude in degrees above equator
  31.  *
  32.  * returns: double
  33.  *  distance between the points in miles on earth's surface
  34.  ******************************************************/
  35. double *
  36. geo_distance(Point *pt1, Point *pt2)
  37. {
  38. double long1,
  39. lat1,
  40. long2,
  41. lat2;
  42. double longdiff;
  43. double    *resultp = palloc(sizeof(double));
  44. /* convert degrees to radians */
  45. long1 = degtorad(pt1->x);
  46. lat1 = degtorad(pt1->y);
  47. long2 = degtorad(pt2->x);
  48. lat2 = degtorad(pt2->y);
  49. /* compute difference in longitudes - want < 180 degrees */
  50. longdiff = fabs(long1 - long2);
  51. if (longdiff > M_PI)
  52. longdiff = TWO_PI - longdiff;
  53. *resultp = EARTH_RADIUS * acos
  54. (sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(longdiff));
  55. return resultp;
  56. }