UTM class

# File lib/ravn/geo/utm.rb, line 58
def self::from_position( pos )
        zone = ( (pos.long + 180).to_f / 6 ).floor + 1
        λ0   = ( ( zone - 1 ) * 6 - 180 + 3 ).to_radians
        band = MGRS_LAT_BANDS[ ( pos.lat.to_f / 8 + 10 ).floor ]

        adjust_inc = -> { zone += 1 && λ0 += 6.to_radians }
        adjust_dec = -> { zone -= 1 && λ0 -= 6.to_radians }

        #
        # Zone exceptions / adjustments
        #

        # Norway
        if zone == 31 && band == 'V' && long >= 3
                adjust_inc.call
        end

        # Svalbard
        if band == 'X'
                if zone == 32
                        long >= 9 ? adjust_inc.call : adjust_dec.call
                end
                if zone == 34
                        long >= 21 ? adjust_inc.call : adjust_dec.call
                end
                if zone == 36
                        long >= 33? adjust_inc.call : adjust_dec.call
                end
        end

        φ = pos.lat.to_radians
        λ = pos.long.to_radians - λ0

        a, f = WGS84.values_at( :a, :f )

        e  = Math.sqrt( f * ( 2 - f ) ) # eccentricity
        n  = f / ( 2 - f ) # 3rd flattening
        n2 = n*n
        n3 = n*n2
        n4 = n*n3
        n5 = n*n4
        n6 = n*n5

        cosλ = Math.cos( λ )
        sinλ = Math.sin( λ )
        tanλ = Math.tan( λ )

        τ = Math.tan( φ )
        σ = Math.sinh( e * Math.atanh( e*τ / Math.sqrt( 1+τ*τ ) ) )

        τʹ = τ * Math.sqrt(1 + σ * σ) - σ * Math.sqrt( 1 + τ * τ)

        ξʹ = Math.atan2( τʹ, cosλ )
        ηʹ = Math.asinh( sinλ / Math.sqrt( τʹ * τʹ + cosλ * cosλ) )

        a_ = a / (1 + n) * (1 + 1.0/4*n2 + 1.0/64*n4 + 1.0/256*n6) # 2πa_ is the circumference of a meridian

        α = [ nil, # note α is one-based array (6th order Krüger expressions)
         1.0/2*n - 2.0/3*n2 + 5.0/16*n3 + 41.0/180*n4 - 127.0/288*n5 + 7891.0/37800*n6,
         13.0/48*n2 - 3.0/5*n3 + 557.0/1440*n4 + 281.0/630*n5 - 1983433.0/1935360*n6,
         61.0/240*n3 - 103.0/140*n4 + 15061.0/26880*n5 + 167603.0/181440*n6,
         49561.0/161280*n4 - 179.0/168*n5 + 6601661.0/7257600*n6,
         34729.0/80640*n5 - 3418889.0/1995840*n6,
         212378941.0/319334400*n6
        ]

        ξ = ξʹ
        (1..6).each {|j| ξ += α[j] * Math.sin(2*j*ξʹ) * Math.cosh(2*j*ηʹ) }

        η = ηʹ
        (1..6).each {|j| η += α[j] * Math.cos(2*j*ξʹ) * Math.sinh(2*j*ηʹ) }

        x = UTM_SCALE * a_ * η
        y = UTM_SCALE * a_ * ξ

        #
        # Convergence
        #

        pʹ = 1
        (1..6).each {|j| pʹ += 2*j*α[j] * Math.cos(2*j*ξʹ) * Math.cosh(2*j*ηʹ) }
        qʹ = 0
        (1..6).each {|j| qʹ += 2*j*α[j] * Math.sin(2*j*ξʹ) * Math.sinh(2*j*ηʹ) }

        γʹ = Math.atan( τʹ / Math.sqrt(1+τʹ*τʹ)*tanλ )
        γʺ = Math.atan2( qʹ, pʹ )
        γ  = γʹ + γʺ

        #
        # Scale
        #

        sinφ = Math.sin( φ )
        kʹ   = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) / Math.sqrt(τʹ*τʹ + cosλ*cosλ)
        kʺ   = a_ / a * Math.sqrt(pʹ*pʹ + qʹ*qʹ)
        k    = UTM_SCALE * kʹ * kʺ

        # shift x/y to false origins
        x = x + FALSE_EASTING
        y = y + FALSE_NORTHING if y < 0

        # round to reasonable precision
        x = x.round # x = x.round( 9 )
        y = y.round # y = y.round( 9 )
        convergence = γ.to_degrees.round( 9 )
        scale = k.round( 12 )

        # hemisphere
        h = pos.lat >= 0 ? 'N' : 'S'

        return new(
                convergence,
                scale,
                zone:       zone,
                hemisphere: h,
                easting:    x,
                northing:   y,
        )
end

# File lib/ravn/geo/utm.rb, line 201
def initialize( convergence=nil, scale=UTM_SCALE, zone:, hemisphere:, easting:, northing: )
        @convergence = convergence
        @scale       = scale
        @zone        = zone
        @hemisphere  = hemisphere
        @easting     = easting
        @northing    = northing

        raise RangeError, "Invalid zone %p" % [ zone ] if zone <= 1 || zone >= 60
        raise RangeError, "Invalid hemisphere, must be 'N' or 'S'" unless hemisphere =~ /^[NS]$/
end

# File lib/ravn/geo/utm.rb, line 181
def self::parse( utm )
        match = utm.match( MATCHER ) or raise ArgumentError, "Invalid UTM: %p" % [ utm ]
        return new(
                zone:       match[ :zone ].to_i,
                hemisphere: match[ :hemisphere ],
                easting:    match[ :easting ].to_i,
                northing:   match[ :northing ].to_i
        )
end