开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > 传真(Fax)编程 > 查看源码
geo.cpp
资源名称：GTIFFSample.zip　[点击查看]
上传用户：lbr_007
上传日期：2019-05-31
资源大小：282k
文件大小：18k
源码类别：

传真(Fax)编程

开发平台：

Visual C++

geo.cpp：源码内容
  1. // geo.cpp
  2. //
  4. #include "stdafx.h"
  6. #pragma warning (disable:4786)
  7. #include "geo.h"
  9. #include <math.h>
  10. #include <UtilityParser.h>
  11. #include <StringUtil.h>
  13. double GeoCalc::EllipseArcLength(double lat1, double lat2, double a, double f)
  14. {
  15. double result = 0.0;
  17.   // if this path crosses the equator (i.e. lat1 & lat2 are different signs)
  18.   // then we have to break the path up into 2 parts, from the southern point
  19.   // up to the equator, then from the equator to the northern point.
  20.   // if this path is entirely on one side of the equator we can do it in one
  21.   // part.
  22. UINT32 numParts = 1;
  23. double multLat = lat1 * lat2;
  24. if (multLat < 0.0)
  25. {
  26. numParts = 2;
  27. }
  29. INT32 steps;
  30. double snLat1, snLat2, twoFF, x1, x2, dx, x, y1, y2, dy;
  31. double aOneF, dydx, adx, a2, oneF;
  32. double dlat1, dlat2;
  34. a2 = a * a;
  35. oneF = 1 - f;
  36. aOneF = a * oneF;
  38. for (UINT32 part = 0; part < numParts; part++){
  39. if (numParts == 2)
  40. {
  41. if (part == 0)
  42. {
  43. dlat1 = (lat1 < lat2) ? lat1 : lat2;
  44. dlat2 = 0.0;
  45. }
  46. else
  47. {
  48. dlat1 = 0.0;
  49. dlat2 = (lat1 > lat2) ? lat1 : lat2;
  50. }
  51. }
  52. else
  53. {
  54. dlat1 = (lat1 < lat2) ? lat1 : lat2;
  55. dlat2 = (lat1 > lat2) ? lat1 : lat2;
  56. }
  58. steps = 100;
  59. steps += 100 * (INT32)(0.50 + (dlat2-dlat1));
  60. steps = (steps > 4000) ? 4000 : steps;
  62. snLat1 = sin(GEO::DE2RA*dlat1);
  63. snLat2 = sin(GEO::DE2RA*dlat2);
  64. twoFF = 2 * f - f * f;
  65. x1 = a * cos(GEO::DE2RA * dlat1) / sqrt(1-twoFF*snLat1*snLat1);
  66. x2 = a * cos(GEO::DE2RA * dlat2) / sqrt(1-twoFF*snLat2*snLat2);
  68. dx = (x2 - x1) / (double)(steps - 1);
  69. x, y1, y2, dy, dydx;
  70. adx = (dx < 0.0) ? -dx : dx; // absolute value of dx
  72. for (INT32 i = 0; i < (steps - 1); i++){
  73. x = x1 + dx * i;
  74. dydx = ((aOneF * sqrt((1.0 - ((x+dx)*(x+dx))/a2))) -
  75. (aOneF * sqrt((1.0 - (x*x)/a2)))) / dx;
  76. result += adx * sqrt(1.0 + dydx*dydx);
  77. }
  78. }
  80. return result;
  81. }
  83. double GeoCalc::GCDistance(double lat1, double lon1, double lat2, double lon2)
  84. {
  85. lat1 *= GEO::DE2RA;
  86. lon1 *= GEO::DE2RA;
  87. lat2 *= GEO::DE2RA;
  88. lon2 *= GEO::DE2RA;
  90. double d = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1 - lon2);
  91. return (GEO::AVG_ERAD * acos(d));
  92. }
  94. void GeoCalc::GCDistanceAzimuth(double lat1, double lon1, double dist, double az,
  95.    double& lat2, double& lon2)
  96. {
  97. double b = dist / GEO::AVG_ERAD;
  98. double sinb = sin(b);
  99. double cosb = cos(b);
  100. double sinc = sin(GEO::DE2RA * (90.0 - lat1));
  101. double cosc = cos(GEO::DE2RA * (90.0 - lat1));
  102. double azrad = GEO::DE2RA * az;
  104. double a = acos(cosb*cosc + sinc*sinb*cos(azrad));
  105. double B = asin(sinb*sin(azrad)/sin(a));
  107. lat2 = GEO::RA2DE * ((GEO::PIE/2.0) - a);
  108. lon2 = GEO::RA2DE * B + lon1;
  109. }
  111. double GeoCalc::GCAzimuth(double lat1, double lon1, double lat2, double lon2)
  112. {
  113. double result = 0.0;
  115. INT32 ilat1 = (INT32)(0.50 + lat1 * 360000.0);
  116. INT32 ilat2 = (INT32)(0.50 + lat2 * 360000.0);
  117. INT32 ilon1 = (INT32)(0.50 + lon1 * 360000.0);
  118. INT32 ilon2 = (INT32)(0.50 + lon2 * 360000.0);
  120. lat1 *= GEO::DE2RA;
  121. lon1 *= GEO::DE2RA;
  122. lat2 *= GEO::DE2RA;
  123. lon2 *= GEO::DE2RA;
  125. if ((ilat1 == ilat2) && (ilon1 == ilon2))
  126. {
  127. return result;
  128. }
  129. else if (ilon1 == ilon2)
  130. {
  131. if (ilat1 > ilat2)
  132. result = 180.0;
  133. }
  134. else
  135. {
  136. double c = acos(sin(lat2)*sin(lat1) + cos(lat2)*cos(lat1)*cos((lon2-lon1)));
  137. double A = asin(cos(lat2)*sin((lon2-lon1))/sin(c));
  138. result = (A * GEO::RA2DE);
  140. if ((ilat2 > ilat1) && (ilon2 > ilon1))
  141. {
  142. }
  143. else if ((ilat2 < ilat1) && (ilon2 < ilon1))
  144. {
  145. result = 180.0 - result;
  146. }
  147. else if ((ilat2 < ilat1) && (ilon2 > ilon1))
  148. {
  149. result = 180.0 - result;
  150. }
  151. else if ((ilat2 > ilat1) && (ilon2 < ilon1))
  152. {
  153. result += 360.0;
  154. }
  155. }
  157. return result;
  158. }
  160. double GeoCalc::ApproxDistance(double lat1, double lon1, double lat2, double lon2)
  161. {
  162. lat1 = GEO::DE2RA * lat1;
  163. lon1 = -GEO::DE2RA * lon1;
  164. lat2 = GEO::DE2RA * lat2;
  165. lon2 = -GEO::DE2RA * lon2;
  167. double F = (lat1 + lat2) / 2.0;
  168. double G = (lat1 - lat2) / 2.0;
  169. double L = (lon1 - lon2) / 2.0;
  171. double sing = sin(G);
  172. double cosl = cos(L);
  173. double cosf = cos(F);
  174. double sinl = sin(L);
  175. double sinf = sin(F);
  176. double cosg = cos(G);
  178. double S = sing*sing*cosl*cosl + cosf*cosf*sinl*sinl;
  179. double C = cosg*cosg*cosl*cosl + sinf*sinf*sinl*sinl;
  180. double W = atan2(sqrt(S),sqrt(C));
  181. double R = sqrt((S*C))/W;
  182. double H1 = (3 * R - 1.0) / (2.0 * C);
  183. double H2 = (3 * R + 1.0) / (2.0 * S);
  184. double D = 2 * W * GEO::ERAD;
  185. return (D * (1 + GEO::FLATTENING * H1 * sinf*sinf*cosg*cosg -
  186. GEO::FLATTENING*H2*cosf*cosf*sing*sing));
  187. }
  189. double GeoCalc::EllipsoidDistance(double lat1, double lon1, double lat2, double lon2)
  190. {
  191. double distance = 0.0;
  192. double  faz, baz;
  193. double  r = 1.0 - GEO::FLATTENING;
  194. double  tu1, tu2, cu1, su1, cu2, x, sx, cx, sy, cy, y, sa, c2a, cz, e, c, d;
  195. double  cosy1, cosy2;
  196. distance = 0.0;
  198. if((lon1 == lon2) && (lat1 == lat2)) return distance;
  199. lon1 *= GEO::DE2RA;
  200. lon2 *= GEO::DE2RA;
  201. lat1 *= GEO::DE2RA;
  202. lat2 *= GEO::DE2RA;
  204. cosy1 = cos(lat1);
  205. cosy2 = cos(lat2);
  207. if(cosy1 == 0.0) cosy1 = 0.0000000001;
  208. if(cosy2 == 0.0) cosy2 = 0.0000000001;
  210. tu1 = r * sin(lat1) / cosy1;
  211. tu2 = r * sin(lat2) / cosy2;
  212. cu1 = 1.0 / sqrt(tu1 * tu1 + 1.0);
  213. su1 = cu1 * tu1;
  214. cu2 = 1.0 / sqrt(tu2 * tu2 + 1.0);
  215. x = lon2 - lon1;
  217. distance = cu1 * cu2;
  218. baz = distance * tu2;
  219. faz = baz * tu1;
  221. do {
  222. sx = sin(x);
  223. cx = cos(x);
  224. tu1 = cu2 * sx;
  225. tu2 = baz - su1 * cu2 * cx;
  226. sy = sqrt(tu1 * tu1 + tu2 * tu2);
  227. cy = distance * cx + faz;
  228. y = atan2(sy, cy);
  229. sa = distance * sx / sy;
  230. c2a = -sa * sa + 1.0;
  231. cz = faz + faz;
  232. if(c2a > 0.0) cz = -cz / c2a + cy;
  233. e = cz * cz * 2. - 1.0;
  234. c = ((-3.0 * c2a + 4.0) * GEO::FLATTENING + 4.0) * c2a * GEO::FLATTENING / 16.0;
  235. d = x;
  236. x = ((e * cy * c + cz) * sy * c + y) * sa;
  237. x = (1.0 - c) * x * GEO::FLATTENING + lon2 - lon1;
  238. } while(fabs(d - x) > GEO::EPS);
  240. x = sqrt((1.0 / r / r - 1.0) * c2a + 1.0) + 1.0;
  241. x = (x - 2.0) / x;
  242. c = 1.0 - x;
  243. c = (x * x / 4.0 + 1.0) / c;
  244. d = (0.375 * x * x - 1.0) * x;
  245. x = e * cy;
  246. distance = 1.0 - e - e;
  247. distance = ((((sy * sy * 4.0 - 3.0) *
  248. distance * cz * d / 6.0 - x) * d / 4.0 + cz) * sy * d + y) * c * GEO::ERAD * r;
  250. return distance;
  251. }
  253. // the first path corresponds to { lat1A, lon1A, lat1B, lon1B }
  254. // second path is: { lat2A, lon2A, lat2B, lon2B }
  255. //
  256. // Any 3 non-coincident points form a plane, so the first plane is formed
  257. // by lat1A, lon1A, lat1B, lon1B, and the center of the earth, the second
  258. // plane is formed by lat2A, lon2A, lat2B, lon2B, and the earth's center.
  259. // You can get the normal to each plane by taking the cross product of
  260. // the vector from the Earth's center to the first path point with the 
  261. // vector from the Earth's center to the second path point. The cross
  262. // product guarantees a vector that is normal (perpendicular) to the
  263. // first two, and any vector normal to either of the respective normals
  264. // must lie in the corresponding plane, so by taking the cross product of
  265. // the two normals (normals to each path plane) you get the vector of the
  266. // line of intersection of the planes. From this vector its easy to deduce
  267. // if the intersection point of the planes is on either of the paths.
  268. //
  269. bool GeoCalc::GCIntersectSegment(double lat1A, double lon1A, double lat1B, double lon1B,
  270.    double lat2A, double lon2A, double lat2B, double lon2B,
  271.    double& lat3A, double& lon3A, double& lat3B, double& lon3B)
  272. {
  273. bool isInside = false;
  275. double v1[3], v2[3], v3a[3], v3b[3], n1[3], n2[3];
  276. double m;
  278. double d1 = ApproxDistance(lat1A, lon1A, lat1B, lon1B);
  279. double d2 = ApproxDistance(lat2A, lon2A, lat2B, lon2B);
  281.   //
  282.   // for path 1, setting up my 2 vectors, v1 is vector
  283.   // from center of the Earth to point A, v2 is vector
  284.   // from center of the Earth to point B.
  285.   //
  286. v1[0] = cos(lat1A * GEO::DE2RA) * cos(lon1A * GEO::DE2RA);
  287. v1[1] = sin(lat1A * GEO::DE2RA);
  288. v1[2] = cos(lat1A * GEO::DE2RA) * sin(lon1A * GEO::DE2RA);
  290. v2[0] = cos(lat1B * GEO::DE2RA) * cos(lon1B * GEO::DE2RA);
  291. v2[1] = sin(lat1B * GEO::DE2RA);
  292. v2[2] = cos(lat1B * GEO::DE2RA) * sin(lon1B * GEO::DE2RA);
  294.   //
  295.   // n1 is the normal to the plane formed by the three points:
  296.   // center of the Earth, point 1A, and point 1B
  297.   //
  298. n1[0] = (v1[1]*v2[2]) - (v1[2]*v2[1]);
  299. n1[1] = (v1[2]*v2[0]) - (v1[0]*v2[2]);
  300. n1[2] = (v1[0]*v2[1]) - (v1[1]*v2[0]);
  302.   //
  303.   // for path 2, setting up my 2 vectors, v1 is vector
  304.   // from center of the Earth to point A, v2 is vector
  305.   // from center of the Earth to point B.
  306.   //
  307. v1[0] = cos(lat2A * GEO::DE2RA) * cos(lon2A * GEO::DE2RA);
  308. v1[1] = sin(lat2A * GEO::DE2RA);
  309. v1[2] = cos(lat2A * GEO::DE2RA) * sin(lon2A * GEO::DE2RA);
  311. v2[0] = cos(lat2B * GEO::DE2RA) * cos(lon2B * GEO::DE2RA);
  312. v2[1] = sin(lat2B * GEO::DE2RA);
  313. v2[2] = cos(lat2B * GEO::DE2RA) * sin(lon2B * GEO::DE2RA);
  315.   //
  316.   // n1 is the normal to the plane formed by the three points:
  317.   // center of the Earth, point 2A, and point 2B
  318.   //
  319. n2[0] = (v1[1]*v2[2]) - (v1[2]*v2[1]);
  320. n2[1] = (v1[2]*v2[0]) - (v1[0]*v2[2]);
  321. n2[2] = (v1[0]*v2[1]) - (v1[1]*v2[0]);
  323.   //
  324.   // v3 is perpendicular to both normal 1 and normal 2, so
  325.   // it lies in both planes, so it must be the line of
  326.   // intersection of the planes. The question is: does it
  327.   // go towards the correct intersection point or away from
  328.   // it.
  329.   //
  330. v3a[0] = (n1[1]*n2[2]) - (n1[2]*n2[1]);
  331. v3a[1] = (n1[2]*n2[0]) - (n1[0]*n2[2]);
  332. v3a[2] = (n1[0]*n2[1]) - (n1[1]*n2[0]);
  334.   //
  335.   // want to make v3a a unit vector, so first have to get
  336.   // magnitude, then each component by magnitude
  337.   //
  338. m = sqrt(v3a[0]*v3a[0] + v3a[1]*v3a[1] + v3a[2]*v3a[2]);
  339. v3a[0] /= m;
  340. v3a[1] /= m;
  341. v3a[2] /= m;
  343.   //
  344.   // calculating intersection points 3A & 3B. A & B are
  345.   // arbitrary designations right now, later we make A
  346.   // the one close to, or within, the paths.
  347.   //
  348. lat3A = asin(v3a[1]);
  349. if ((lat3A > GEO::EPS) || (-lat3A > GEO::EPS))
  350. lon3A = asin(v3a[2]/cos(lat3A));
  351. else
  352. lon3A = 0.0;
  354. v3b[0] = (n2[1]*n1[2]) - (n2[2]*n1[1]);
  355. v3b[1] = (n2[2]*n1[0]) - (n2[0]*n1[2]);
  356. v3b[2] = (n2[0]*n1[1]) - (n2[1]*n1[0]);
  358. m = sqrt(v3b[0]*v3b[0] + v3b[1]*v3b[1] + v3b[2]*v3b[2]);
  359. v3b[0] /= m;
  360. v3b[1] /= m;
  361. v3b[2] /= m;
  363. lat3B = asin(v3b[1]);
  364. if ((lat3B > GEO::EPS) || (-lat3B > GEO::EPS))
  365. lon3B = asin(v3b[2]/cos(lat3B));
  366. else
  367. lon3B = 0.0;
  369.   //
  370.   // get the distances from the path endpoints to the two
  371.   // intersection points. these values will be used to determine
  372.   // which intersection point lies on the paths, or which one
  373.   // is closer.
  374.   //
  375. double d1a3a = ApproxDistance(lat1A, lon1A, lat3A, lon3A);
  376. double d1b3a = ApproxDistance(lat1B, lon1B, lat3A, lon3A);
  377. double d2a3a = ApproxDistance(lat2A, lon2A, lat3A, lon3A);
  378. double d2b3a = ApproxDistance(lat2B, lon2B, lat3A, lon3A);
  380. double d1a3b = ApproxDistance(lat1A, lon1A, lat3B, lon3B);
  381. double d1b3b = ApproxDistance(lat1B, lon1B, lat3B, lon3B);
  382. double d2a3b = ApproxDistance(lat2A, lon2A, lat3B, lon3B);
  383. double d2b3b = ApproxDistance(lat2B, lon2B, lat3B, lon3B);
  385. if ((d1a3a < d1) && (d1b3a < d1) && (d2a3a < d2) && (d2b3a < d2))
  386. {
  387. isInside = true;
  388. }
  389. else if ((d1a3b < d1) && (d1b3b < d1) && (d2a3b < d2) && (d2b3b < d2))
  390. {
  391.   //
  392.   // 3b is inside the two paths, so swap 3a & 3b
  393.   //
  394. isInside = true;
  396. m = lat3A;
  397. lat3A = lat3B;
  398. lat3B = m;
  399. m = lon3A;
  400. lon3A = lon3B;
  401. lon3B = m;
  402. }
  403. else
  404. {
  405.   //
  406.   // figure out which one is closer to the path
  407.   //
  408. d1 = d1a3a + d1b3a + d2a3a + d2b3a;
  409. d2 = d1a3b + d1b3b + d2a3b + d2b3b;
  411. if (d1 > d2)
  412. {
  413.   //
  414.   // ok, we are here because 3b {lat3B,lon3B} is closer to the
  415.   // paths, so we need to swap 3a & 3b. the other case is already
  416.   // the way 3a & 3b are organized, so no need to swap
  417.   //
  418. m = lat3A;
  419. lat3A = lat3B;
  420. lat3B = m;
  421. m = lon3A;
  422. lon3A = lon3B;
  423. lon3B = m;
  424. }
  425. }
  427.   //
  428.   // convert the intersection points to degrees
  429.   //
  430. lat3A *= GEO::RA2DE;
  431. lon3A *= GEO::RA2DE;
  432. lat3B *= GEO::RA2DE;
  433. lon3B *= GEO::RA2DE;
  435. return isInside;
  436. }
  438. void GeoCalc::GCSatellite(double lat1, double lon1, double satellite_lon,
  439.    double& az, double& elev)
  440. {
  441. double lat2 = 0.0;
  442. double lon2 = satellite_lon;
  444. // before we convert to radians, lets just get the azimuth
  445. az = GCAzimuth(lat1,lon1,lat2,lon2);
  447. lat1 *= GEO::DE2RA;
  448. lon1 *= GEO::DE2RA;
  449. lon2 *= GEO::DE2RA;
  451. double R = GEO::AVG_ERAD;
  452. double h = GEO::GEOSTATIONARY_ALT;
  454. // calculate "b" just like we did in GCDistance
  455. double b = acos(cos(lon2-lon1)*cos(lat1));//sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1 - lon2);
  457. // calculate "d", the distance to the satellite in kilometers
  458. double d = sqrt(R*R + (R+h)*(R+h) - 2*R*(R+h)*cos(b));
  460. // now calculate the elevation angle above the horizon
  461. elev = GEO::RA2DE * acos((R+h)*sin(b)/d);
  462. }
  464. double GeoCalc::DecodeLatitude(const char * txt)
  465. {
  466. double result = 0.0;
  468. UtilityParser p;
  469. StringUtil su;
  471. std::string arg1, arg2, arg3, arg4;
  473. bool is_negative = false;
  475. std::string tmp = txt;
  476. su.StripSpaces(tmp);
  478. if (su.StartsWith(tmp, "-"))
  479. {
  480. is_negative = true;
  481. }
  483. p.SetDelimiters(" -:,");
  484. p.ParseString(txt);
  486. int size = p.GetSize();
  487. if (size > 0)
  488. {
  489. if (size == 1)
  490. {
  491. if (p.GetArgument(0,arg1))
  492. {
  493. if (p.IsNumber(arg1))
  494. {
  495. result = atof(arg1.c_str());
  496. }
  497. }
  498. }
  499. else if (size == 2)
  500. {
  501. if (p.GetArgument(0,arg1) && p.GetArgument(1,arg2))
  502. {
  503. if (p.IsNumber(arg1))
  504. {
  505. result = atof(arg1.c_str());
  506. if (p.IsNumber(arg2))
  507. {
  508. result += atof(arg2.c_str())/60.0;
  509. }
  510. else
  511. {
  512. su.Uppercase(arg2);
  513. su.StripControlCharacters(arg2);
  514. su.StripSpaces(arg2);
  516. if (arg2 == "N")
  517. {
  518. result = fabs(result);
  519. }
  520. else if (arg2 == "S")
  521. {
  522. result = -1.0000000 * fabs(result);
  523. }
  524. else if (arg2.length() > 0)
  525. {
  526. result = 0;
  527. }
  528. }
  529. }
  530. }
  531. }
  532. else if (size > 2)
  533. {
  534. p.GetArgument(0,arg1);
  535. p.GetArgument(1,arg2);
  536. p.GetArgument(2,arg3);
  537. p.GetArgument(3,arg4);
  539. if (p.IsNumber(arg1) && p.IsNumber(arg2) && p.IsNumber(arg3))
  540. {
  541. result = atof(arg1.c_str()) + atof(arg2.c_str())/60.0 + atof(arg3.c_str())/3600.0;
  542. su.Uppercase(arg4);
  543. su.StripControlCharacters(arg4);
  544. su.StripSpaces(arg4);
  545. if (arg4 == "N")
  546. {
  547. result = fabs(result);
  548. }
  549. else if (arg4 == "S")
  550. {
  551. result = -1.0000000 * fabs(result);
  552. }
  553. }
  554. else if (p.IsNumber(arg1) && p.IsNumber(arg2))
  555. {
  556. result = atof(arg1.c_str()) + atof(arg2.c_str())/60.0;
  557. su.Uppercase(arg3);
  558. su.StripControlCharacters(arg3);
  559. su.StripSpaces(arg3);
  560. if (arg3 == "N")
  561. {
  562. result = fabs(result);
  563. }
  564. else if (arg3 == "S")
  565. {
  566. result = -1.0000000 * fabs(result);
  567. }
  568. }
  569. }
  570. }
  572. if (is_negative)
  573. {
  574. result = -fabs(result);
  575. }
  577. return result;
  578. }
  580. double GeoCalc::DecodeLongitude(const char * txt)
  581. {
  582. double result = 0.0;
  584. UtilityParser p;
  585. StringUtil su;
  587. std::string arg1, arg2, arg3, arg4;
  589. bool is_negative = false;
  591. std::string tmp = txt;
  592. su.StripSpaces(tmp);
  594. if (su.StartsWith(tmp, "-"))
  595. {
  596. is_negative = true;
  597. }
  599. p.SetDelimiters(" -:,");
  600. p.ParseString(txt);
  602. int size = p.GetSize();
  603. if (size > 0)
  604. {
  605. if (size == 1)
  606. {
  607. if (p.GetArgument(0,arg1))
  608. {
  609. if (p.IsNumber(arg1))
  610. {
  611. result = atof(arg1.c_str());
  612. }
  613. }
  614. }
  615. else if (size == 2)
  616. {
  617. if (p.GetArgument(0,arg1) && p.GetArgument(1,arg2))
  618. {
  619. if (p.IsNumber(arg1))
  620. {
  621. result = atof(arg1.c_str());
  622. if (p.IsNumber(arg2))
  623. {
  624. result += atof(arg2.c_str())/60.0;
  625. }
  626. else
  627. {
  628. su.Uppercase(arg2);
  629. su.StripControlCharacters(arg2);
  630. su.StripSpaces(arg2);
  632. if (arg2 == "E")
  633. {
  634. result = fabs(result);
  635. }
  636. else if (arg2 == "W")
  637. {
  638. result = -1.0000000 * fabs(result);
  639. }
  640. else if (arg2.length() > 0)
  641. {
  642. result = 0;
  643. }
  644. }
  645. }
  646. }
  647. }
  648. else if (size > 2)
  649. {
  650. p.GetArgument(0,arg1);
  651. p.GetArgument(1,arg2);
  652. p.GetArgument(2,arg3);
  653. p.GetArgument(3,arg4);
  655. if (p.IsNumber(arg1) && p.IsNumber(arg2) && p.IsNumber(arg3))
  656. {
  657. result = atof(arg1.c_str()) + atof(arg2.c_str())/60.0 + atof(arg3.c_str())/3600.0;
  658. su.Uppercase(arg4);
  659. su.StripControlCharacters(arg4);
  660. su.StripSpaces(arg4);
  661. if (arg4 == "E")
  662. {
  663. result = fabs(result);
  664. }
  665. else if (arg4 == "W")
  666. {
  667. result = -1.0000000 * fabs(result);
  668. }
  669. }
  670. else if (p.IsNumber(arg1) && p.IsNumber(arg2))
  671. {
  672. result = atof(arg1.c_str()) + atof(arg2.c_str())/60.0;
  673. su.Uppercase(arg3);
  674. su.StripControlCharacters(arg3);
  675. su.StripSpaces(arg3);
  676. if (arg3 == "E")
  677. {
  678. result = fabs(result);
  679. }
  680. else if (arg3 == "W")
  681. {
  682. result = -1.0000000 * fabs(result);
  683. }
  684. }
  685. }
  686. }
  688. if (is_negative)
  689. {
  690. result = -fabs(result);
  691. }
  693. return result;
  694. }
  696. void GeoCalc::CoordinateToDMS(double coord, INT32& deg, INT32& minutes, double& secs)
  697. {
  698. double dlat, d, m, s;
  700. dlat = coord + (double)0.000001;
  702. d = (double)(long)dlat;
  703. m = (double)((dlat - d)*60.0);
  704. s = (double)((m - (double)(long)(m))*60.0);
  706. if (s >= 60.00)
  707. {
  708. if (m == 59)
  709. {
  710. d += 1;
  711. if (d < 0.0)
  712. d -= 1;
  713. else
  714. d += 1;
  716. m = 0;
  717. s -= 60.0;
  718. }
  719. else
  720. {
  721. s -= 60.00;
  722. s += 1;
  723. }
  724. }
  726. deg = (int)d;
  727. minutes = (int)m;
  728. secs = s;
  729. }
  731. #define IABS(a) ((a>0) ? a : -a)
  733. void GeoCalc::EncodeLatitude(double lat, std::string& s)
  734. {
  735. char txt[24];
  736. INT32 d, m;
  737. double seconds;
  739. CoordinateToDMS(lat,d,m,seconds);
  740. sprintf(txt,"%d-%.2d-%2.1f %c",IABS(d),IABS(m),fabs(seconds),(lat>0.0)?'N':'S');
  741. s = txt;
  742. }
  744. void GeoCalc::EncodeLongitude(double lon, std::string& s)
  745. {
  746. char txt[24];
  747. INT32 d, m;
  748. double seconds;
  750. CoordinateToDMS(lon,d,m,seconds);
  751. sprintf(txt,"%d-%.2d-%2.1f %c",IABS(d),IABS(m),fabs(seconds),(lon>0.0)?'E':'W');
  752. s = txt;
  753. }