mapproj.tcl

Go to the documentation of this file.

/*  mapproj.tcl --*/
/* */
/*  Package for map projections.*/
/* */
/*  Copyright (c) 2007 by Kevin B. Kenny.*/
/* */
/*  See the file "license.terms" for information on usage and redistribution*/
/*  of this file, and for a DISCLAIMER OF ALL WARRANTIES.*/
/* */
/*  RCS: @(#) $Id: mapproj.tcl,v 1.1 2007/08/24 22:36:35 kennykb Exp $*/
/* ------------------------------------------------------------------------------*/

package require Tcl 8.4
package require math::interpolate 1.0
package require math::special 0.2.1

package provide mapproj 1.0

/*  ::mapproj --*/
/* */
/*  Namespace holding the procedures and values.*/

namespace ::mapproj {

    /*  Cylindrical projections*/

    namespace export toPlateCarree fromPlateCarree
    namespace export toCassini fromCassini
    namespace export toCylindricalEqualArea fromCylindricalEqualArea
    namespace export toMercator fromMercator

    /*  Named cylindric equal-area projections with specific*/
    /*  standard parallels*/

    namespace export toLambertCylindricalEqualArea \
    fromLambertCylindricalEqualArea
    namespace export toBehrmann fromBehrmann
    namespace export toTrystanEdwards fromTrystanEdwards
    namespace export toHoboDyer fromHoboDyer
    namespace export toGallPeters fromGallPeters
    namespace export toBalthasart fromBalthasart

    /*  Pseudocylindrical projections - equal area*/

    namespace export toSinusoidal fromSinusoidal
    namespace export toMillerCylindrical fromMillerCylindrical
    namespace export toMollweide fromMollweide
    namespace export toEckertIV fromEckertIV toEckertVI fromEckertVI

    /*  Pseudocylindrical projections - compromise*/

    namespace export toRobinson fromRobinson

    /*  Azimuthal projections*/

    namespace export toAzimuthalEquidistant fromAzimuthalEquidistant
    namespace export toLambertAzimuthalEqualArea fromLambertAzimuthalEqualArea
    namespace export toStereographic fromStereographic
    namespace export toOrthographic fromOrthographic
    namespace export toGnomonic fromGnomonic

    /*  Pseudo-azimuthal projections*/

    namespace export toHammer fromHammer

    /*  Conic projections*/

    namespace export toConicEquidistant fromConicEquidistant
    namespace export toLambertConformalConic fromLambertConformalConic
    namespace export toAlbersEqualAreaConic fromAlbersEqualAreaConic

    /*  Miscellaneous projections*/

    namespace export toPeirceQuincuncial fromPeirceQuincuncial

    /*  Fundamental constants*/

    variable pi [expr {acos(-1.0)}]
    variable twopi [expr {2.0 * $pi}]
    variable halfpi [expr {0.5 * $pi}]
    variable quarterpi [expr {0.25 * $pi}]
    variable threequarterpi [expr {0.75 * $pi}]
    variable mquarterpi [expr {-0.25 * $pi}]
    variable mthreequarterpi [expr {-0.75 * $pi}]
    variable radian [expr {180. / $pi}]
    variable degree [expr {$pi / 180.}]
    variable sqrt2 [expr {sqrt(2.)}]
    variable sqrt8 [expr {2. * $sqrt2}]
    variable halfSqrt3 [expr {sqrt(3.) / 2.}]
    variable halfSqrt2 [expr {sqrt(2.) / 2.}]
    variable EckertIVK1 [expr {2.0 / sqrt($pi * (4.0 + $pi))}]
    variable EckertIVK2 [expr {2.0 * sqrt($pi / (4.0 + $pi))}]
    variable EckertVIK1 [expr {sqrt(2.0 + $pi)}]
    variable PeirceQuincuncialScale 3.7081493546027438 ;/*  2*K(1/2)*/
    variable PeirceQuincuncialLimit 1.8540746773013719 ;/*  K(1/2)*/

    /*  Table of parallel length and distance from equator for the*/
    /*  Robinson projection*/

    variable RobinsonLatitude {
    -90.0   -85.0   -80.0   -75.0   -70.0   -65.0
    -60.0   -55.0   -50.0   -45.0   -40.0   -35.0
    -30.0   -25.0   -20.0   -15.0   -10.0   -5.0
    0.0 5.0 10.0    15.0    20.0    25.0    30.0
    35.0    40.0    45.0    50.0    55.0    60.0
    65.0    70.0    75.0    80.0    85.0    90.0
    }
    variable RobinsonPLEN {
    0.5322  0.5722  0.6213  0.6732  0.7186  0.7597
    0.7986  0.8350  0.8679  0.8962  0.9216  0.9427
    0.9600  0.9730  0.9822  0.9900  0.9954  0.9986
    1.0000  0.9986  0.9954  0.9900  0.9822  0.9730  0.9600
    0.9427  0.9216  0.8962  0.8679  0.8350  0.7986
    0.7597  0.7186  0.6732  0.6213  0.5722  0.5322
    }
    variable RobinsonPDFE {
    -1.0000 -0.9761 -0.9394 -0.8936 -0.8435 -0.7903
    -0.7346 -0.6769 -0.6176 -0.5571 -0.4958 -0.4340
    -0.3720 -0.3100 -0.2480 -0.1860 -0.1240 -0.0620
    0.0000  0.0620  0.1240  0.1860  0.2480  0.3100  0.3720
    0.4340  0.4958  0.5571  0.6176  0.6769  0.7346
    0.7903  0.8435  0.8936  0.9394  0.9761  1.0000
    }

    /*  Interpolation tables for Robinson*/

    variable RobinsonSplinePLEN \
    [math::interpolate::prepare-cubic-splines \
         $RobinsonLatitude $RobinsonPLEN]
    variable RobinsonSplinePDFE \
    [math::interpolate::prepare-cubic-splines \
         $RobinsonLatitude $RobinsonPDFE]
    variable RobinsonM [expr {0.5072 * $pi}]

    namespace import ::math::special::cn

}

/*  ::mapproj::ellF - */
/* */
/*  Computes the Legendre incomplete elliptic integral of the*/
/*  first kind:*/
/* */
/*  F(phi, k) = \integral_0^phi dtheta/sqrt(1 - k**2 sin**2 theta)*/
/* */
/* */
/*  Parameters:*/
/*  phi -- Limit of integration; angle around the ellipse*/
/*  k -- Eccentricity*/
/* */
/*  Results:*/
/*  Returns F(phi, k)*/
/* */
/*  Notes:*/
/*  We compute this integral in terms of the Carlson elliptic integral*/
/*  ellRF(x, y, z).*/

ret  ::mapproj::ellF (type phi , type k) {
    return [ellFaux [expr {cos($phi)}] [expr {sin($phi)}] $k]
}

/*  ::mapproj::ellFaux -*/
/* */
/*  Computes the Legendre incomplete elliptic integral of the*/
/*  first kind when circular functions of the 'phi' argument*/
/*  are already available.*/
/* */
/*  Parameters:*/
/*  cos_phi - Cosine of the argument*/
/*  sin_phi - Sine of the argument*/
/*  k - Parameter*/
/* */
/*  Results:*/
/*  Returns F(atan(sin_phi/cos_phi), k)*/

