home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / qtawk / geodh.exp < prev    next >
Text File  |  1990-04-23  |  5KB  |  137 lines
  1. # QTAwk utility to compare iterative and exact solutions for
  2. # (C) Copyright 1990 Terry D. Boldt, All Rights Reserved
  3. #
  4. # Geodetic Latitude from x, y, z versus iterative solution of DSM
  5. # Method:
  6. # 1) Input Geodetic Latitude, Longitude
  7. # 2) Compute x, y, z (earth centered)
  8. # 3) from x, y, z compute Geodetic Latitude, Longitude:
  9. #    a) from iterative solution of DSM
  10. #    b) from exact solution
  11. # 4) compare iterative solution for Geodetic Latitude with input Latitude
  12. # 5) compare exact     solution for Geodetic Latitude with input Latitude
  13. # 6) Compute two new x, y, z positions:
  14. #    a) from iterative solution for Geodetic Latitude
  15. #    a) from exact     solution for Geodetic Latitude
  16. # 7) compare positions from 6) with positions from 2)
  17. #
  18. # Iterative solution for Geodetic Latitude: (from DSM)
  19. #   Input:
  20. #     1) x, y, z position co-ordinates
  21. #     2) Re, equatorial radius of earth == 6,378,137 m (from 84 WGS)
  22. #     3) f, earth flattening == 1/298.257223563 (from 84 WGS)
  23. #     4) conv, convergence criteria == 1e-300 (strigent guess)
  24. #
  25. #   Algorithm:
  26. #   1) compute ecentricity of earth ellipsoid == f*(2-f) == ecen ^ 2
  27. #   2) initial guess on zi = -ecen2 * z
  28. #   3) compute following until absolute value of iteration difference for
  29. #      zi less than convergence criteria:
  30. #    zib = z - zi
  31. #    N = sqrt(x^2 + y^2 + zib^2)
  32. #    sp = zib/N
  33. #    N = Re/sqrt(1 - ecen2 * sp^2)
  34. #    zi = -N * ecen2 * sp
  35. #   4) compute iterative Geodetic Latitude:
  36. #    latz = arcsine(zib/N)
  37. #   5) compute Longitude:
  38. #    longz = atan2(y,x)
  39. #
  40. # Exact solution for Geodetic Latitude:
  41. #   Inputs:
  42. #     1) x, y, z position co-ordinates
  43. #     2) Re, equatorial radius of earth == 6378137 m (from 84 WGS)
  44. #     3) f, earth flattening == 1/298.257223563 (from 84 WGS)
  45. #
  46. #   Algorithm:
  47. #   1) compute rho = sqrt(x^2 + y^2)
  48. #   2) compute Geodetic Latitude:
  49. #    late = atan2(z,rho * (1 - ecen2) )
  50. #   3) compute Longitude:
  51. #    longz = atan2(y,x)
  52. #
  53.  
  54. BEGIN {
  55.     OFMT = "%.14g";
  56.     Re = 6378137;
  57.     flat = 1/298.257223563;
  58.     ecen2 = flat*(2 - flat);
  59.     ecen = sqrt(ecen2);
  60.     conv = 1e-300;
  61.     DEGREES = TRUE;
  62. }
  63.  
  64.     {
  65.     lat = $1;
  66.     long = $2;
  67.     print "Input Lat = "lat"°, Long = "long"°";
  68.     N = Re/sqrt(1 - ecen2*sin(lat)^2);
  69.     rho = N * cos(lat);
  70.     x = rho * cos(long);
  71.     y = rho * sin(long);
  72.     z = N * (1 - ecen2) * sin(lat);
  73.     print "x = "addcomma(x)" m, y = "addcomma(y)" m, z = "addcomma(z)" m";
  74.     rho2 = x^2 + y^2;
  75.     for ( zi = -ecen2 * z , i = 0 ; conv < abs(zi - zl) || i >= 1000 ; i++ ) {
  76.     zl = zi;
  77.     zib = z - zi;
  78.     N = sqrt(rho2 + zib^2);
  79.     sp = zib/N;
  80.     N = Re/sqrt(1 - ecen2 * sp^2);
  81.     zi = -N * ecen2 * sp;
  82.     }
  83.     if ( i >= 1000 ) print "No Convergence";
  84.       else print "No iterations: "i"\nδzi = "abs(zi-zl);
  85.     latz = asin(zib/N);
  86.     longz = atan2(y,x);
  87.     print "Iterative Derived Lat/Long\nLat = "latz"°, Long = "longz"°";
  88.     print "Delta δLat = "(latz - lat);
  89.     rho = sqrt(rho2);
  90.     late = atan2(z,rho * (1 - ecen2));
  91.     print "Exact Derived Lat/Long\nLat = "late"°, Long = "longz"°";
  92.     print "Delta δLat = "(late - lat)"°";
  93.  
  94.     h = sqrt(z^2 + rho2 * (1 - ecen2)^2)
  95.      - (Re * (1 - ecen2))/sqrt(1 - ecen2 * sin(lat)^2);
  96.  
  97.     print "Height above ellipsoid (h): "h;
  98.  
  99.     N = Re/sqrt(1 - ecen2*sin(latz)^2);
  100.     rho = N * cos(latz);
  101.     xz = rho * cos(longz);
  102.     yz = rho * sin(longz);
  103.     zz = N * (1 - ecen2) * sin(latz);
  104.     RVz = sqrt((xz-x)^2 + (yz-y)^2 + (zz-z)^2);
  105.     print "Position from Iterative\nx = "addcomma(xz)" m, y = "addcomma(yz)" m, z = "addcomma(zz)" m";
  106.     print "Delta: δx = "(xz-x)" m, δy = "(yz-y)" m, δz = "(zz-z)" m\n|R| = "RVz" m";
  107.     N = Re/sqrt(1 - ecen2*sin(late)^2);
  108.     rho = N * cos(late);
  109.     xe = rho * cos(longz);
  110.     ye = rho * sin(longz);
  111.     ze = N * (1 - ecen2) * sin(late);
  112.     RVe = sqrt((xe-x)^2 + (ye-y)^2 + (ze-z)^2);
  113.     RVd = sqrt((xe-xz)^2 + (ye-yz)^2 + (ze-zz)^2);
  114.     print "Position from Exact\nx = "addcomma(xe)" m, y = "addcomma(ye)" m, z = "addcomma(ze)" m";
  115.     print "Delta: δx = "(xe-x)" m, δy = "(ye-y)" m, δz = "(ze-z)" m\n|R| = "RVe" m";
  116.     print "Position Delta\nδx = "(xe-xz)" m, δy = "(ye-yz)" m, δz = "(ze-zz)" m\n|δR| = "RVd" m\f";
  117. }
  118.  
  119. function abs(x) {
  120.     return x > 0 ? x : -x;
  121. }
  122.  
  123. # function to add commas to numbers
  124. function addcomma(x) {
  125.     local num;
  126.     local spat;
  127.     local bnum = /{_d}{3,3}([,.]|$)/;
  128.  
  129.     if ( x < 0 ) return "-" addcomma(-x);
  130.     num = sprintf("%.14g",x);        # num is dddddd.dd
  131. #    if ( fract(x) ) spat = /{_d}{4,4}[.,]/; else spat = /{_d}{4,4}(,|$)/;
  132. #    spat = fract(x) ? /{_d}{4,4}[.,]/ : /{_d}{4,4}(,|$)/;
  133.     spat = num ~~ /\./ ? /{_d}{4,4}[.,]/ : /{_d}{4,4}(,|$)/;
  134.     while ( num ~~ spat ) sub(bnum,",&",num);
  135.     return num;
  136. }
  137.