FUNCTION AngNormalize, x ;+ ; NAME: ; AngNormalize ; ; PURPOSE: ; Reduce angle to range [-180, 180) ;- IF x GE 180D THEN x = x - 360D IF x LT -180D THEN x = x + 360D RETURN, x END FUNCTION SUMX, u, v ;+ ; NAME: ; SUMX ; ; PURPOSE: ; Produces an error free sum. ; SUM = SUMX(u, v) returns the rounded sum U+V in SUM[0] and the error in ; SUM[1], such that SUM[0]+SUM[1]=U+V exactly. ;- s = u + v up = s - v vpp = s - up up = up - u vpp = vpp - v t = -(up + vpp) SUM = [s, t] RETURN, SUM END FUNCTION AngDiff, x, y ;+ ; NAME: ; AngDiff ; ; PURPOSE: ; Computes Y-X and reduces the result to (-180,180]. ; X and Y must be in [-180,180]. ;- SUM = SUMX(-x, y) IF SUM[0] + SUM[1] LE 180D AND SUM[0] + SUM[1] GT -180D THEN BEGIN diff = SUM[0] + SUM[1] ENDIF ELSE BEGIN IF (SUM[0] - 180D) + SUM[1] GT 0D THEN BEGIN diff = (SUM[0] - 360D) + SUM[1] ENDIF IF (SUM[0] + 180D) + SUM[1] LE 0D THEN BEGIN diff = (SUM[0] + 360D) + SUM[1] ENDIF ENDELSE RETURN, diff END FUNCTION I1, sigma, epsilon ;+ ; NAME: ; I1 ; ; PURPOSE: ; To solve the Fourier series in equation 15, using equations 17 & 18 ;- A1 = (1D $ + (1D/4D) * epsilon^2D $ + (1D/64D) * epsilon^4D $ + (1D/256D) * epsilon^6D) / (1D - epsilon) ; Eqn 17 C11 = - (1D/2D) * epsilon $ + (3D/16D) * epsilon^3D $ - (1D/32D) * epsilon^5D C12 = - (1D/16D) * epsilon^2D $ + (1D/32D) * epsilon^4D $ - (9D/2048D) * epsilon^6D C13 = - (1D/48D) * epsilon^3D $ + (3D/256D) * epsilon^5D C14 = - (5D/512D) * epsilon^4D $ + (3D/512D) * epsilon^6D C15 = - (7D/1280D) * epsilon^5D ; Eqn 18 RETURN, A1*(sigma + C11*SIN(2D*sigma) $ + C12*SIN(4D*sigma) $ + C13*SIN(6D*sigma) $ + C14*SIN(8D*sigma) $ + C15*SIN(10D*sigma)) ; Eqn 15 END FUNCTION I2, sigma, epsilon ;+ ; NAME: ; I2 ; ; PURPOSE: ; To solve the Fourier series in equation 41, using equations 42 & 43 ;- A2 = (1D $ + (1D/4D) * epsilon^2D $ + (9D/64D) * epsilon^4D $ + (25D/256D) * epsilon^6D) * (1D - epsilon) ; Eqn 42 C21 = (1D/2D) * epsilon $ + (1D/16D) * epsilon^3D $ + (1D/32D) * epsilon^5D C22 = (3D/16D) * epsilon^2D $ + (1D/32D) * epsilon^4D $ + (35D/2048D) * epsilon^6D C23 = (5D/48D) * epsilon^3D $ + (5D/256D) * epsilon^5D C24 = + (35D/512D) * epsilon^4D $ + (7D/512D) * epsilon^6D C25 = (63D/1280D) * epsilon^5D C26 = (77D/2048D) * epsilon^6D ; Eqn 43 RETURN, A2*(sigma + C21*SIN(2D*sigma) $ + C22*SIN(4D*sigma) $ + C23*SIN(6D*sigma) $ + C24*SIN(8D*sigma) $ + C25*SIN(10D*sigma)) ; Eqn 41 END FUNCTION I3, sigma, epsilon, n ;+ ; NAME: ; I3 ; ; PURPOSE: ; To solve the Fourier series in equation 23, using equations 24 & 25 ;- A3 = 1D $ - ((1D/2D) - (1D/2D)*n) * epsilon $ - ((1D/4D) + (1D/8D)*n - (3D/8D)*n^2D) * epsilon^2D $ - ((1D/16D) + (3D/16D)*n + (1D/16D)*n^2D) * epsilon^3D $ - ((3D/64D) + (1D/32D)*n) * epsilon^4D $ - (3D/128D) * epsilon^5D C31 = ((1D/4D) - (1D/4D)*n) * epsilon $ + ((1D/8D) - (1D/8D)*n^2D) * epsilon^2D $ + ((3D/64D) + (3D/64D)*n - (1D/64D)*n^2D) * epsilon^3D $ + ((5D/128D)+(1D/64D)*n) * epsilon^4D $ + (3D/128D) * epsilon^5D C32 = ((1D/16D) - (3D/32D)*n + (1D/32D)*n^2D) * epsilon^2D $ + ((3D/64D) - (1D/32D)*n - (3D/64D)*n^2D) * epsilon^3D $ + ((3D/128D) + (1D/128D)*n) * epsilon^4D $ + (5D/256D) * epsilon^5D C33 = ((5D/192D) - (3D/64D)*n + (5D/192D)*n^2D) * epsilon^3D $ + ((3D/128D) - (5D/192D)*n) * epsilon^4D $ + (7D/512D) * epsilon^5D C34 = ((7D/512D) - (7D/256D)*n) * epsilon^4D $ + (7D/512D) * epsilon^5D C35 = (21D/2560D) * epsilon^5D RETURN, A3*(sigma + C31*SIN(2D*sigma) $ + C32*SIN(4D*sigma) $ + C33*SIN(6D*sigma) $ + C34*SIN(8D*sigma) $ + C35*SIN(10D*sigma)) ; Eqn 23 END ;;;; ;+ ; NAME: ; inverse_geodesic ; ; PURPOSE: ; Find the distance and azimuths between two points on ; WGS84 ellipsoid. ; ; INPUTS: ; lat1, lng1 = Point 1 geodetic long and lat (degrees) ; lat2, lng2 = Point 2 geodetic long and lat (degrees) ; For the most accurate results enter the inputs as double precision numbers. ; baz = boolean (optional) set if you require the heading at the ; second point to be the back azimuth rather than the default ; forawrd azimuth. ; ; OUTPUTS: ; dist = Distance along surface of ellipsoid (metres) ; azi1 = forward geodetic azimuth (pt1 to pt2, degrees) ; azi1 = reverse geodetic azimuth (pt2 to pt1, degrees) ; ; NOTES: ; This algorithm follows the arguments of Charles F. F. Karney's ; 2013 paper 'Algorithms for geodesics' (J Geod (2013) 87:43-55 ; DOI 10.1007/s00190-012-0578-z). ; Karney uses phi to represent latitudes, and lambda for longs. ; ; AUTHORS: ; Sean Elvidge ; E-mail: sean@seanelvidge.com ; http://seanelvidge.com ; ; Chris Mannix ; E-mail: c.mannix@poynting.bham.ac.uk ; http://chrismannix.com ; ; COPYRIGHT ; Copyright (c) 2013, Poynting Research Institute, University of Birmingham ; ; MODIFICATION HISTORY: ; Written by Elvidge & Mannix, March 2013 ;- PRO INVERSE_GEODESIC, lat1, lon1, lat2, lon2, dist, azi1, azi2, BAZ=BAZ IF lat1 EQ lat2 AND lon1 EQ lon2 THEN BEGIN ; If same two points are entered, return 0's dist = 0D azi1 = 0D azi2 = 0D ENDIF ELSE BEGIN ; Default is to return the forward azimuth at point 2 IF N_ELEMENTS(BAZ) EQ 0 THEN baz = 0 ;-------- WGS84 constants (www.wgs84.com) -------------- f = 1.D0/298.257223563D0 ; Flattening factor a = 6378137.0D ; Equatorial radius b = 6356752.314245D ; Polar semi-axis e2=0.00669437999014132D ; Eccentricity squared ep2 = 0.00673949674227643D ; Second eccentricity squared n = 0.00167922038638370D ; Third flattening ;--------------------------------------------------------- dtor = !DPi/180D ; Used for converting from degrees to radians (more precise than !DtoR) eps = 0.5D-12 ; Tolerance of the iteration ; Put the points in canonical config (Eqn 44) l120 = AngDiff(AngNormalize(lon1),AngNormalize(lon2)) ; Make longitude difference positive IF l120 GE 0D THEN lonsign = 1D ELSE lonsign = -1D l120 = l120 * lonsign ; Swap points so that point with higher (absolute) latitude is point 1 IF ABS(lat1) GE ABS(lat2) THEN swapp = 1D ELSE swapp = -1D IF swapp LT 0D THEN BEGIN lonsign = lonsign * (-1D) SWAP, lat1, lat2 ENDIF ; Make lat1 <= 0 IF lat1 LT 0D THEN latsign = 1D ELSE latsign = -1D lat1 = lat1 * latsign lat2 = lat2 * latsign ; Convert lat, lons and l12 to radians glon1 = lon1*dtor glat1 = lat1*dtor glon2 = lon2*dtor glat2 = lat2*dtor ; We adopt the following notation here, throughout the programme: ; tb1 stands for tan(b1) and, similarly, sb1 is sin(b1). ; b1 and b2 represent beta 1 and beta 2 from Karney tb1 = (1D - f)*TAN(glat1) ; Eqn 6 tb2 = (1D - f)*TAN(glat2) ; Eqn 6 cb1 = 1D/SQRT(1D + tb1*tb1) cb2 = 1D/SQRT(1D + tb2*tb2) sb1 = cb1*tb1 sb2 = cb2*tb2 ; Check if the two points are "nearly" antipodal IF ABS(lat2+lat1) LT 0.5D AND ABS(ABS(lon1)+ABS(lon2)-180D) LT 0.5D THEN antipodal=1 ELSE antipodal=-1 ;IF ABS(ABS(lon1)+ABS(lon2)-180D) GT 179D OR ABS(ABS(lon1)+ABS(lon2)-180D) LT 1D THEN antipodal=1 IF l120 GT 179D THEN antipodal=1 ; Make an initial guess of alpha1, the initial guess depends ; on whether the points are nearly antipodal or not. IF antipodal EQ -1 THEN BEGIN ; If the points are not antipodal we proceed as in the beginning of ; Section 5 "Starting point for Newton's method" l120=l120*dtor ; Convert l120 to radians omegabar = SQRT(1D - e2*((cb1+cb2)/2D)^2D) ; Eqn 48 omega12 = l120/omegabar ; Just before Eqn 49 z1 = COMPLEX((cb1*sb2 - sb1*cb2*COS(omega12)),cb2*SIN(omega12),/DOUBLE) z2 = COMPLEX((-sb1*cb2 + cb1*sb2*COS(omega12)),cb1*SIN(omega12),/DOUBLE) sigma12 = ATAN(ABS(z1),(sb1*sb2+cb1*cb2*COS(omega12)),/PHASE) ; Eqn 51 alpha1 = ATAN(z1,/PHASE) ; Eqn 49 alpha2 = ATAN(z2,/PHASE) ; Eqn 50 ENDIF ELSE BEGIN ; If the points are nearly antipodal we run the following: l120=l120*dtor ; lambda(1,2) in radians del = f*a*!DPi*cb1*cb1 ; Just above Eqn 53 x = ((l120)-!DPi)*a*cb1/del ; Eqn 53 y = (ATAN(tb1)+ATAN(tb2))*a/del ; Eqn 53 IF y NE 0 THEN BEGIN coeffs = [-Y^2D, -2D*y^2D, (1D - x^2D - y^2D), 2D, 1D] ; Eqn 55 roots = FZ_ROOTS(coeffs,/DOUBLE, EPS = eps) posroot = WHERE(REAL_PART(roots) GT 0D) posroot = REAL_PART(roots[posroot]) alpha1 = ATAN(COMPLEX(y/posroot,-x/(1D + posroot),/DOUBLE), $ /PHASE) ; Eqn 56 ENDIF ELSE BEGIN IF x GT 0 THEN neg = -1D ELSE neg = 1D alpha1 = ATAN(COMPLEX(SQRT(MAX([0D,(1D - x^2D)]))*neg,-x, $ /DOUBLE),/PHASE) ; Eqn 57 ENDELSE ENDELSE dell12 = 1 ; To start the while loop off loop = 1 WHILE ABS(dell12) GT eps AND loop LE 30 DO BEGIN ;---- SOLVE TRIANGLE NEA ------- ca1 = COS(alpha1) sa1 = SIN(alpha1) alpha0 = ATAN(COMPLEX(ABS(COMPLEX(ca1,sa1*sb1,/DOUBLE)),sa1*cb1, $ /DOUBLE),/PHASE) ; Eqn 10 ca0 = COS(alpha0) sa0 = SIN(alpha0) sigma1 = ATAN(COMPLEX(ca1*cb1,sb1,/DOUBLE),/PHASE) ; Eqn 11 cs1 = COS(sigma1) ss1 = SIN(sigma1) omega1 = ATAN(COMPLEX(cs1,sa0*ss1,/DOUBLE),/PHASE) ; Eqn 12 ;---- SOLVE TRIANGLE NEB ------- ca2 = SQRT(MAX([(ca1*ca1*cb1*cb1 + cb2*cb2 - cb1*cb1),0D])) / cb2 ; Eqn 45 sigma2 = ATAN(COMPLEX(ca2*cb2,sb2,/DOUBLE),/PHASE) ; Eqn 11 cs2 = COS(sigma2) ss2 = SIN(sigma2) omega2 = ATAN(COMPLEX(cs2,sa0*ss2,/DOUBLE),/PHASE) ; Eqn 12 ;---- DETERMINE l12 ------- k = SQRT(ep2)*ca0 ; Eqn 9 k2 = k^2D epsilon = (SQRT(1D + k2)-1D) / (SQRT(1D + k2) + 1D) ; To find l1 (lambda 1) in Eqn 8 we perform the integration using ; Eqns 23-25, which are built as seperate functions I1, I2 and I3. I3s1 = I3(sigma1,epsilon,n) l1 = omega1 - f*sa0*I3s1 ; Eqn 8 I3s2 = I3(sigma2,epsilon,n) l2 = omega2 - f*sa0*I3s2 ; Eqn 8 ;l12 = AngDiff(l1,l2) ; New l12 value l12 = AngDiff(AngNormalize(l1),AngNormalize(l2)) ; Make longitude difference positive IF l12 GE 0D THEN lonsignn = 1D ELSE lonsignn = -1D l12 = l12 * lonsignn ;---- UPDATE alpha1 ------- dell12 = l12-l120 ; Difference between calced. l12 and original value I1s1 = I1(sigma1,epsilon) I1s2 = I1(sigma2,epsilon) I2s1 = I2(sigma1,epsilon) I2s2 = I2(sigma2,epsilon) Js1 = I1S1 - I2S1 ; Eqn 40 Js2 = I1S2 - I2S2 ; Eqn 40 m12 = b*(SQRT(1D + k2*ss2*ss2)*cs1*ss2 - $ SQRT(1D + k2*ss1*ss1)*ss1*cs2 - $ cs1*cs2*(Js2-Js1)) ; Eqn 38 dl12_da1 = (m12/a) / (ca2*cb2) ; Eqn 46 delalpha1 = -dell12/dl12_da1 ; Update the value of alpha1 alpha1 = alpha1 + delalpha1 ; If dell12 is greater than eps perform the iteration again ; typically only 2 to 4 iterations are required, in a small ; number of cases up to 16 may be required. ; Loop forces iterations to be less than 30. In a small number of ; cases dell12 remains slightly above eps indefinitely. loop = loop + 1 ENDWHILE ;---- SOLVE TRIANGLE NEA ------- ca1 = COS(alpha1) sa1 = SIN(alpha1) alpha0 = ATAN(COMPLEX(ABS(COMPLEX(ca1,sa1*sb1,/DOUBLE)),sa1*cb1, $ /DOUBLE),/PHASE) ; Eqn 10 ca0 = COS(alpha0) sa0 = SIN(alpha0) sigma1 = ATAN(COMPLEX(ca1*cb1,sb1,/DOUBLE),/PHASE) ; Eqn 11 cs1 = COS(sigma1) ss1 = SIN(sigma1) omega1 = ATAN(COMPLEX(cs1,sa0*ss1,/DOUBLE),/PHASE) ; Eqn 12 ;---- SOLVE TRIANGLE NEB ------- ca2 = SQRT(ca1*ca1*cb1*cb1 + cb2*cb2 - cb1*cb1) / cb2 ; Eqn 45 sa2 = sa0/cb2 sigma2 = ATAN(COMPLEX(ca2*cb2,sb2,/DOUBLE),/PHASE) ; Eqn 11 cs2 = COS(sigma2) ss2 = SIN(sigma2) omega2 = ATAN(COMPLEX(cs2,sa0*ss2,/DOUBLE),/PHASE) ; Eqn 12 k = SQRT(ep2)*ca0 ; Eqn 9 k2 = k^2D epsilon = (SQRT(1D + k2)-1D) / (SQRT(1D + k2) + 1D) s1 = b * I1(sigma1,epsilon) s2 = b * I1(sigma2,epsilon) s12 = s2-s1 I3s1 = I3(sigma1,epsilon,n) l1 = omega1 - f*sa0*I3s1 ; Eqn 8 I3s2 = I3(sigma2,epsilon,n) l2 = omega2 - f*sa0*I3s2 ; Eqn 8 l12 = l2 - l1 ; New l12 value ; Set up output values dist = s12 IF swapp LT 0D THEN BEGIN SWAP, sa1, sa2 SWAP, ca1, ca2 ENDIF sa1 = sa1 * swapp * lonsign sa2 = sa2 * swapp * lonsign ca1 = ca1 * swapp * latsign ca2 = ca2 * swapp * latsign azi1 = 0D - ATAN(-sa1, ca1, /PHASE)/dtor IF BAZ EQ 0 THEN BEGIN azi2 = 0D - ATAN(-sa2, ca2, /PHASE)/dtor ENDIF ELSE BEGIN azi2 = (0D - ATAN(-sa2, ca2, /PHASE)/dtor) - 180D ENDELSE ; azi1 and azi2 should be bearings, so we want them to be between ; 0 and 360 degrees IF azi1 LT 0D THEN azi1 = azi1 + 360D IF azi2 LT 0D THEN azi2 = azi2 + 360D ENDELSE END ;;;;;;