ret  ::mapproj::ellFaux (type cos_, type phi , type sin_, type phi , type k) {
    set rf [ellRF [expr {$cos_phi * $cos_phi}] \
           [expr {1.0 - $k * $k * $sin_phi * $sin_phi}] \
           1.0]
    return [expr {$sin_phi * $rf}]
}

/*  ::mapproj::ellRF --*/
/* */
/*  Computes the Carlson incomplete elliptic integral of the*/
/*  first kind:*/
/* */
/*  RF(x, y, z) = 1/2 * integral_0^inf dt/sqrt((t+x)*(t+y)*(t+z))*/
/* */
/*  Parameters:*/
/*  x, y, z -- Interchangeable parameters of the integral*/
/* */
/*  Results:*/
/*  Returns the value of RF*/

ret  ::mapproj::ellRF (type x , type y , type z) {
    if {$x < 0.0 || $y < 0.0 || $z < 0.0} {
    return -code error "Negative argument to Carlson's ellRF" \
        -errorCode "ellRF negArgument"
    }
    set delx 1.0; set dely 1.0; set delz 1.0
    while {abs($delx) > 0.0025 || abs($dely) > 0.0025 || abs($delz) > 0.0025} {
    set sx [expr {sqrt($x)}]
    set sy [expr {sqrt($y)}]
    set sz [expr {sqrt($z)}]
    set len [expr {$sx * ($sy + $sz) + $sy * $sz}]
    set x [expr {0.25 * ($x + $len)}]
    set y [expr {0.25 * ($y + $len)}]
    set z [expr {0.25 * ($z + $len)}]
    set mean [expr {($x + $y + $z) / 3.0}]
    set delx [expr {($mean - $x) / $mean}]
    set dely [expr {($mean - $y) / $mean}]
    set delz [expr {($mean - $z) / $mean}]
    }
    set e2 [expr {$delx * $dely - $delz * $delz}]
    set e3 [expr {$delx * $dely * $delz}]
    return [expr {(1.0 + ($e2 / 24.0 - 0.1 - 3.0 * $e3 / 44.0) * $e2
           + $e3 / 14.) / sqrt($mean)}]
}

/*  ::mapproj::toPlateCarree --*/
/* */
/*  Project a latitude and longitude onto the plate carrée.*/
/* */
/*  Parameters:*/
/*  phi_0 -- Latitude of the center of the sheet in degrees*/
/*  lambda_0 -- Longitude of the center of sheet in degrees*/
/*  lambda -- Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates.  Units of x and y are Earth radii;*/
/*  scale is true at the Equator.*/

ret  ::mapproj::toPlateCarree (type lambda_0 , type phi_0 , type lambda , type phi) {
    variable degree
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set x [expr {$lambda * $degree}]
    set y [expr {($phi - $phi_0) * $degree}]
    return [list $x $y]
}

/*  ::mapproj::fromPlateCarree --*/
/* */
/*  Solve a plate carrée projection for the*/
/*  latitude and longitude represented by a point on the map.*/
/* */
/*  Parameters:*/
/*  phi_0 -- Latitude of the center of the sheet in degrees*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromPlateCarree (type phi_0 , type lambda_0 , type x , type y) {
    variable radian
    set lambda [expr {$lambda_0 + $x * $radian + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$y * $radian + $phi_0}]
    return [list $lambda $phi]
}

/*  mapproj::toCylindricalEqualArea --*/
/* */
/*  Project a latitude and longitude into cylindrical equal-area*/
/*  co-ordinates.*/
/* */
/*  Parameters:*/
/*  phi_1 --    Standard latitude in degrees*/
/*  phi_0 -- Latitude of the center of the sheet in degrees*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates.  Units of x and y are Earth radii;*/
/*  scale is true at the reference latitude*/

ret  ::mapproj::toCylindricalEqualArea (type phi_1 , type lambda_0 , type phi_0 , type lambda , type phi) {
    variable degree
    set cos_phi_s [expr {cos($phi_1 * $degree)}]
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set x [expr {$lambda * $degree * $cos_phi_s}]
    set y0 [expr {sin($phi_0 * $degree) / $cos_phi_s}]
    set y [expr {sin($phi * $degree) / $cos_phi_s}]
    return [list $x [expr {$y - $y0}]]
}

/*  ::mapproj::fromCylindricalEqualArea --*/
/* */
/*  Solve a cylindrical equal area map projection for the*/
/*  latitude and longitude represented by a point on the map.*/
/* */
/*  Parameters:*/
/*  phi_1 -- Standard latitude in degrees*/
/*  phi_0 -- Latitude of the center of the sheet in degrees*/
/*  lambda_0 -- Longitude of the center of sheet in degrees*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromCylindricalEqualArea (type phi_1 , type lambda_0 , type phi_0 , type x , type y) {
    variable degree
    variable radian
    set cos_phi_s [expr {cos($phi_1 * $degree)}]
    set lambda [expr {$lambda_0 + $x / $cos_phi_s * $radian + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set y0 [expr {sin($phi_0 * $degree) / $cos_phi_s}]
    set phi [expr {asin(($y + $y0) * $cos_phi_s) * $radian}]
    return [list $lambda $phi]
}

/*  ::mapproj::toMercator --*/
/* */
/*  Project a latitude and longitude into the Mercator projection*/
/*  co-ordinates.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  phi_0 -- Latitude of the center of sheet in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates in Earth radii. Scale is true*/
/*  at the Equator and increases without bounds toward the Poles.*/

ret  ::mapproj::toMercator (type lambda_0 , type phi_0 , type lambda , type phi) {
    variable trace; if {[info exists trace]} { puts [info level 0] }
    variable degree
    variable quarterpi
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set x [expr {$lambda * $degree}]
    set y [expr {log(tan($quarterpi + 0.5 * $phi * $degree))}]
    set y0 [expr {log(tan($quarterpi + 0.5 * $phi_0 * $degree))}]
    if {[info exists trace]} { puts "[info level 0] -> $x [expr {$y - $y0}]" }
    return [list $x [expr {$y - $y0}]]
}

/*  ::mapproj::fromMercator --*/
/* */
/*  Converts Mercator map co-ordinates to latitude and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  phi_0 -- Latitude of the center of sheet in degrees*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromMercator (type lambda_0 , type phi_0 , type x , type y) {
    variable quarterpi
    variable degree
    variable radian
    set y0 [expr {log(tan($quarterpi + 0.5 * $phi_0 * $degree))}]
    set lambda [expr {$lambda_0 + $x  * $radian + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$radian * 2.0 * atan(exp($y + $y0)) - 90.0}]
    return [list $lambda $phi]
}

/*  ::mapproj::toMillerCylindrical --*/
/* */
/*  Project a latitude and longitude into the Miller Cylindrical projection*/
/*  co-ordinates.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates.  The x co-ordinate ranges from*/
/*  -$pi to pi.*/

ret  ::mapproj::toMillerCylindrical (type lambda_0 , type lambda , type phi) {
    foreach {x y} [toMercator $lambda_0 0.0 \
               $lambda [expr {0.8 * $phi}]] break
    set y [expr {1.25 * $y}]
    return [list $x $y]
}

