Category : Miscellaneous Language Source Code
Archive   : QTAWKU42.ZIP
Filename : GEODH.EXP

# QTAwk utility to compare iterative and exact solutions for
# (C) Copyright 1990 Terry D. Boldt, All Rights Reserved
#
# Geodetic Latitude from x, y, z versus iterative solution of DSM
# Method:
# 1) Input Geodetic Latitude, Longitude
# 2) Compute x, y, z (earth centered)
# 3) from x, y, z compute Geodetic Latitude, Longitude:
# a) from iterative solution of DSM
# b) from exact solution
# 4) compare iterative solution for Geodetic Latitude with input Latitude
# 5) compare exact solution for Geodetic Latitude with input Latitude
# 6) Compute two new x, y, z positions:
# a) from iterative solution for Geodetic Latitude
# a) from exact solution for Geodetic Latitude
# 7) compare positions from 6) with positions from 2)
#
# Iterative solution for Geodetic Latitude: (from DSM)
# Input:
# 1) x, y, z position co-ordinates
# 2) Re, equatorial radius of earth == 6,378,137 m (from 84 WGS)
# 3) f, earth flattening == 1/298.257223563 (from 84 WGS)
# 4) conv, convergence criteria == 1e-300 (strigent guess)
#
# Algorithm:
# 1) compute ecentricity of earth ellipsoid == f*(2-f) == ecen ^ 2
# 2) initial guess on zi = -ecen2 * z
# 3) compute following until absolute value of iteration difference for
# zi less than convergence criteria:
# zib = z - zi
# N = sqrt(x^2 + y^2 + zib^2)
# sp = zib/N
# N = Re/sqrt(1 - ecen2 * sp^2)
# zi = -N * ecen2 * sp
# 4) compute iterative Geodetic Latitude:
# latz = arcsine(zib/N)
# 5) compute Longitude:
# longz = atan2(y,x)
#
# Exact solution for Geodetic Latitude:
# Inputs:
# 1) x, y, z position co-ordinates
# 2) Re, equatorial radius of earth == 6378137 m (from 84 WGS)
# 3) f, earth flattening == 1/298.257223563 (from 84 WGS)
#
# Algorithm:
# 1) compute rho = sqrt(x^2 + y^2)
# 2) compute Geodetic Latitude:
# late = atan2(z,rho * (1 - ecen2) )
# 3) compute Longitude:
# longz = atan2(y,x)
#

BEGIN {
OFMT = "%.14g";
Re = 6378137;
flat = 1/298.257223563;
ecen2 = flat*(2 - flat);
ecen = sqrt(ecen2);
conv = 1e-300;
DEGREES = TRUE;
}

{
lat = $1;
long = $2;
print "Input Lat = "lat"ø, Long = "long"ø";
N = Re/sqrt(1 - ecen2*sin(lat)^2);
rho = N * cos(lat);
x = rho * cos(long);
y = rho * sin(long);
z = N * (1 - ecen2) * sin(lat);
print "x = "addcomma(x)" m, y = "addcomma(y)" m, z = "addcomma(z)" m";
rho2 = x^2 + y^2;
for ( zi = -ecen2 * z , i = 0 ; conv < abs(zi - zl) || i >= 1000 ; i++ ) {
zl = zi;
zib = z - zi;
N = sqrt(rho2 + zib^2);
sp = zib/N;
N = Re/sqrt(1 - ecen2 * sp^2);
zi = -N * ecen2 * sp;
}
if ( i >= 1000 ) print "No Convergence";
else print "No iterations: "i"\nëzi = "abs(zi-zl);
latz = asin(zib/N);
longz = atan2(y,x);
print "Iterative Derived Lat/Long\nLat = "latz"ø, Long = "longz"ø";
print "Delta ëLat = "(latz - lat);
rho = sqrt(rho2);
late = atan2(z,rho * (1 - ecen2));
print "Exact Derived Lat/Long\nLat = "late"ø, Long = "longz"ø";
print "Delta ëLat = "(late - lat)"ø";

h = sqrt(z^2 + rho2 * (1 - ecen2)^2)
- (Re * (1 - ecen2))/sqrt(1 - ecen2 * sin(lat)^2);

print "Height above ellipsoid (h): "h;

N = Re/sqrt(1 - ecen2*sin(latz)^2);
rho = N * cos(latz);
xz = rho * cos(longz);
yz = rho * sin(longz);
zz = N * (1 - ecen2) * sin(latz);
RVz = sqrt((xz-x)^2 + (yz-y)^2 + (zz-z)^2);
print "Position from Iterative\nx = "addcomma(xz)" m, y = "addcomma(yz)" m, z = "addcomma(zz)" m";
print "Delta: ëx = "(xz-x)" m, ëy = "(yz-y)" m, ëz = "(zz-z)" m\n|R| = "RVz" m";
N = Re/sqrt(1 - ecen2*sin(late)^2);
rho = N * cos(late);
xe = rho * cos(longz);
ye = rho * sin(longz);
ze = N * (1 - ecen2) * sin(late);
RVe = sqrt((xe-x)^2 + (ye-y)^2 + (ze-z)^2);
RVd = sqrt((xe-xz)^2 + (ye-yz)^2 + (ze-zz)^2);
print "Position from Exact\nx = "addcomma(xe)" m, y = "addcomma(ye)" m, z = "addcomma(ze)" m";
print "Delta: ëx = "(xe-x)" m, ëy = "(ye-y)" m, ëz = "(ze-z)" m\n|R| = "RVe" m";
print "Position Delta\nëx = "(xe-xz)" m, ëy = "(ye-yz)" m, ëz = "(ze-zz)" m\n|ëR| = "RVd" m\f";
}

function abs(x) {
return x > 0 ? x : -x;
}

# function to add commas to numbers
function addcomma(x) {
local num;
local spat;
local bnum = /{_d}{3,3}([,.]|$)/;

if ( x < 0 ) return "-" addcomma(-x);
num = sprintf("%.14g",x); # num is dddddd.dd
# if ( fract(x) ) spat = /{_d}{4,4}[.,]/; else spat = /{_d}{4,4}(,|$)/;
# spat = fract(x) ? /{_d}{4,4}[.,]/ : /{_d}{4,4}(,|$)/;
spat = num ~~ /\./ ? /{_d}{4,4}[.,]/ : /{_d}{4,4}(,|$)/;
while ( num ~~ spat ) sub(bnum,",&",num);
return num;
}