open System
module AzimuthalEquidistantProjection =

    let inline degToRad d = 0.0174532925199433 * d; // (1.0/180.0 * Math.PI) * d

    let project centerlon centerlat lon lat =
        // http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
        // http://www.radicalcartography.net/?projectionref
        let t:float = degToRad lat
        let l:float = degToRad lon
        let t1 = degToRad centerlat // latitude center of projection
        let l0 = degToRad centerlon // longitude center of projection
        let c = Math.Acos ((sin t1) * (sin t) + (cos t1) * (cos t) * (cos (l-l0)))
        let k = c / (sin c)
        let x = k * (cos t) * (sin (l-l0))
        let y = k * (cos t1) * (sin t) - (sin t1) * (cos t) * (cos (l-l0))
        (x, y)

    let project_optimized l0 t1 cost1 lon lat =
        let t:float = degToRad lat
        let l:float = degToRad lon

        let costcosll0 = (cos t) * (cos (l-l0))
        let sint = sin t
        let sint1 = sin t1
        let c = Math.Acos ((sint1) * (sint) + (cost1) * costcosll0)
        let k = c / (sin c)
        let x = k * (cos t) * (sin (l-l0))
        let y = k * (cost1) * (sint) - (sint1) * costcosll0
        (x, y)

    let buildProjection centerLon centerLat = 
        let t1 = degToRad centerLat
        let l0 = degToRad centerLon
        let cost1 = cos t1
        project_optimized l0 t1 cost1