/*  ::mapproj::fromMillerCylindrical --*/
/* */
/*  Converts Miller Cylindrical projected  map co-ordinates */
/*  to latitude and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromMillerCylindrical (type lambda_0 , type x , type y) {
    foreach {lambda phi} [fromMercator $lambda_0 0.0 \
                  $x [expr {0.8 * $y}]] break
    return [list $lambda [expr {1.25 * $phi}]]
}

/*  ::mapproj::toSinusoidal --*/
/* */
/*  Project a latitude and longitude into the sinusoidal*/
/*  (Sanson-Flamsteed) projection.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  phi_0 -- Latitude of the center of the sheet, in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates, in Earth radii. */
/*  Scale is true along the Equator and central meridian.*/

ret  ::mapproj::toSinusoidal (type lambda_0 , type phi_0 , type lambda , type phi) {
    variable degree
    variable quarterpi
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.) * $degree}]
    set phi [expr {$phi * $degree}]
    set x [expr {$lambda * cos($phi)}]
    set phi [expr {$phi - $phi_0 * $degree}]
    return [list $x $phi]
}

/*  ::mapproj::fromSinusoidal --*/
/* */
/*  Converts sinusoidal (Sanson-Flamsteed) map co-ordinates */
/*  to latitude and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  phi_0 -- Latitude of the center of the sheet, in degrees*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromSinusoidal (type lambda_0 , type phi_0 , type x , type y) {
    variable degree
    variable radian
    set y [expr {$y + $phi_0 * $degree}]
    set phi [expr {$y * $radian}]
    set lambda [expr {180. + $lambda_0 + $radian * $x / cos($y)}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda $phi]
}

/*  ::mapproj::toMollweide --*/
/* */
/*  Project a latitude and longitude into the Mollweide projection*/
/*  co-ordinates.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates, in Earth radii. */
/*  Scale is true along the 40 deg 44 min parallels*/

ret  ::mapproj::toMollweide (type lambda_0 , type lambda , type phi) {
    variable degree
    variable pi
    variable halfpi
    variable sqrt2
    variable sqrt8
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.) * $degree}]
    set phi [expr {$phi * $degree}]
    set theta [expr {2.0 * asin(2.0 * $phi / $pi)}]
    set diff 1.0
    set pisinphi [expr {$pi * sin($phi)}]
    while {abs($diff) >= 1.0e-4} {
    set diff [expr {($theta + sin($theta) - $pisinphi)
            / (1.0 + cos($theta))}]
    set theta [expr {$theta - $diff}]
    }
    set theta [expr {0.5 * $theta}]
    set x [expr {$sqrt8 * $lambda * cos($theta) / $pi}]
    set y [expr {$sqrt2 * sin($theta)}]
    return [list $x $y]
}

/*  ::mapproj::fromMollweide --*/
/* */
/*  Converts Mollweide projected  map co-ordinates to latitude */
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromMollweide (type lambda_0 , type x , type y) {
    variable pi
    variable radian
    variable degree
    variable halfpi
    variable sqrt2
    variable sqrt8
    set theta [expr {asin($y / $sqrt2)}]
    set lambda [expr {$lambda_0 + $radian * $pi * $x / 
              ($sqrt8 * cos($theta)) + 180.}]
    set phi [expr {asin((2.0 * $theta + sin(2.0 * $theta)) / $pi)}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda [expr {$phi * $radian}]]
}

/*  ::mapproj::toEckertIV --*/
/* */
/*  Project a latitude and longitude into the Eckert IV projection*/
/*  co-ordinates.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates.  Scale is true along the 40 deg 30 min*/
/*  parallels.*/

ret  ::mapproj::toEckertIV (type lambda_0 , type lambda , type phi) {
    variable degree
    variable pi
    variable halfpi
    variable EckertIVK1
    variable EckertIVK2
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.) * $degree}]
    set phi [expr {$phi * $degree}]
    set theta [expr {$phi / 2}]
    set diff 1.0
    set A [expr {(2.0 + $halfpi) * sin($phi)}]
    while {abs($diff) >= 1.0e-4} {
    set costheta [expr {cos($theta)}]
    set sintheta [expr {sin($theta)}]
    set diff \
        [expr {($theta + $sintheta * $costheta + 2.0 * sin($theta) - $A)
           / (2.0 * $costheta * (1.0 + $costheta))}]
    set theta [expr {$theta - $diff}]
    }
    set x [expr {$EckertIVK1 * $lambda * (1.0 + cos($theta))}]
    set y [expr {$EckertIVK2 * sin($theta)}]
    return [list $x $y]
}

/*  ::mapproj::fromEckertIV --*/
/* */
/*  Converts Eckert IV projected  map co-ordinates to latitude */
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromEckertIV (type lambda_0 , type x , type y) {
    variable pi
    variable radian
    variable degree
    variable halfpi
    variable sqrt2
    variable EckertIVK1
    variable EckertIVK2
    set sintheta [expr {$y / $EckertIVK2}]
    set costheta [expr {sqrt(1.0 - $sintheta * $sintheta)}]
    set theta [expr {atan2($sintheta, $costheta)}]
    set phi [expr {asin(($theta + $sintheta*$costheta + 2.*$sintheta)
            / (2. + $halfpi))}]
    set lambda [expr {180.0 + $lambda_0
              + $radian / $EckertIVK1 * $x / (1.0 + $costheta)}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda [expr {$phi * $radian}]]
}

/*  ::mapproj::toEckertVI --*/
/* */
/*  Project a latitude and longitude into the Eckert IV projection*/
/*  co-ordinates.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates.  Scale is true along the 40 deg 30 min*/
/*  parallels.*/

ret  ::mapproj::toEckertVI (type lambda_0 , type lambda , type phi) {
    variable degree
    variable pi
    variable halfpi
    variable EckertVIK1
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.) * $degree}]
    set phi [expr {$phi * $degree}]
    set theta [expr {$phi / 2}]
    set diff 1.0
    set A [expr {(1.0 + $halfpi) * sin($phi)}]
    while {abs($diff) >= 1.0e-4} {
    set costheta [expr {cos($theta)}]
    set sintheta [expr {sin($theta)}]
    set diff \
        [expr {($theta + $sintheta - $A)
           / (1.0 + $costheta)}]
    set theta [expr {$theta - $diff}]
    }
    set x [expr {$lambda * (1.0 + cos($theta)) / $EckertVIK1}]
    set y [expr {2.0 * $theta / $EckertVIK1}]
    return [list $x $y]
}

/*  ::mapproj::fromEckertVI --*/
/* */
/*  Converts Eckert IV projected  map co-ordinates to latitude */
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromEckertVI (type lambda_0 , type x , type y) {
    variable pi
    variable radian
    variable degree
    variable halfpi
    variable sqrt2
    variable EckertVIK1
    puts [info level 0]
    set theta [expr {0.5 * $EckertVIK1 * $y}]
    puts [list theta = $theta]
    set phi [expr {asin(($theta + sin($theta)) / (1.0 + $halfpi))}]
    puts [list phi = $phi]
    set lambda [expr {180.0 + $lambda_0 + $radian * $EckertVIK1 * $x
              / (1 + cos($theta))}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda [expr {$phi * $radian}]]
}

