<?php # -*- PHP -*-
#
# Thy (english) reference is the original Fortran code!
#
#
# Tämä aliohjelma laskee ns. "geodeettisen
# käänteisfunktion ongelman" tavalla, jonka
# ensimmäinen esittäjä on T.Vincenty.
# Tämä on muunnelma Rainsford:n menetelmästä
# käyttäen Helmert:n elliptisiä termejä.
#
# Tämä toimii kaikkiin suuntiin ja etäisyyksiin
# lukuunottamatta ns. antipodal tilannetta, eli
# että pisteet ovat täsmälleen vastakkaisilla
# puolilla. Lähde-/kohdepiste ei saa myöskään
# olla navalla.
#
#
# Input: $WGS1[La] Source Latitude, WGS84, radians
# $WGS1[Lo] Source Longitude, WGS84, radians
# $WGS2[La] Destination Latitude, WGS84, radians
# $WGS2[Lo] Destination Longitude, WGS84, radians
# $EL[a] Reference ellipsoid semimajor axis
# $EL[f] Reference ellipsoid flattening
#
# Example ellipsoid:
# $EL[a] = 6378137.0
# $EL[f] = 1.0/298.25722210088
# a.k.a. 'GRS80'
# $EL[a] = 6378137.0
# $EL[f] = 1.0/298.257223563
# a.k.a. 'WGS84'
#
# Output: $R[Az] Az from source to dest, north 0, clockwise, 0..2*Pi
# $R[rAz] Az from dest to source, ... (in radians)
# $R[dist] Distance in ellipsoidal units (Depends on what
# units you use for $EL[a]: meters/miles... )
#
function ellipsoidal_great_circle($WGS1, $WGS2, $EL) {
global $Debug;
$iter = 40;
$epsilon = 0.5E-13;
$r = 1.0 - $EL[f];
if (abs($WGS1[La] - $WGS2[La]) < $epsilon &&
abs($WGS1[Lo] - $WGS2[Lo]) < $epsilon) {
$R[Az] = 0.0;
$R[rAz] = 0.0;
$R[dist] = 0.0;
return $R;
}
#
# Formulate initial estimates
#
$tu1 = $r * sin($WGS1[La])/cos($WGS1[La]);
$tu2 = $r * sin($WGS2[La])/cos($WGS2[La]);
$cu1 = 1.0 / sqrt($tu1 * $tu1 + 1.0);
$su1 = $cu1 * $tu1;
$cu2 = 1.0 / sqrt($tu2 * $tu2 + 1.0);
$s = $cu1 * $cu2;
$baz = $s * $tu2;
$faz = $baz * $tu1;
$x = $WGS2[Lo] - $WGS1[Lo];
#
# Now iterate until the error is below 'epsilon'.
#
while (--$iter > 0) { # Max limit, usually 3-4 rounds..
$sx = sin($x);
$cx = cos($x);
$tu1 = $cu2 * $sx;
$tu2 = $baz - $su1 * $cu2 * $cx;
$sy = sqrt($tu1 * $tu1 + $tu2 * $tu2);
$cy = $s * $cx + $faz;
$y = atan2($sy,$cy);
$sa = $s * $sx / $sy;
$c2a = 1.0 - $sa * $sa;
$cz = $faz + $faz;
if ($c2a > 0.0)
$cz = $cy - $cz/$c2a;
$e = $cz * $cz * 2.0 - 1.0;
$c = ((4.0 - 3.0 * $c2a) * $EL[f] + 4.0)
* $c2a * $EL[f] / 16.0;
$d = $x; # Previous version of X
$x = (($e * $cy * $c + $cz) * $sy * $c + $y) * $sa;
$x = (1.0 - $c) * $x * $EL[f] + $WGS2[Lo] - $WGS1[Lo];
if ($Debug)
printf("ellip-GCD: x-x'/epsilon = %f\n",($d-$x)/$epsilon);
if (abs($d-$x) < $epsilon) break;
}
$faz = atan2($tu1,$tu2);
$baz = atan2($cu1 * $sx, $baz * $cx - $su1 * $cu2) + M_PI;
$x = sqrt(((1.0 / ($r * $r)) - 1.0) * $c2a + 1.0) + 1.0;
$x = ($x - 2.0) / $x;
$c = 1.0 - $x;
$xx = $x * $x;
$c = ($xx / 4.0 + 1.0) / $c;
$d = (0.375 * $xx - 1.0) * $x;
$x = $e * $cy;
$s = 1.0 - $e - $e;
$s = (((($sy * $sy * 4.0 - 3.0)
* $s * $cz * $d / 6.0 - $x)
* $d / 4.0 + $cz)
* $sy * $d + $y)
* $c * $EL[a] * $r;
# Normalize into 0..2Pi range
if ($faz < 0) $faz += M_PI + M_PI;
if ($baz < 0) $baz += M_PI + M_PI;
$R[dist] = $s;
$R[Az] = $faz;
$R[rAz] = $baz;
return $R;
}
#
# SUBROUTINE INVER1(GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S)
# C
# C *** SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY
# C *** MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
# C *** EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
# C *** STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
# C
# C *** A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
# C *** F IS THE FLATTENING (NOT RECIPROCAL) OF THE REFERNECE ELLIPSOID
# C *** LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST
# C *** FORWARD AZIMUTHS AT BOTH POINTS RETURNED IN RADIANS FROM NORTH
# C
# C *** PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 18FEB75
# C *** MODIFIED FOR IBM SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 7507
# C
# IMPLICIT REAL*8 (A-H,O-Z)
# COMMON/CONST/PI,RAD
# COMMON/ELIPSOID/A,F
# DATA EPS/0.5D-13/
# R=1.D0-F
# TU1=R*DSIN(GLAT1)/DCOS(GLAT1)
# TU2=R*DSIN(GLAT2)/DCOS(GLAT2)
# CU1=1./DSQRT(TU1*TU1+1.)
# SU1=CU1*TU1
# CU2=1./DSQRT(TU2*TU2+1.)
# S=CU1*CU2
# BAZ=S*TU2
# FAZ=BAZ*TU1
# X=GLON2-GLON1
# 100 SX=DSIN(X)
# CX=DCOS(X)
# TU1=CU2*SX
# TU2=BAZ-SU1*CU2*CX
# SY=DSQRT(TU1*TU1+TU2*TU2)
# CY=S*CX+FAZ
# Y=DATAN2(SY,CY)
# SA=S*SX/SY
# C2A=-SA*SA+1.
# CZ=FAZ+FAZ
# IF(C2A.GT.0.)CZ=-CZ/C2A+CY
# E=CZ*CZ*2.-1.
# C=((-3.*C2A+4.)*F+4.)*C2A*F/16.
# D=X
# X=((E*CY*C+CZ)*SY*C+Y)*SA
# X=(1.-C)*X*F+GLON2-GLON1
# IF(DABS(D-X).GT.EPS) GO TO 100
# FAZ=DATAN2(TU1,TU2)
# BAZ=DATAN2(CU1*SX,BAZ*CX-SU1*CU2)+PI
# X=DSQRT((1./R/R-1.)*C2A+1.)+1.
# X=(X-2.)/X
# C=1.-X
# C=(X*X/4.+1.)/C
# D=(0.375*X*X-1.)*X
# X=E*CY
# S=1.-E-E
# S=((((SY*SY*4.-3.)*S*CZ*D/6.-X)*D/4.+CZ)*SY*D+Y)*C*A*R
# RETURN
# END
#
if ($_REQUEST["SOURCE"] != "" || $_REQUEST["source"] != "") {
printf("\n<FONT SIZE=1>\n");
# highlight_file("geodinver.php");
show_source("geodinver.php");
printf("</FONT><P><HR><P>\n");
include("../include/sign-oh2mqk.inc");
printf("<P></TR></TABLE>\n");
include("../include/base.inc");
include("../include/foot.inc");
exit;
}
?>
Matti Aarnio <matti.aarnio@zmailer.org>; OH2MQK
