import * as L from "leaflet" ; export interface GeodesicOptions extends L . PolylineOptions { wrap ? : boolean ; steps ? : number ; radius ? : number ; } export interface WGS84Vector extends L . LatLngLiteral { bearing : number ; } export interface GeoDistance { distance : number ; initialBearing : number ; finalBearing : number ; } export class GeodesicCore { readonly options : GeodesicOptions = { wrap : true , steps : 3 } ; readonly ellipsoid = { a : 6378137 , b : 6356752.3142 , f : 1 / 298.257223563 } ; constructor ( options ? ) { this . options = { ... this . options , ... options } ; } toRadians ( degree ) { return degree * Math . PI / 180 ; } toDegrees ( radians ) { return radians * 180 / Math . PI ; } param n param p return mod ( n , p ) { const r = n % p ; return r < 0 ? r + p : r ; } param degrees return wrap360 ( degrees ) { if ( 0 <= degrees && degrees < 360 ) { return degrees ; } else { return this . mod ( degrees , 360 ) ; } } param degrees param max return wrap ( degrees , max = 360 ) { if ( - max <= degrees && degrees <= max ) { return degrees ; } else { return this . mod ( ( degrees + max ) , 2 * max ) - max ; } } param start param bearing param distance param maxInterations return direct ( start , bearing , distance , maxInterations = 100 ) { const φ1 = this . toRadians ( start . lat ) ; const λ1 = this . toRadians ( start . lng ) ; const α1 = this . toRadians ( bearing ) ; const s = distance ; const ε = Number . EPSILON * 1000 ; const { a , b , f } = this . ellipsoid ; const sinα1 = Math . sin ( α1 ) ; const cosα1 = Math . cos ( α1 ) ; const tanU1 = ( 1 - f ) * Math . tan ( φ1 ) , cosU1 = 1 / Math . sqrt ( ( 1 + tanU1 * tanU1 ) ) , sinU1 = tanU1 * cosU1 ; const σ1 = Math . atan2 ( tanU1 , cosα1 ) ; const sinα = cosU1 * sinα1 ; const cosSqα = 1 - sinα * sinα ; const uSq = cosSqα * ( a * a - b * b ) / ( b * b ) ; const A = 1 + uSq / 16384 * ( 4096 + uSq * ( - 768 + uSq * ( 320 - 175 * uSq ) ) ) ; const B = uSq / 1024 * ( 256 + uSq * ( - 128 + uSq * ( 74 - 47 * uSq ) ) ) ; let σ = s / ( b * A ) , sinσ = null , cosσ = null , Δσ = null ; let cos2σₘ = null ; let σʹ = null , iterations = 0 ; do { cos2σₘ = Math . cos ( 2 * σ1 + σ ) ; sinσ = Math . sin ( σ ) ; cosσ = Math . cos ( σ ) ; Δσ = B * sinσ * ( cos2σₘ + B / 4 * ( cosσ * ( - 1 + 2 * cos2σₘ * cos2σₘ ) - B / 6 * cos2σₘ * ( - 3 + 4 * sinσ * sinσ ) * ( - 3 + 4 * cos2σₘ * cos2σₘ ) ) ) ; σʹ = σ ; σ = s / ( b * A ) + Δσ ; } while ( Math . abs ( σ - σʹ ) > ε && ++ iterations < maxInterations ) ; if ( iterations >= maxInterations ) { throw new EvalError ( `Direct vincenty formula failed to converge after ${ maxInterations } iterations 
                (start=${ start . lat }/${ start . lng }; bearing=${ bearing }; distance=${ distance })` ) ; } const x = sinU1 * sinσ - cosU1 * cosσ * cosα1 ; const φ2 = Math . atan2 ( sinU1 * cosσ + cosU1 * sinσ * cosα1 , ( 1 - f ) * Math . sqrt ( sinα * sinα + x * x ) ) ; const λ = Math . atan2 ( sinσ * sinα1 , cosU1 * cosσ - sinU1 * sinσ * cosα1 ) ; const C = f / 16 * cosSqα * ( 4 + f * ( 4 - 3 * cosSqα ) ) ; const dL = λ - ( 1 - C ) * f * sinα * ( σ + C * sinσ * ( cos2σₘ + C * cosσ * ( - 1 + 2 * cos2σₘ * cos2σₘ ) ) ) ; const λ2 = λ1 + dL ; const α2 = Math . atan2 ( sinα , - x ) ; return { lat : this . toDegrees ( φ2 ) , lng : this . toDegrees ( λ2 ) , bearing : this . wrap360 ( this . toDegrees ( α2 ) ) } ; } param start param dest return inverse ( start , dest , maxInterations = 100 , mitigateConvergenceError = true ) { const p1 = start , p2 = dest ; const φ1 = this . toRadians ( p1 . lat ) , λ1 = this . toRadians ( p1 . lng ) ; const φ2 = this . toRadians ( p2 . lat ) , λ2 = this . toRadians ( p2 . lng ) ; const π = Math . PI ; const ε = Number . EPSILON ; const { a , b , f } = this . ellipsoid ; const dL = λ2 - λ1 ; const tanU1 = ( 1 - f ) * Math . tan ( φ1 ) , cosU1 = 1 / Math . sqrt ( ( 1 + tanU1 * tanU1 ) ) , sinU1 = tanU1 * cosU1 ; const tanU2 = ( 1 - f ) * Math . tan ( φ2 ) , cosU2 = 1 / Math . sqrt ( ( 1 + tanU2 * tanU2 ) ) , sinU2 = tanU2 * cosU2 ; const antipodal = Math . abs ( dL ) > π / 2 || Math . abs ( φ2 - φ1 ) > π / 2 ; let λ = dL , sinλ = null , cosλ = null ; let σ = antipodal ? π : 0 , sinσ = 0 , cosσ = antipodal ? - 1 : 1 , sinSqσ = null ; let cos2σₘ = 1 ; let sinα = null , cosSqα = 1 ; let C = null ; let λʹ = null , iterations = 0 ; do { sinλ = Math . sin ( λ ) ; cosλ = Math . cos ( λ ) ; sinSqσ = ( cosU2 * sinλ ) * ( cosU2 * sinλ ) + ( cosU1 * sinU2 - sinU1 * cosU2 * cosλ ) * ( cosU1 * sinU2 - sinU1 * cosU2 * cosλ ) ; if ( Math . abs ( sinSqσ ) < ε ) { break ; } sinσ = Math . sqrt ( sinSqσ ) ; cosσ = sinU1 * sinU2 + cosU1 * cosU2 * cosλ ; σ = Math . atan2 ( sinσ , cosσ ) ; sinα = cosU1 * cosU2 * sinλ / sinσ ; cosSqα = 1 - sinα * sinα ; cos2σₘ = ( cosSqα !== 0 ) ? ( cosσ - 2 * sinU1 * sinU2 / cosSqα ) : 0 ; C = f / 16 * cosSqα * ( 4 + f * ( 4 - 3 * cosSqα ) ) ; λʹ = λ ; λ = dL + ( 1 - C ) * f * sinα * ( σ + C * sinσ * ( cos2σₘ + C * cosσ * ( - 1 + 2 * cos2σₘ * cos2σₘ ) ) ) ; const iterationCheck = antipodal ? Math . abs ( λ ) - π : Math . abs ( λ ) ; if ( iterationCheck > π ) { throw new EvalError ( 'λ > π' ) ; } } while ( Math . abs ( λ - λʹ ) > 1e-12 && ++ iterations < maxInterations ) ; if ( iterations >= maxInterations ) { if ( mitigateConvergenceError ) { return this . inverse ( start , new L . LatLng ( dest . lat , dest . lng - 0.01 ) , maxInterations , mitigateConvergenceError ) ; } else { throw new EvalError ( `Inverse vincenty formula failed to converge after ${ maxInterations } iterations 
                    (start=${ start . lat }/${ start . lng }; dest=${ dest . lat }/${ dest . lng })` ) ; } } const uSq = cosSqα * ( a * a - b * b ) / ( b * b ) ; const A = 1 + uSq / 16384 * ( 4096 + uSq * ( - 768 + uSq * ( 320 - 175 * uSq ) ) ) ; const B = uSq / 1024 * ( 256 + uSq * ( - 128 + uSq * ( 74 - 47 * uSq ) ) ) ; const Δσ = B * sinσ * ( cos2σₘ + B / 4 * ( cosσ * ( - 1 + 2 * cos2σₘ * cos2σₘ ) - B / 6 * cos2σₘ * ( - 3 + 4 * sinσ * sinσ ) * ( - 3 + 4 * cos2σₘ * cos2σₘ ) ) ) ; const s = b * A * ( σ - Δσ ) ; const α1 = Math . abs ( sinSqσ ) < ε ? 0 : Math . atan2 ( cosU2 * sinλ , cosU1 * sinU2 - sinU1 * cosU2 * cosλ ) ; const α2 = Math . abs ( sinSqσ ) < ε ? π : Math . atan2 ( cosU1 * sinλ , - sinU1 * cosU2 + cosU1 * sinU2 * cosλ ) ; return { distance : s , initialBearing : Math . abs ( s ) < ε ? NaN : this . wrap360 ( this . toDegrees ( α1 ) ) , finalBearing : Math . abs ( s ) < ε ? NaN : this . wrap360 ( this . toDegrees ( α2 ) ) } as GeoDistance ; } param firstPos param firstBearing param secondPos param secondBearing intersection ( firstPos , firstBearing , secondPos , secondBearing ) { const φ1 = this . toRadians ( firstPos . lat ) ; const λ1 = this . toRadians ( firstPos . lng ) ; const φ2 = this . toRadians ( secondPos . lat ) ; const λ2 = this . toRadians ( secondPos . lng ) ; const θ13 = this . toRadians ( firstBearing ) ; const θ23 = this . toRadians ( secondBearing ) ; const Δφ = φ2 - φ1 , Δλ = λ2 - λ1 ; const π = Math . PI ; const ε = Number . EPSILON ; const δ12 = 2 * Math . asin ( Math . sqrt ( Math . sin ( Δφ / 2 ) * Math . sin ( Δφ / 2 ) + Math . cos ( φ1 ) * Math . cos ( φ2 ) * Math . sin ( Δλ / 2 ) * Math . sin ( Δλ / 2 ) ) ) ; if ( Math . abs ( δ12 ) < ε ) { return firstPos ; } const cosθa = ( Math . sin ( φ2 ) - Math . sin ( φ1 ) * Math . cos ( δ12 ) ) / ( Math . sin ( δ12 ) * Math . cos ( φ1 ) ) ; const cosθb = ( Math . sin ( φ1 ) - Math . sin ( φ2 ) * Math . cos ( δ12 ) ) / ( Math . sin ( δ12 ) * Math . cos ( φ2 ) ) ; const θa = Math . acos ( Math . min ( Math . max ( cosθa , - 1 ) , 1 ) ) ; const θb = Math . acos ( Math . min ( Math . max ( cosθb , - 1 ) , 1 ) ) ; const θ12 = Math . sin ( λ2 - λ1 ) > 0 ? θa : 2 * π - θa ; const θ21 = Math . sin ( λ2 - λ1 ) > 0 ? 2 * π - θb : θb ; const α1 = θ13 - θ12 ; const α2 = θ21 - θ23 ; if ( Math . sin ( α1 ) === 0 && Math . sin ( α2 ) === 0 ) { return null ; } if ( Math . sin ( α1 ) * Math . sin ( α2 ) < 0 ) { return null ; } const cosα3 = - Math . cos ( α1 ) * Math . cos ( α2 ) + Math . sin ( α1 ) * Math . sin ( α2 ) * Math . cos ( δ12 ) ; const δ13 = Math . atan2 ( Math . sin ( δ12 ) * Math . sin ( α1 ) * Math . sin ( α2 ) , Math . cos ( α2 ) + Math . cos ( α1 ) * cosα3 ) ; const φ3 = Math . asin ( Math . min ( Math . max ( Math . sin ( φ1 ) * Math . cos ( δ13 ) + Math . cos ( φ1 ) * Math . sin ( δ13 ) * Math . cos ( θ13 ) , - 1 ) , 1 ) ) ; const Δλ13 = Math . atan2 ( Math . sin ( θ13 ) * Math . sin ( δ13 ) * Math . cos ( φ1 ) , Math . cos ( δ13 ) - Math . sin ( φ1 ) * Math . sin ( φ3 ) ) ; const λ3 = λ1 + Δλ13 ; return new L . LatLng ( this . toDegrees ( φ3 ) , this . toDegrees ( λ3 ) ) ; } midpoint ( start , dest ) { const φ1 = this . toRadians ( start . lat ) ; const λ1 = this . toRadians ( start . lng ) ; const φ2 = this . toRadians ( dest . lat ) ; const Δλ = this . toRadians ( dest . lng - start . lng ) ; const A = { x : Math . cos ( φ1 ) , y : 0 , z : Math . sin ( φ1 ) } ; const B = { x : Math . cos ( φ2 ) * Math . cos ( Δλ ) , y : Math . cos ( φ2 ) * Math . sin ( Δλ ) , z : Math . sin ( φ2 ) } ; const C = { x : A . x + B . x , y : A . y + B . y , z : A . z + B . z } ; const φm = Math . atan2 ( C . z , Math . sqrt ( C . x * C . x + C . y * C . y ) ) ; const λm = λ1 + Math . atan2 ( C . y , C . x ) ; return new L . LatLng ( this . toDegrees ( φm ) , this . toDegrees ( λm ) ) ; } }