/*  ::mapproj::toRobinson --*/
/* */
/*  Project a latitude and longitude into the Robinson projection.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates, in Earth radii. */
/*  Scale is true along the Equator.*/

ret  ::mapproj::toRobinson (type lambda_0 , type lambda , type phi) {
    variable RobinsonLatitude
    variable RobinsonSplinePLEN
    variable RobinsonSplinePDFE
    variable RobinsonM
    variable pi
    variable degree
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set y [math::interpolate::interp-cubic-splines $RobinsonSplinePDFE $phi]
    set y [expr {$RobinsonM * $y}]
    set s [math::interpolate::interp-cubic-splines $RobinsonSplinePLEN $phi]
    set x [expr {$degree * $s * $lambda}]
    return [list $x $y]
}

/*  ::mapproj::fromRobinson --*/
/* */
/*  Solve the Robinson projection for the*/
/*  latitude and longitude represented by a point on the map.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromRobinson (type lambda_0 , type x , type y) {
    variable RobinsonLatitude
    variable RobinsonPDFE
    variable RobinsonSplinePLEN
    variable RobinsonSplinePDFE
    variable RobinsonM
    variable radian

    # We know that Robinson latitudes are equally spaced from [-90..90]
    # at 5-degree intervals.  Find the values for RobinsonPDFE that
    # bracket the y co-ordinate.

    set y [expr {$y / $RobinsonM}]
    set l 0
    set u [expr {[llength $RobinsonPDFE] - 1}]
    while {$l < $u} {
    set m [expr {($l + $u + 1) / 2}]
    if {$y >= [lindex $RobinsonPDFE $m]} {
        set l $m
    } else {
        set u [expr {$m - 1}]
    }
    }
    set u [lindex $RobinsonLatitude [expr {$l+1}]]
    set l [lindex $RobinsonLatitude $l]
    for {set i 0} {$i < 12} {incr i} {
    set m [expr {0.5 * ($u + $l)}]
    set ystar [math::interpolate::interp-cubic-splines \
               $RobinsonSplinePDFE $m]
    if {$ystar < $y} {
        set l $m
    } else {
        set u $m
    }
    }
    puts "latitude $m"
    set s [math::interpolate::interp-cubic-splines $RobinsonSplinePLEN $m]
    puts "parallel length $s"
    set lambda [expr {180.0 + $lambda_0 + $radian * $x / $s}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda $m]
}

/*  ::mapproj::toCassini --*/
/* */
/*  Project a latitude and longitude into the Cassini projection.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection in degrees*/
/*  phi_0 -- Latitude of the center of the sheet*/
/*  lambda --   Longitude of the point to be projected in degrees*/
/*  phi -- Latitude of the point to be projected in degrees*/
/* */
/*  Results:*/
/*  Returns x and y co-ordinates, in Earth radii. */
/*  Scale is true along the central meridian.*/

ret  ::mapproj::toCassini (type lambda_0 , type phi_0 , type lambda , type phi) {
    variable degree
    variable pi
    variable twopi
    variable quarterpi
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.) * $degree}]
    set phi [expr {$phi * $degree}]
    set x [expr {asin(cos($phi) * sin($lambda))}]
    set y [expr {atan2(tan($phi), cos($lambda)) - $degree * $phi_0}]
    if {$y < -$pi} {
    set y [expr {$y + $twopi}]
    } elseif {$y > $pi} {
    set y [expr {$y - $twopi}]
    }
    return [list $x $y]
}

/*  ::mapproj::fromCassini --*/
/* */
/*  Converts Cassini map co-ordinates to latitude and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  phi_0 -- Latitude of the center of the sheet*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromCassini (type lambda_0 , type phi_0 , type x , type y) {
    variable degree
    variable radian
    set y [expr {$y + $degree * $phi_0}]
    set phi [expr {$radian * asin(cos($x) * sin($y))}]
    set lambda [expr {180. + $lambda_0 + $radian * atan2(tan($x), cos($y))}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda $phi]
}

/*  ::mapproj::toPeirceQuincuncial*/
/* */
/*  Converts geodetic co-ordinates to the Peirce Quincuncial projection.*/
/* */
/*  Parameters:*/
/*  lambda_0 - Longitude of the central meridian.  (Conventionally, 20.0).*/
/*  lambda - Longitude of the point to be projected in degrees*/
/*  phi - Latitude of the point to be projected in degrees.*/
/* */
/*  Results:*/
/*  Returns a list of the x and y co-ordinates.*/

ret  ::mapproj::toPeirceQuincuncial (type lambda_0 , type lambda , type phi) {
    variable degree
    variable halfSqrt2
    variable pi
    variable quarterpi
    variable mquarterpi
    variable threequarterpi
    variable mthreequarterpi
    variable PeirceQuincuncialScale

    # Convert latitude and longitude to radians relative to the
    # central meridian

    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.) * $degree}]
    set phi [expr {$phi * $degree}]

    # Compute the auxiliary quantities 'm' and 'n'. Set 'm' to match
    # the sign of 'lambda' and 'n' to be positive if |lambda| > pi/2

    set cos_phiosqrt2 [expr {$halfSqrt2 * cos($phi)}]
    set cos_lambda [expr {cos($lambda)}]
    set sin_lambda [expr {sin($lambda)}]
    set cos_a [expr {$cos_phiosqrt2 * ($sin_lambda + $cos_lambda)}]
    set cos_b [expr {$cos_phiosqrt2 * ($sin_lambda - $cos_lambda)}]
    set sin_a [expr {sqrt(1.0 - $cos_a * $cos_a)}]
    set sin_b [expr {sqrt(1.0 - $cos_b * $cos_b)}]
    set cos_a_cos_b [expr {$cos_a * $cos_b}]
    set sin_a_sin_b [expr {$sin_a * $sin_b}]
    set sin2_m [expr {1.0 + $cos_a_cos_b - $sin_a_sin_b}]
    set sin2_n [expr {1.0 - $cos_a_cos_b - $sin_a_sin_b}]
    if {$sin2_m < 0.0} {set sin2_m 0.0}
    set sin_m [expr {sqrt($sin2_m)}]
    if {$sin2_m > 1.0} { set sin2_m 1.0 }
    set cos_m [expr {sqrt(1.0 - $sin2_m)}]
    if {$sin_lambda < 0.0} {
    set sin_m [expr {-$sin_m}]
    }
   if {$sin2_n < 0.0} { set sin2_n 0.0 }
    set sin_n [expr {sqrt($sin2_n)}]
    if {$sin2_n > 1.0} { set sin2_n 1.0 }
    set cos_n [expr {sqrt(1.0 - $sin2_n)}]
    if {$cos_lambda > 0.0} {
    set sin_n [expr {-$sin_n}]
    }

    # Compute elliptic integrals to map the disc to the square

    set x [ellFaux $cos_m $sin_m $halfSqrt2]
    set y [ellFaux $cos_n $sin_n $halfSqrt2]

    # Reflect the Southern Hemisphere outward

    if {$phi < 0} {
    if {$lambda < $mthreequarterpi} {
        set y [expr {$PeirceQuincuncialScale - $y}]
    } elseif {$lambda < $mquarterpi} {
        set x [expr {-$PeirceQuincuncialScale - $x}]
    } elseif {$lambda < $quarterpi} {
        set y [expr {-$PeirceQuincuncialScale - $y}]
    } elseif {$lambda < $threequarterpi} {
        set x [expr {$PeirceQuincuncialScale - $x}]
    } else {
        set y [expr {$PeirceQuincuncialScale - $y}]
    }
    }

    # Rotate the square by 45 degrees to fit the screen better

    set X [expr {($x - $y) * $halfSqrt2}]
    set Y [expr {($x + $y) * $halfSqrt2}]

    return [list $X $Y]
}

/*  ::mapproj::fromPeirceQuincuncial --*/
/* */
/*  Converts Peirce Quincuncial map co-ordinates to latitude and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0 -- Longitude of the center of projection*/
/*  x,y -- normalized x and y co-ordinates of a point on the map*/
/* */
/*  Results:*/
/*  Returns a list consisting of the longitude and latitude in degrees.*/

ret  ::mapproj::fromPeirceQuincuncial (type lambda_0 , type x , type y) {
    variable halfSqrt2
    variable radian
    variable pi
    variable halfpi
    variable quarterpi
    variable PeirceQuincuncialScale
    variable PeirceQuincuncialLimit

    # Rotate x and y 45 degrees

    set X [expr {($x + $y) * $halfSqrt2}]
    set Y [expr {($y - $x) * $halfSqrt2}]

    # Reflect Southern Hemisphere into the Northern

    set southern 0
    if {$X < -$PeirceQuincuncialLimit} {
    set X [expr {-$PeirceQuincuncialScale - $X}]
    set southern 1
    } elseif {$X > $PeirceQuincuncialLimit} {
    set X [expr {$PeirceQuincuncialScale - $X}]
    set southern 1
    } elseif {$Y < -$PeirceQuincuncialLimit} {
    set Y [expr {-$PeirceQuincuncialScale - $Y}]
    set southern 1
    } elseif {$Y > $PeirceQuincuncialLimit} {
    set Y [expr {$PeirceQuincuncialScale - $Y}]
    set southern 1
    }

    # Now we know that latitude will be positive.  If X is negative, then
    # longitude will be negative; reflect the Western Hemisphere into the
    # Eastern.

    set western 0
    if {$X < 0.0} {
    set western 1
    set X [expr {-$X}]
    }

    # If Y is positive, the point is in the back hemisphere.  Reflect
    # it to the front.

    set back 0
    if {$Y > 0.0} {
    set back 1
    set Y [expr {-$Y}]
    }

    # Finally, constrain longitude to be less than pi/4, by reflecting across
    # the 45 degree meridian.

    set complement 0
    if {$X > -$Y} {
    set complement 1
    set t [expr {-$X}]
    set X [expr {-$Y}]
    set Y $t
    }

    # Compute the elliptic functions to map the plane onto the sphere

    set cnx [cn $X $halfSqrt2]
    set cny [cn $Y $halfSqrt2]

    # Undo the mapping to latitude and longitude

    set a1 [expr {acos(-$cnx * $cnx)}]
    set a2 [expr {acos($cny * $cny)}]
    set b [expr {0.5 * ($a1 + $a2)}]
    set a [expr {0.5 * ($a1 - $a2)}]
    set cos_a [expr {cos($a)}]
    set cos_b [expr {-cos($b)}]
    set lambda [expr {$quarterpi - atan2($cos_b, $cos_a)}]
    set phi [expr {acos(hypot($cos_b, $cos_a))}]

    # Undo the reflections that were done above, to get correct latitude
    # and longitude

    if {$complement} {
    set lambda [expr {$halfpi - $lambda}]
    }
    if {$back} {
    set lambda [expr {$pi - $lambda}]
    }
    if {$western} {
    set lambda [expr {-$lambda}]
    }
    if {$southern} {
    set phi [expr {-$phi}]
    }

    # Convert latitude and longitude to degrees

    set lambda [expr {$lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::RotateCartesianY*/
/* */
/*  Rotates Cartesian co-ordinates about the y axis*/
/* */
/*  Parameters:*/
/*  phi - Angle (in degrees) about which to rotate.*/
/*  x,y,z - Cartesian co-ordinates*/
/* */
/*  Results:*/
/*  Returns a three-element list giving the rotated co-ordinates.*/

ret  ::mapproj::RotateCartesianY (type phi , type x , type y , type z) {
    variable degree
    set phi [expr {$degree * $phi}]
    set cos_phi [expr {cos($phi)}]
    set sin_phi [expr {sin($phi)}]
    return [list [expr {$x * $cos_phi - $z * $sin_phi}] \
        $y \
        [expr {$z * $cos_phi + $x * $sin_phi}]]
}

/*  ::mapproj::ToCartesian --*/
/* */
/*  Converts geodetic co-ordinates to Cartesian*/
/* */
/*  Parameters:*/
/*  lambda - Longitude of the point to be projected, in degrees*/
/*  phi - Latitude of the point to be projected, in degrees*/
/* */
/*  Results:*/
/*  Returns a three-element list, x, y, z where x is the component*/
/*  in the direction of longitude 0, latitude 0, y is the component*/
/*  in the direction of longitude 90 East, latitude 0, and*/
/*  z is the component in the direction of the North Pole*/
/* */
/*  Auxiliary procedure used in several projections to convert*/
/*  geodetic coordinates to Cartesian range and bearing.*/

ret  ::mapproj::ToCartesian (type lambda , type phi) {
    variable degree
    set lambda [expr {$degree * $lambda}]
    set phi [expr {$degree * $phi}]
    set cos_phi [expr cos($phi)]
    return [list [expr {$cos_phi * cos($lambda)}] \
        [expr {$cos_phi * sin($lambda)}] \
        [expr {sin($phi)}]]
}

/*  ::mapproj::CartesianToRangeAndBearing*/
/* */
/*  Transforms view-relative Cartesian co-ordinates to range and*/
/*  bearing.*/
/* */
/*  Parameters:*/
/*  x,y,z - Cartesian co-ordinates relative to center of Earth;*/
/*  +x points to the viewer and +z to the "view-up" direction.*/
/* */
/*  Results:*/
/*  Returns a three-element list containing, in order,*/
/*  the cosine (easting) of the bearing, the sine (northing) of the*/
/*  bearing, and the range.*/

ret  ::mapproj::CartesianToRangeAndBearing (type x , type y , type z) {
    set c [expr {hypot($z, $y)}]
    if {$c == 0} {
    set cos_b 1.0
    set sin_b 0.0
    } else {
    set cos_b [expr {$y / $c}]
    set sin_b [expr {$z / $c}]
    }
    set range [expr {atan2($c, $x)}]
    return [list $cos_b $sin_b $range]
}

/*  ::mapproj::RangeAndBearingToCartesian --*/
/* */
/*  Converts range and bearing to Cartesian co-ordinates.*/
/* */
/*  Parameters:*/
/*  cos_b, sin_b -- Cosine (easting) and sine (northing) of the bearing*/
/*  range - Range, in Earth radii*/
/* */
/*  Results:*/
/*  Returns Cartesian co-ordinates relative to center of Earth.*/
/*  x is toward the station, and z is "view up"*/

ret  ::mapproj::RangeAndBearingToCartesian (type cos_, type b , type sin_, type b , type range) {
    set c [expr {sin($range)}]
    set x [expr {cos($range)}]
    set y [expr {$cos_b * $c}]
    set z [expr {$sin_b * $c}]
    return [list $x $y $z]
}
    

/*  ::mapproj::CartesianToSpherical --*/
/* */
/*  Transforms Cartesian x, y, z to spherical co-ordinates*/
/* */
/*  Parameters:*/
/*  x, y, z -- Coordinates of a point on the surface of the Earth,*/
/*             in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two-element list comprising longitude and latitude*/
/*  in radians*/

ret  ::mapproj::CartesianToSpherical (type x , type y , type z) {
    return [list [expr {atan2($y, $x)}] [expr {atan2($z, hypot($y, $x))}]]
}

/*  ::mapproj::toOrthographic --*/
/* */
/*  Transforms latitude and longitude to x and y co-ordinates*/
/*  on an orthographic projection.  Scale is true only at the*/
/*  point of projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/*             in degrees*/
/* */
/*  Results:*/
/*  Returns map x and y co-ordinates, in Earth radii.*/

ret  ::mapproj::toOrthographic (type lambda_0 , type phi_0 , type lambda , type phi) {
    foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \
    break
    foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \
    break
    if {$x < 0} {
    return {}
    } else {
    return [list $y $z]
    }
}

/*  ::mapproj::fromOrthographic --*/
/* */
/*  Transforms x and y on an orthographic projection to latitude*/
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  x, y -- Co-ordinates of the projected point, in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two element list containing longitude and latitude*/
/*  in degrees.*/

ret  ::mapproj::fromOrthographic (type lambda_0 , type phi_0 , type x , type y) {
    variable radian
    set r [expr {hypot($x, $y)}]
    set alpha [expr {asin($r)}]
    set z [expr {sqrt(1.0 - $r*$r)}]
    foreach {x y z} [RotateCartesianY $phi_0 $z $x $y] break
    foreach {lambda phi} [CartesianToSpherical $x $y $z] break

    # Convert latitude and longitude to degrees

    set lambda [expr {$lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::toStereographic --*/
/* */
/*  Transforms latitude and longitude to x and y co-ordinates*/
/*  on an orthographic projection.  Scale is true only at the*/
/*  point of projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/*             in degrees*/
/* */
/*  Results:*/
/*  Returns map x and y co-ordinates, in Earth radii.*/

ret  ::mapproj::toStereographic (type lambda_0 , type phi_0 , type lambda , type phi) {
    foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \
    break
    foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \
    break
    if {$x < -0.5} {
    return {}
    } else {
    set y [expr {2. * $y / (1. + $x)}]
    set z [expr {2. * $z / (1. + $x)}]
    return [list $y $z]
    }
}

/*  ::mapproj::fromStereographic --*/
/* */
/*  Transforms x and y on an orthographic projection to latitude*/
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  x, y -- Co-ordinates of the projected point, in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two element list containing longitude and latitude*/
/*  in degrees.*/

ret  ::mapproj::fromStereographic (type lambda_0 , type phi_0 , type x , type y) {
    variable radian
    variable halfpi
    set denom [expr {4.0 + $x*$x + $y*$y}]
    foreach {x y z} [list \
             [expr {(4.0 - $x*$x - $y*$y) / $denom}] \
             [expr {4. * $x / $denom}] \
             [expr {4. * $y / $denom}]] break
    
    foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break
    foreach {lambda phi} [CartesianToSpherical $x $y $z] break

    # Convert latitude and longitude to degrees

    set lambda [expr {$lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::toGnomonic --*/
/* */
/*  Transforms latitude and longitude to x and y co-ordinates*/
/*  on an orthographic projection.  Scale is true only at the*/
/*  point of projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/*             in degrees*/
/* */
/*  Results:*/
/*  Returns map x and y co-ordinates, in Earth radii.*/

ret  ::mapproj::toGnomonic (type lambda_0 , type phi_0 , type lambda , type phi) {
    foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \
    break
    foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \
    break
    if {$x < 0.01} {
    return {}
    } else {
    set y [expr {$y / $x}]
    set z [expr {$z / $x}]
    return [list $y $z]
    }
}

/*  ::mapproj::fromGnomonic --*/
/* */
/*  Transforms x and y on an orthographic projection to latitude*/
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  x, y -- Co-ordinates of the projected point, in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two element list containing longitude and latitude*/
/*  in degrees.*/

ret  ::mapproj::fromGnomonic (type lambda_0 , type phi_0 , type x , type y) {
    variable radian
    variable halfpi
    set denom [expr {hypot(1.0, hypot($x, $y))}]
    foreach {x y z} [list \
             [expr {1.0 / $denom}] \
             [expr {$x / $denom}] \
             [expr {$y / $denom}]] break
    
    foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break
    foreach {lambda phi} [CartesianToSpherical $x $y $z] break

    # Convert latitude and longitude to degrees

    set lambda [expr {$lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::toAzimuthalEquidistant --*/
/* */
/*  Transforms latitude and longitude to x and y co-ordinates*/
/*  on an orthographic projection.  Scale is true only at the*/
/*  point of projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/*             in degrees*/
/* */
/*  Results:*/
/*  Returns map x and y co-ordinates, in Earth radii.*/

ret  ::mapproj::toAzimuthalEquidistant (type lambda_0 , type phi_0 , type lambda , type phi) {
    foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \
    break
    foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \
    break
    foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break
    return [list [expr {$cs * $range}] [expr {$sn * $range}]]
}

/*  ::mapproj::fromAzimuthalEquidistant --*/
/* */
/*  Transforms x and y on an orthographic projection to latitude*/
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  x, y -- Co-ordinates of the projected point, in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two element list containing longitude and latitude*/
/*  in degrees.*/

ret  ::mapproj::fromAzimuthalEquidistant (type lambda_0 , type phi_0 , type x , type y) {
    variable radian
    variable halfpi

    set range [expr {hypot($y, $x)}]
    set cos_b [expr {$x / $range}]
    set sin_b [expr {$y / $range}]
    foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break
    foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break
    foreach {lambda phi} [CartesianToSpherical $x $y $z] break

    # Convert latitude and longitude to degrees

    set lambda [expr {$lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::toLambertAzimuthalEqualArea --*/
/* */
/*  Transforms latitude and longitude to x and y co-ordinates*/
/*  on an orthographic projection.  Scale is true only at the*/
/*  point of projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/*             in degrees*/
/* */
/*  Results:*/
/*  Returns map x and y co-ordinates, in Earth radii.*/

ret  ::mapproj::toLambertAzimuthalEqualArea (type lambda_0 , type phi_0 , type lambda , type phi) {
    foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \
    break
    foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \
    break
    foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break
    set range [expr {2.0 * sin(0.5 * $range)}]
    return [list [expr {$cs * $range}] [expr {$sn * $range}]]
}

/*  ::mapproj::fromLambertAzimuthalEqualArea --*/
/* */
/*  Transforms x and y on an orthographic projection to latitude*/
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  x, y -- Co-ordinates of the projected point, in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two element list containing longitude and latitude*/
/*  in degrees.*/

ret  ::mapproj::fromLambertAzimuthalEqualArea (type lambda_0 , type phi_0 , type x , type y) {
    variable radian
    variable halfpi

    set range [expr {hypot($y, $x)}]
    set cos_b [expr {$x / $range}]
    set sin_b [expr {$y / $range}]
    set range [expr {2.0 * asin(0.5 * $range)}]
    foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break
    foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break
    foreach {lambda phi} [CartesianToSpherical $x $y $z] break

    # Convert latitude and longitude to degrees

    set lambda [expr {$lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::toHammer --*/
/* */
/*  Transforms latitude and longitude to x and y co-ordinates*/
/*  on an orthographic projection.  Scale is true only at the*/
/*  point of projection.*/
/* */
/*  Parameters:*/
/*  lambda_0-- Longitude of the center of projection*/
/*             in degrees*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/*             in degrees*/
/* */
/*  Results:*/
/*  Returns map x and y co-ordinates, in Earth radii.*/

ret  ::mapproj::toHammer (type lambda_0 , type lambda , type phi) {
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.0}]
    foreach {x y z} [ToCartesian [expr {$lambda/2.}] $phi] \
    break
    foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break
    set range [expr {2.0 * sin(0.5 * $range)}]
    return [list [expr {2.0 * $cs * $range}] [expr {$sn * $range}]]
}

/*  ::mapproj::fromHammer --*/
/* */
/*  Transforms x and y on an orthographic projection to latitude*/
/*  and longitude.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of projection*/
/*             in degrees*/
/*  x, y -- Co-ordinates of the projected point, in Earth radii*/
/* */
/*  Results:*/
/*  Returns a two element list containing longitude and latitude*/
/*  in degrees.*/

ret  ::mapproj::fromHammer (type lambda_0 , type x , type y) {
    variable radian
    variable halfpi

    set x [expr {0.5 * $x}]
    set range [expr {hypot($y, $x)}]
    set cos_b [expr {$x / $range}]
    set sin_b [expr {$y / $range}]
    set range [expr {2.0 * asin(0.5 * $range)}]
    foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break
    foreach {lambda phi} [CartesianToSpherical $x $y $z] break

    # Convert latitude and longitude to degrees

    set lambda [expr {2.0 * $lambda * $radian + 180. + $lambda_0}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    set phi [expr {$phi * $radian}]

    return [list $lambda $phi]
}

/*  ::mapproj::toConicEquidistant*/
/* */
/*  Converts latitude and longitude to map co-ordinates on a*/
/*  conic equidistant projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of the sheet.*/
/*  phi_1, phi_2 -- Latitudes of the two standard parallels at which scale*/
/*          is true*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/* */
/*  Results:*/
/*  Returns a list of map x and y measured in Earth radii.*/

ret  ::mapproj::toConicEquidistant (type lambda_0 , type phi_0 , type phi_1 , type phi_2 , type lambda , type phi) {
    variable degree
    set phi_0 [expr {$phi_0 * $degree}]
    set phi_1 [expr {$phi_1 * $degree}]
    set phi_2 [expr {$phi_2 * $degree}]
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.0) * $degree}]
    set phi [expr {$phi * $degree}]
    set cos_phi_1 [expr {cos($phi_1)}]
    set n [expr {($cos_phi_1 - cos($phi_2)) / ($phi_2 - $phi_1)}]
    set G [expr {$cos_phi_1 / $n + $phi_1}]
    set rho_0 [expr {$G - $phi_0}]
    set theta [expr {$n * $lambda}]
    set rho [expr {$G - $phi}]
    set x [expr {$rho * sin($theta)}]
    set y [expr {$rho_0 - $rho * cos($theta)}]
    return [list $x $y]
}

/*  ::mapproj::fromConicEquidistant --*/
/* */
/*  Unprojects map x and y in a conic equidistant projection to*/
/*  latitude and longitude*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of the sheet.*/
/*  phi_1, phi_2 -- Latitudes of the two standard parallels at which scale*/
/*          is true*/
/*  x, y -- Map co-ordinates in Earth radii*/
/* */
/*  Results:*/
/*  Returns a list of longitude and latitude in degrees.*/

ret  ::mapproj::fromConicEquidistant (type lambda_0 , type phi_0 , type phi_1 , type phi_2 , type x , type y) {
    variable degree
    variable radian
    set phi_0 [expr {$phi_0 * $degree}]
    set phi_1 [expr {$phi_1 * $degree}]
    set phi_2 [expr {$phi_2 * $degree}]
    set cos_phi_1 [expr {cos($phi_1)}]
    set n [expr {($cos_phi_1 - cos($phi_2)) / ($phi_2 - $phi_1)}]
    set G [expr {$cos_phi_1 / $n + $phi_1}]
    set rho_0 [expr {$G - $phi_0}]
    set rho_0my [expr {$rho_0 - $y}]
    set theta [expr {atan2($x, $rho_0my)}]
    set rho [expr {sqrt($x*$x + $rho_0my * $rho_0my)}]
    if {$n < 0.0} {set rho [expr {-$rho}]}
    set phi [expr {($G - $rho) * $radian}]
    set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda $phi]
}

/*  ::mapproj::toAlbersEqualAreaConic*/
/* */
/*  Converts latitude and longitude to map co-ordinates on a*/
/*  conic equal-area projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of the sheet.*/
/*  phi_1, phi_2 -- Latitudes of the two standard parallels at which scale*/
/*          is true*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/* */
/*  Results:*/
/*  Returns a list of map x and y measured in Earth radii.*/

ret  ::mapproj::toAlbersEqualAreaConic (type lambda_0 , type phi_0 , type phi_1 , type phi_2 , type lambda , type phi) {
    variable degree
    set phi_0 [expr {$phi_0 * $degree}]
    set phi_1 [expr {$phi_1 * $degree}]
    set phi_2 [expr {$phi_2 * $degree}]
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.0) * $degree}]
    set phi [expr {$phi * $degree}]
    set cos_phi_1 [expr {cos($phi_1)}]
    set sin_phi_1 [expr {sin($phi_1)}]
    set n [expr {0.5 * ($sin_phi_1 + sin($phi_2))}]
    set theta [expr {$n * $lambda}]
    set C [expr {$cos_phi_1 * $cos_phi_1 + 2.0 * $n * $sin_phi_1}]
    set rho [expr {sqrt($C - 2.0 * $n * sin($phi)) / $n}]
    set rho_0 [expr {sqrt($C - 2.0 * $n * sin($phi_0)) / $n}]
    set x [expr {$rho * sin($theta)}]
    set y [expr {$rho_0 - $rho * cos($theta)}]
    return [list $x $y]
}

/*  ::mapproj::fromAlbersEqualAreaConic --*/
/* */
/*  Unprojects map x and y in a conic equal-area projection to*/
/*  latitude and longitude*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of the sheet.*/
/*  phi_1, phi_2 -- Latitudes of the two standard parallels at which scale*/
/*          is true*/
/*  x, y -- Map co-ordinates in Earth radii*/
/* */
/*  Results:*/
/*  Returns a list of longitude and latitude in degrees.*/

ret  ::mapproj::fromAlbersEqualAreaConic (type lambda_0 , type phi_0 , type phi_1 , type phi_2 , type x , type y) {
    variable degree
    variable radian
    variable twopi
    set phi_0 [expr {$phi_0 * $degree}]
    set phi_1 [expr {$phi_1 * $degree}]
    set phi_2 [expr {$phi_2 * $degree}]
    set cos_phi_1 [expr {cos($phi_1)}]
    set sin_phi_1 [expr {sin($phi_1)}]
    set n [expr {0.5 * ($sin_phi_1 + sin($phi_2))}]
    set C [expr {$cos_phi_1 * $cos_phi_1 + 2.0 * $n * $sin_phi_1}]
    set rho_0 [expr {sqrt($C - 2.0 * $n * sin($phi_0)) / $n}]
    set theta [expr {atan2($x, $rho_0 - $y)}]
    set rho [expr {hypot($x, $rho_0 - $y)}]
    set phi [expr {$radian * asin(($C - $rho*$rho*$n*$n) / (2.0 * $n))}]
    set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda $phi]
}

/*  ::mapproj::toLambertConformalConic*/
/* */
/*  Converts latitude and longitude to map co-ordinates on a*/
/*  conformal conic projection.*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of the sheet.*/
/*  phi_1, phi_2 -- Latitudes of the two standard parallels at which scale*/
/*          is true*/
/*  lambda, phi -- Longitude and latitude of the point to be projected*/
/* */
/*  Results:*/
/*  Returns a list of map x and y measured in Earth radii.*/

ret  ::mapproj::toLambertConformalConic (type lambda_0 , type phi_0 , type phi_1 , type phi_2 , type lambda , type phi) {
    variable degree
    variable quarterpi
    set phi_0 [expr {$phi_0 * $degree}]
    set phi_1 [expr {$phi_1 * $degree}]
    set phi_2 [expr {$phi_2 * $degree}]
    set lambda [expr {$lambda - $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {($lambda - 180.0) * $degree}]
    set phi [expr {$phi * $degree}]
    set cos_phi_1 [expr {cos($phi_1)}]
    set sin_phi_1 [expr {sin($phi_1)}]
    set tan1 [expr {tan($quarterpi + 0.5 * $phi_1)}]
    set n [expr {log($cos_phi_1 / cos($phi_2))
         / log(tan($quarterpi + 0.5 * $phi_2) / $tan1)}]
    set F [expr {$cos_phi_1 * pow($tan1, $n) / $n}]
    set rho [expr {$F * pow(tan($quarterpi + 0.5 * $phi), -$n)}]
    set rho_0 [expr {$F * pow(tan($quarterpi + 0.5 * $phi_0), -$n)}]
    set x [expr {$rho * sin($n * $lambda)}]
    set y [expr {$rho_0 - $rho * cos($n * $lambda)}]
    return [list $x $y]
}

/*  ::mapproj::fromLambertConformalConic --*/
/* */
/*  Unprojects map x and y in a conformal conic projection to*/
/*  latitude and longitude*/
/* */
/*  Parameters:*/
/*  lambda_0, phi_0 -- Longitude and latitude of the center of the sheet.*/
/*  phi_1, phi_2 -- Latitudes of the two standard parallels at which scale*/
/*          is true*/
/*  x, y -- Map co-ordinates in Earth radii*/
/* */
/*  Results:*/
/*  Returns a list of longitude and latitude in degrees.*/

ret  ::mapproj::fromLambertConformalConic (type lambda_0 , type phi_0 , type phi_1 , type phi_2 , type x , type y) {
    variable degree
    variable radian
    variable quarterpi
    set phi_0 [expr {$phi_0 * $degree}]
    set phi_1 [expr {$phi_1 * $degree}]
    set phi_2 [expr {$phi_2 * $degree}]
    set cos_phi_1 [expr {cos($phi_1)}]
    set sin_phi_1 [expr {sin($phi_1)}]
    set tan1 [expr {tan($quarterpi + 0.5 * $phi_1)}]
    set n [expr {log($cos_phi_1 / cos($phi_2))
         / log(tan($quarterpi + 0.5 * $phi_2) / $tan1)}]
    set F [expr {$cos_phi_1 * pow($tan1, $n) / $n}]
    set rho_0 [expr {$F * pow(tan($quarterpi + 0.5 * $phi_0), -$n)}]
    set y [expr {$rho_0 - $y}]
    set rho [expr {sqrt($x*$x + $y*$y)}]
    if {$n < 0} { set rho [expr {-$rho}] }
    set theta [expr {atan2($x, $y)}]
    set phi [expr {$radian * 2 * atan(pow($F / $rho, 1.0 / $n)) - 90.}]
    set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}]
    if {$lambda < 0.0 || $lambda > 360.0} {
    set lambda [expr {$lambda - 360. * floor($lambda / 360.)}]
    }
    set lambda [expr {$lambda - 180.}]
    return [list $lambda $phi]
}

/*  Define commonly used cylindrical equal-area projections*/

ret  ::mapproj::toLambertCylindricalEqualArea (type lambda_0 , type phi_0 , type lambda , type phi) {
    toCylindricalEqualArea 0.0 $lambda_0 $phi_0 $lambda $phi
}
ret  ::mapproj::fromLambertCylindricalEqualArea (type lambda_0 , type phi_0 , type x , type y) {
    fromCylindricalEqualArea 0.0 $lambda_0 $phi_0 $x $y
}
ret  ::mapproj::toBehrmann (type lambda_0 , type phi_0 , type lambda , type phi) {
    toCylindricalEqualArea 30.0 $lambda_0 $phi_0 $lambda $phi
}
ret  ::mapproj::fromBehrmann (type lambda_0 , type phi_0 , type x , type y) {
    fromCylindricalEqualArea 30.0 $lambda_0 $phi_0 $x $y
}
ret  ::mapproj::toTrystanEdwards (type lambda_0 , type phi_0 , type lambda , type phi) {
    toCylindricalEqualArea 37.4 $lambda_0 $phi_0 $lambda $phi
}
ret  ::mapproj::fromTrystanEdwards (type lambda_0 , type phi_0 , type x , type y) {
    fromCylindricalEqualArea 37.4 $lambda_0 $phi_0 $x $y
}
ret  ::mapproj::toHoboDyer (type lambda_0 , type phi_0 , type lambda , type phi) {
    toCylindricalEqualArea 37.5 $lambda_0 $phi_0 $lambda $phi
}
ret  ::mapproj::fromHoboDyer (type lambda_0 , type phi_0 , type x , type y) {
    fromCylindricalEqualArea 37.5 $lambda_0 $phi_0 $x $y
}
ret  ::mapproj::toGallPeters (type lambda_0 , type phi_0 , type lambda , type phi) {
    toCylindricalEqualArea 45.0 $lambda_0 $phi_0 $lambda $phi
}
ret  ::mapproj::fromGallPeters (type lambda_0 , type phi_0 , type x , type y) {
    fromCylindricalEqualArea 45.0 $lambda_0 $phi_0 $x $y
}
ret  ::mapproj::toBalthasart (type lambda_0 , type phi_0 , type lambda , type phi) {
    toCylindricalEqualArea 50.0 $lambda_0 $phi_0 $lambda $phi
}
ret  ::mapproj::fromBalthasart (type lambda_0 , type phi_0 , type x , type y) {
    fromCylindricalEqualArea 50.0 $lambda_0 $phi_0 $x $y
}