#
#  gpsman --- GPS Manager: a manager for GPS receiver data
#
# Copyright (c) 1998-2009 Miguel Filgueiras mig@ncc.up.pt Universidade do Porto
#
#    This program is free software; you can redistribute it and/or modify
#      it under the terms of the GNU General Public License as published by
#      the Free Software Foundation; either version 2 of the License, or
#      (at your option) any later version.
#
#      This program is distributed in the hope that it will be useful,
#      but WITHOUT ANY WARRANTY; without even the implied warranty of
#      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#      GNU General Public License for more details.
#
#      You should have received a copy of the GNU General Public License
#      along with this program.
#
#  File: projs_main.tcl
#  Last change:  17 June 2009
#

# Includes contributions by
#  - Sandor Laky (laky.sandor_AT_freemail.hu) marked "SL contribution"
#  - Miguel Filgueiras (to code by others) marked "MF contribution/change"

### most of this code used to be at the end of file projections.tcl

### general initialization procedure for all main projections
# see also the comments before the procedures for UTM in projections.tcl

proc ProjInit {proj data datum ps} {
    # changes here must be reflected on proc BadProjSetup below
    # initialize projection parameters
    #  $ps is a list of latd,longd,datum to be projected
    # should call proc ProjParams if its parameters may be changed by
    #  the user and if $ASKPROJPARAMS is true
    # should call proc to compute auxiliary parameters if any
    # must return projection of first position in $ps
    global $data TMPARAM LNTFPARAMS ASKPROJPARAMS

    foreach "latd longd posdatum" [lindex $ps 0] { break }
    if { $posdatum != $datum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $posdatum $datum] { break }
    }
    set ${data}(datum) $datum
    set askauxcomp 1
    switch $proj {
	UTM {
	    # only the first position is used
	    # no parameters can be changed by the user
	    set askauxcomp 0
	    set p [DegreesToUTM $latd $longd $datum]
	    scan [lindex $p 0] %0d ze
	    set zn [lindex $p 1]
	    set ${data}(UTMzone) [list $ze $zn]
	    set ${data}(m_0) [expr -183+6.0*$ze]
	    set res [list [lindex $p 2] [lindex $p 3]]
	}
	TM {
	    # take averages of latitudes and of longitudes as central
	    #  coordinates, and use k0=0.9996
	    foreach "${data}(lat0) ${data}(long0)" \
		    [AverageLatLong $latd $longd $datum $ps] {}
	    set ${data}(k0) 0.9996
	    if { $ASKPROJPARAMS } {
		ProjParams change TM $data
	    }
	    set res [ConvToTM $latd $longd [set ${data}(lat0)] \
		          [set ${data}(long0)] [set ${data}(k0)] $datum]
	    set askauxcomp 0
	}
	BMN -  CTR -  GKK -  KKJP -  TWG {
	    # special cases of TM
	    # take averages of longitudes to find long0, and use lat0=0
	    set long0 [AverageLong $longd $datum $ps]
	    foreach "pr lmin lmax delta k0" $TMPARAM($proj) {}
	    if { $long0 < $lmin } {
		set long0 $lmin
	    } elseif { $long0 > $lmax } {
		set long0 $lmax
	    } else {
		set long0 [expr int(1.0*($long0-$lmin+0.5*$delta)/ \
			            $delta)*$delta+$lmin]
	    }
	    set ${data}(${pr}long0) $long0
	    set ${data}(m_0) $k0
	    if { $ASKPROJPARAMS } {
		ProjParams change $proj $data
	    }
	    set res [ConvToTM $latd $longd 0 [set ${data}(${pr}long0)] \
		              $k0 $datum]
	    set askauxcomp 0
	}
	LCC1 -  Stereogr -  SOM {
	    # Lambert Conformal Conic single parallel,
	    #  Stereographic and Swiss Oblique Mercator projections
	    # take averages of latitudes and of longitudes as central
	    #  coordinates, and use k0=1
	    # in Stereographic central latitude determines the case: polar,
	    #  oblique, equatorial
	    foreach "${data}(lat0) ${data}(long0)" \
		    [AverageLatLong $latd $longd $datum $ps] {}
	    set ${data}(k0) 1.0
	}
	LCC2 {
	    # Lambert Conformal Conic projection, two parallels
	    # take extremes of latitudes as latitudes of standard parallels,
	    #  averages of latitudes and longitudes as false origin coordinates
	    foreach "${data}(latF) ${data}(longF) \
		     ${data}(lat1) ${data}(lat2)" \
		    [AverageExtrLatLong $latd $longd $datum $ps] {}
	}
	LambNTF {
	    # as LCC2 but with zones with specific parameters
	    # only the first position is used
	    # no parameters can be changed by the user
	    set askauxcomp 0	    
	    set p [DegreesToLambNTF $latd $longd $datum]
	    if { [set z [lindex $p 0]] == "--" } {
		set z II
	    }
	    global LNTFz$z
	    array set $data [array get LNTFz$z]
	    set res [ProjLambNTFPoint $data $latd $longd $datum]
	}
	AEA {
	    # Albers Equal Area Conic projection
	    # take extremes of latitudes as latitudes of standard parallels,
	    #  average of longitudes as central longitude, 0 for central
	    #  latitude
	    foreach "x ${data}(long0) lamax lamin" \
		    [AverageExtrLatLong $latd $longd $datum $ps] {}
	    if { abs($lamax) < abs($lamin) } {
		set x $lamax ; set lamax $lamin ; set lamin $x
	    }
	    set ${data}(lat1) $lamax ; set ${data}(lat2) $lamin
	    set ${data}(lat0) 0
	}
	LEAC {
	    # Lambert Equal Area Conic projection
	    # take extreme of latitudes as latitude of standard parallel,
	    #  average of longitudes as central longitude, North polar
	    #  aspect if extreme latitude is nonnegative, 0 for central
	    #  latitude
	    foreach "x ${data}(long0) lamax lamin" \
		    [AverageExtrLatLong $latd $longd $datum $ps] {}
	    if { abs($lamax) < abs($lamin) } {
		set lamax $lamin
	    }
	    set ${data}(lat1) $lamax
	    if { $lamax < 0 } {
		set ${data}(polasp) south
	    } else { set ${data}(polasp) north }
	    set ${data}(lat0) 0
	}
	Merc1 {
	    # Mercator projection, single parallel
	    # take average of longitudes as central longitude, and use k0=1
	    set ${data}(long0) [AverageLong $longd $datum $ps]
	    set ${data}(k0) 1.0
	}
	Merc2 {
	    # Mercator projection, two parallels
	    # take average of longitudes as central longitude and extreme
	    #  latitude (in absolute value) as first standard parallel
	    foreach "x ${data}(long0) lamin lamax" \
		    [AverageExtrLatLong $latd $longd $datum $ps] {}
	    if { abs($lamin) > $lamax } { set lamax [expr abs($lamin)] }
	    set ${data}(lat1) $lamax
	}
	CS {
	    # Cassini-Soldner projection
	    # take averages of latitudes and longitudes as central coordinates
	    foreach "${data}(lat0) ${data}(long0)" \
		    [AverageLatLong $latd $longd $datum $ps] {}
	}
	APOLY {
	    # American Polyconic projection
	    # take averages of latitudes as standard latitude
	    set ${data}(lat0) [AverageLat $latd $datum $ps]
	}
	Schreiber {
	    # Schreiber projection
	    # all parameters are fixed, except the datum that must use the
	    #  "Bessel 1841" ellipsoid
	    if { [EllipsdOf $datum] != "Bessel 1841" } {
		GMMessage "Datum: Rijks Driehoeksmeting"
		set datum "Rijks Driehoeksmeting"
	    }
	    ProjSchreiberComputeAux $data $datum
	    set res [ProjSchreiberPoint $data $latd $longd $datum]
	    set askauxcomp 0
	}
	EOV {
	    # Hungarian grid projection: no parameters needed
	    set res [ProjEOVPoint {} $latd $longd $datum]
	}
	default {
	    BUG "no initialization proc for projection"
	    set askauxcomp 0
	}
    }
    if { $askauxcomp } {
	if { $ASKPROJPARAMS } {
	    ProjParams change $proj $data
	}
	Proj${proj}ComputeAux $data $datum
	set res [Proj${proj}Point $data $latd $longd $datum]
    }
    return $res
}

proc BadProjSetup {data proj latd longd datum} {
    # this must be kept in line with proc ProjInit (above)
    # initialize projection for use in command mode
    #  $data is array for the parameters with external parameters and
    #    datum defined, except in the case of "UTM" and "EOV" that only
    #    have datum
    # return 0 on success, 1 on error
    global ASKPROJPARAMS $data

    set ASKPROJPARAMS 0
    set aux 0
    switch $proj {
	UTM {
	    if { [catch {set ${data}(UTMzone)}] } {
		set p [DegreesToUTM $latd $longd $datum]
		scan [lindex $p 0] %0d ze
		set zn [lindex $p 1]
		set ${data}(UTMzone) [list $ze $zn]
		set ${data}(m_0) [expr -183+6.0*$ze]
	    }
	}
	LambNTF {
	    # as LCC2 but with zones with specific parameters
	    # no parameters can be changed by the user
	    set p [DegreesToLambNTF $latd $longd $datum]
	    if { [set z [lindex $p 0]] == "--" } {
		set z II
	    }
	    global LNTFz$z
	    catch {unset $data}
	    array set $data [array get LNTFz$z]
	}
	Schreiber {
	    # Schreiber projection
	    # all parameters are fixed, except the datum that must use the
	    #  "Bessel 1841" ellipsoid
	    if { [EllipsdOf [set ${data}(datum)]] != "Bessel 1841" } {
		set ${data}(datum) "Rijks Driehoeksmeting"
	    }
	}
	TM -  BMN -  CTR -  EOV -  GKK -  KKJP -  TWG {
	}
	LCC1 -  Stereogr -  SOM -  LCC2 -  AEA -
	LEAC -  Merc1 -  Merc2 -  CS -  APOLY {
	    set aux 1
	}
	default {
	    return 1
	}
    }
    if { $aux } {
	Proj${proj}ComputeAux $data [set ${data}(datum)]
    }
    return 0
}

proc AverageExtrLatLong {latd longd datum ps} {
    # compute the averages of latitudes and longitudes and the minimum and
    #  maximum latitude for $latd,$longd and all positions in $ps except
    #  the first
    # return list with the two averages followed by the maximum and the minimum

    set n 1
    set n 1
    set sla [set lamx [set lamn $latd]] ; set slo $longd
    foreach p [lreplace $ps 0 0] {
	foreach "la lo posdatum" $p { break }
	if { $posdatum != $datum } {
	    foreach "la lo" \
		    [ToDatum $la $lo $posdatum $datum] { break }
	}
	set sla [expr $sla+$la] ; set slo [expr $slo+$lo]
	if { $la > $lamx } { set lamx $la }
	if { $la < $lamn } { set lamn $la }
	incr n
    }
    return [list [expr 1.0*$sla/$n] [expr 1.0*$slo/$n] $lamx $lamn]
}

proc AverageLatLong {latd longd datum ps} {
    # compute the averages of latitudes and longitudes for $latd,$longd and
    #  all positions in $ps except the first
    # return list with the two averages

    set n 1
    foreach p [lreplace $ps 0 0] {
	foreach "la lo posdatum" $p { break }
	if { $posdatum != $datum } {
	    foreach "la lo" \
		    [ToDatum $la $lo $posdatum $datum] { break }
	}
	set latd [expr $latd+$la] ; set longd [expr $longd+$lo]
	incr n
    }
    return [list [expr 1.0*$latd/$n] [expr 1.0*$longd/$n]]
}

proc AverageLat {latd datum ps} {
    # compute the average of latitudes for $latd and
    #  all positions in $ps except the first

    set n 1
    foreach p [lreplace $ps 0 0] {
	set la [lindex $p 0] ; set pdatum [lindex $p 2]
	if { $pdatum != $datum } {
	    set la [lindex [ToDatum $la [lindex $p 1] $pdatum $datum] 0]
	}
	set latd [expr $la+$latd]
	incr n
    }
    return [expr 1.0*$latd/$n]
}

proc AverageLong {longd datum ps} {
    # compute the average of longitudes for $longd and
    #  all positions in $ps except the first

    set n 1
    foreach p [lreplace $ps 0 0] {
	set lo [lindex $p 1] ; set pdatum [lindex $p 2]
	if { $pdatum != $datum } {
	    set lo [lindex [ToDatum [lindex $p 0] $lo $pdatum $datum] 1]
	}
	set longd [expr $longd+$lo]
	incr n
    }
    return [expr 1.0*$longd/$n]
}

## converting longitudes to -180, 180 (inclusive) or -360, 360 (inclusive)

proc NormalLong {long} {
    # convert longitude in decimal degrees to the -180,180 (inclusive) range

    if { abs([set x [expr abs($long)*0.00555555555555555556]])-1 > 1e-14 } {
	set x [expr 0.5*($x+1)]
	set x [expr (($x-floor($x))-0.5)*360]
    } else { set x [expr abs($long)] }
    if { $long < 0 && abs($x-180) > 1e-14 } { return [expr -$x] }
    return $x
}

proc NormalLongCentred {long longc} {
    # convert longitude in decimal degrees to the closest value to $longc in
    #  the -360,360 (inclusive) range

    set long [NormalLong $long]
    if { $long <= $longc-180 } {
	set long [expr $long+360]
    } elseif { $long > $longc+180 } {
	set long [expr $long-360]
    }
    return $long
}

## projection procedures for UTM projection

# there may be auxiliary parameters to be computed from the main ones

proc ProjUTMComputeAux {data args} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    #  $args is not used but is needed because of uniform way of calling
    #   these procedures
    global $data 

    set z [set ${data}(UTMzone)]
    # set the zone to a list of number and letter if it is a string
    if { [regexp {^([0-9]+)([A-Z])$} $z x ze zn] } {
	set ${data}(UTMzone) [list $ze $zn]
    } else { set ze [lindex $z 0] }
    set ${data}(m_0) [expr -183+6.0*$ze]
    return
}

# projection proc is a function of latd,longd,datum to planar Cartesian
#  coordinates x,y in the terrain

proc ProjUTMPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    foreach "pze pzn px py" [DegreesToUTM $latd $longd $prdatum] {}
    foreach "mze mzn" [set ${data}(UTMzone)] {}
    if { $mze!=$pze || $mzn != $pzn } {
	foreach "px py" [CompUTMOnZone $latd $longd [set ${data}(m_0)]] {
	    break
	}
    }
    return [list $px $py]
}

# inverse projection proc is a function of planar Cartesian coordinates in the
#  terrain to list with latitude and longitude in the datum of the projection

proc ProjUTMInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    foreach "ze zn" [set ${data}(UTMzone)] {}
    set datum [set ${data}(datum)]
    set p [UTMToDegrees $ze $zn $x $y $datum]
    return $p
}

## Transverse Mercator projection

proc ProjTMPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    return [ConvToTM $latd $longd [set ${data}(lat0)] \
	    [set ${data}(long0)] [set ${data}(k0)] $prdatum]
}

proc ProjTMInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    set datum [set ${data}(datum)]
    set p [ConvFromTM $x $y [set ${data}(lat0)] [set ${data}(long0)] \
	    [set ${data}(k0)] $datum]
    return $p
}

## German Grid (Gauss-Krueger-Koordinatensystem) projection
# Information provided by Andreas Lange (Andreas.C.Lange@GMX.de)
# This is a Transverse Mercator projection with zone in [0-6], lat0=0,
#  lon0=zone*3 (0, 3, 6, 9, 12, 15E), and scale factor at central meridian
#  of k0=1.0

# Basic Finnish Grid (KKJP) projection
# Similar with zone in [1-4], lon0=zone*3+18 (21, 24, 27, 30E)

# Taiwan Grid projection
# Information provided by Dan Jacobson (jidanni@yahoo.com.tw)
# Similar with zone in [1-6], lon0=zone*2+113 (115, 117, ..., 125E),
#  and k0=0.9999

# Carta Tecnica Regionale (CTR), Italian projection
# "Le carte topografiche CTR ed il loro uso GPS"
# (http://www.gpscomefare.com/guide/tutorialgps/mapdatum.htm)
# May 2003 (information kindly sent by Alessandro Palmas)
# Similar with zone in [1-2], lon0=zone*6+3 (9, 15), k0=0.9996

# Austrian BMN grid projection
# information kindly sent by Alessandro Palmas, July 2003
# Similar with zones M28, M31, M34, lon0=zone_index*3+10.333333, k0=1

proc ProjGKKComputeAux {data args} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    #  $args is not used but is needed because of uniform way of calling
    #   these procedures
    global $data 

    set ${data}(m_0) 1.0
    return
}

proc ProjKKJPComputeAux {data args} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    #  $args is not used but is needed because of uniform way of calling
    #   these procedures
    global $data 

    set ${data}(m_0) 1.0
    return
}

proc ProjTWGComputeAux {data args} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    #  $args is not used but is needed because of uniform way of calling
    #   these procedures
    global $data 

    set ${data}(m_0) 0.9999
    return
}

proc ProjCTRComputeAux {data args} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    #  $args is not used but is needed because of uniform way of calling
    #   these procedures
    global $data 

    set ${data}(m_0) 0.9996
    return
}

proc ProjBMNComputeAux {data args} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    #  $args is not used but is needed because of uniform way of calling
    #   these procedures
    global $data 

    set ${data}(m_0) 1.0
    return
}

proc ProjGKKPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position

    return [Proj_TMZ_Point gkk $data $latd $longd $datum]
}

proc ProjKKJPPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position

    return [Proj_TMZ_Point kkjp $data $latd $longd $datum]
}

proc ProjTWGPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position

    return [Proj_TMZ_Point twg $data $latd $longd $datum]
}

proc ProjCTRPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position

    return [Proj_TMZ_Point ctr $data $latd $longd $datum]
}

proc ProjBMNPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position

    return [Proj_TMZ_Point bmn $data $latd $longd $datum]
}

proc Proj_TMZ_Point {pr data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    return [ConvToTM $latd $longd 0 [set ${data}(${pr}long0)] \
	             [set ${data}(m_0)] $prdatum]
}

proc ProjGKKInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain

    return [Proj_TMZ_Invert gkk $data $x $y]
}

proc ProjKKJPInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain

    return [Proj_TMZ_Invert kkjp $data $x $y]
}

proc ProjTWGInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain

    return [Proj_TMZ_Invert twg $data $x $y]
}

proc ProjCTRInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain

    return [Proj_TMZ_Invert ctr $data $x $y]
}

proc ProjBMNInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain

    return [Proj_TMZ_Invert bmn $data $x $y]
}

proc Proj_TMZ_Invert {pr data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    set datum [set ${data}(datum)]
    set p [ConvFromTM $x $y 0 [set ${data}(${pr}long0)] \
	              [set ${data}(m_0)] $datum]
    return $p
}

## Lambert Conic Conformal projections

# single standard parallel

proc ProjLCC1ComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data 

    set d [EllipsdData $datum]
    set a [lindex $d 0] ; set f [lindex $d 1]
    set es [expr $f*(2-$f)]
    set e [set ${data}(m_e) [expr sqrt($es)]]
    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
    if { [set n [set ${data}(lcc_n) [expr sin($phi0)]]] < 0 } {
	set ${data}(lcc_sn) -1
    } else { set ${data}(lcc_sn) 1 }
    set m0 [expr cos($phi0)/sqrt(1-$es*$n*$n)]
    set t0n [expr pow([ExpRedLat $phi0 $n $e],$n)]
    set F [expr $m0/($n*$t0n)]
    set aFk0 [set ${data}(m_a) [expr $a*$F*[set ${data}(k0)]]]
    set ${data}(lcc_rho0F) [expr $aFk0*$t0n]
    return
}

proc ProjLCC1Point {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set long0 [set ${data}(long0)]
    set longd [NormalLongCentred $longd $long0]
    set phi [expr $latd*0.01745329251994329576]
    set theta [expr [set ${data}(lcc_n)]* \
	            ($longd-$long0)*0.01745329251994329576]
    set t [ExpRedLat $phi [expr sin($phi)] [set ${data}(m_e)]]
    set rho [expr [set ${data}(m_a)]* \
	          pow($t,[set ${data}(lcc_n)])]
    set x [expr $rho*sin($theta)]
    set y [expr [set ${data}(lcc_rho0F)]-$rho*cos($theta)]
    return [list $x $y]
}

proc ProjLCC1Invert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    set z [expr [set ${data}(lcc_rho0F)]-$y]
    set thetap [expr atan(1.0*$x/$z)]
    set rhop [expr [set ${data}(lcc_sn)]*[Hypot $x $z]]
    set tp [expr pow($rhop/([set ${data}(m_a)]), \
	             1.0/[set ${data}(lcc_n)])]
    set longd [expr $thetap/[set ${data}(lcc_n)]*57.29577951308232087684+ \
	            [set ${data}(long0)]]
    set longd [NormalLong $longd]
    set latd [expr [LatFromRedLat $tp [set ${data}(m_e)]]* \
	           57.29577951308232087684]
    return [list $latd $longd]
}

# two standard parallels

proc ProjLCC2ComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data MESS

    set d [EllipsdData $datum]
    set a [lindex $d 0] ; set f [lindex $d 1]
    set es [expr $f*(2-$f)]
    set e [set ${data}(m_e) [expr sqrt($es)]]
    set phi1 [expr [set ${data}(lat1)]*0.01745329251994329576]
    set phi2 [expr [set ${data}(lat2)]*0.01745329251994329576]
    set phiF [expr [set ${data}(latF)]*0.01745329251994329576]
    set s1 [expr sin($phi1)] ; set s2 [expr sin($phi2)]
    set t1 [ExpRedLat $phi1 $s1 $e]
    set t2 [ExpRedLat $phi2 $s2 $e]
    if { [catch {set m1 [expr cos($phi1)/sqrt(1-$es*$s1*$s1)]}] || \
	    [catch {set m2 [expr cos($phi2)/sqrt(1-$es*$s2*$s2)]}] || \
	    [catch {set n [expr (log($m1)-log($m2))/(log($t1)-log($t2))]}] } {
	GMMessage $MESS(badProjargs)
	ProjParams change LCC2 $data
	ProjLCC2ComputeAux $data $datum
	return
    }
    set ${data}(lcc_n) $n
    if { $n < 0 } {
	set ${data}(lcc_sn) -1
    } else { set ${data}(lcc_sn) 1 }
    set F [expr $m1/($n*pow($t1,$n))]
    set aF [set ${data}(m_a) [expr $a*$F]]
    set tF [ExpRedLat $phiF [expr sin($phiF)] $e]
    set ${data}(lcc_rho0F) [expr $aF*pow($tF,$n)]
    return
}

proc ProjLCC2Point {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set longF [set ${data}(longF)]
    set longd [NormalLongCentred $longd $longF]
    set phi [expr $latd*0.01745329251994329576]
    set theta [expr [set ${data}(lcc_n)]* \
	            ($longd-$longF)*0.01745329251994329576]
    set t [ExpRedLat $phi [expr sin($phi)] [set ${data}(m_e)]]
    set rho [expr [set ${data}(m_a)]* \
	          pow($t,[set ${data}(lcc_n)])]
    set x [expr $rho*sin($theta)]
    set y [expr [set ${data}(lcc_rho0F)]-$rho*cos($theta)]
    return [list $x $y]
}

proc ProjLCC2Invert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    set z [expr [set ${data}(lcc_rho0F)]-$y]
    set thetap [expr atan(1.0*$x/$z)]
    set rhop [expr [set ${data}(lcc_sn)]*[Hypot $x $z]]
    set tp [expr pow($rhop/([set ${data}(m_a)]), \
	             1.0/[set ${data}(lcc_n)])]
    set longd [expr $thetap/[set ${data}(lcc_n)]*57.29577951308232087684+ \
	            [set ${data}(longF)]]
    set longd [NormalLong $longd]
    set latd [expr [LatFromRedLat $tp \
	             [set ${data}(m_e)]]*57.29577951308232087684]
    return [list $latd $longd]
}

# French NTF projection

proc ProjLambNTFComputeAux {data datum} {
    global $data

    set z [set $data(NTFzone)]
    global LNTFz$z
    array set $data [array get LNTFz$z]
    return
}

proc ProjLambNTFPoint {data latd longd datum} {

    return [ProjLCC2Point $data $latd $longd $datum]
}

proc ProjLambNTFInvert {data x y} {

    return [ProjLCC2Invert $data $x $y]
}

# ancillary procs for LCC projections

proc ExpRedLat {phi sinphi e} {
    # compute exp of reduced latitude for $phi
    #  $e is first eccentricity of ellipsoid

    set se [expr $sinphi*$e]
    return [expr tan((1.5707963267948966-$phi)/2.0) / \
	    pow((1.0-$se)/(1.0+$se), $e/2.0)]
}

proc LatFromRedLat {ts e} {
    # compute latitude from reduced latitude

    set eccnth [expr $e/2.0]
    set phi [expr 1.5707963267948966-2*atan($ts)]
    set i 16 ; set d 1
    while { abs($d) > 1e-10 && [incr i -1] } {
	set con [expr $e*sin($phi)]
	set d [expr 1.5707963267948966- \
		    2*atan($ts*pow((1-$con)/(1+$con),$eccnth))-$phi]
	set phi [expr $phi+$d]
    }
    return $phi
}

## Albers Equal Area Conic and Lambert Equal Area Conic projection
# adapted from PROJ4.3.3

proc ProjAEAComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data

    set phi1 [expr [set ${data}(lat1)]*0.01745329251994329576]
    set phi2 [expr [set ${data}(lat2)]*0.01745329251994329576]
    ComputeAux_AEA_LEAC AEA $data $datum $phi1 $phi2
    return
}

proc ProjAEAPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data MESS

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set long0 [set ${data}(long0)]
    set longd [NormalLongCentred $longd $long0]
    set a [set ${data}(m_a)] ; set n [set ${data}(m_1)]
    set phi [expr $latd*0.01745329251994329576]
    set lam [expr $n*($longd-$long0)*0.01745329251994329576]
    set rho [expr [set ${data}(m_3)]-$n* \
	            [SmallQ [expr sin($phi)] [set ${data}(m_e)] \
		            [set ${data}(m_0)]]]
    if { $rho < 0 } {
	GMMessage $MESS(outrngproj)
	return [list 0 0]
    }
    set rho [expr [set ${data}(m_4)]*sqrt($rho)]
    set x [expr $a*$rho*sin($lam)]
    set y [expr $a*([set ${data}(m_5)]-$rho*cos($lam))]
    return [list $x $y]
}

proc ProjAEAInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data MESS TXT

    set a [set ${data}(m_a)] ; set n [set ${data}(m_1)]
    set datum [set ${data}(datum)]
    set x [expr $x/$a]
    set y [expr [set ${data}(m_5)]-$y/$a]
    if { [set rho [Hypot $x $y]] == 0 } {
	set lam 0
	if { $n > 0 } {
	    set phi 1.5707963267948966
	} else { set phi -1.5707963267948966 }
    } else {
	if { $n < 0 } {
	    set rho [expr -$rho] ; set x [expr -$x] ; set y [expr -$y]
	}
	set phi [expr 1.0*$rho/[set ${data}(m_4)]]
	set phi [expr 1.0*([set ${data}(m_3)]-$phi*$phi)/$n]
	if { abs([set ${data}(m_2)]-abs($phi)) > 1e-7 } {
	    set phi [LatAngle $phi [set ${data}(m_e)] \
		              [set ${data}(m_0)]]
	    if { $phi > 1e19 } {
		## return 0 180 on error...
		GMMessage [format $MESS(badinvproj) $TXT(PRJAEA)/$TXT(PRJLEAC)]
		return [list 0 180]
	    }
	} elseif { $phi < 0 } {
	    set phi -1.5707963267948966
	} else { set phi 1.5707963267948966 }
	set lam [expr 1.0*atan2($x,$y)/$n]
    }
    set latd [expr $phi*57.29577951308232087684]
    set longd [expr ($lam*57.29577951308232087684)+[set ${data}(long0)]]
    return [list $latd [NormalLong $longd]]
}

proc ProjLEACComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data

    set phi2 [expr [set ${data}(lat1)]*0.01745329251994329576]
    if { [set ${data}(polasp)] == "south" } {
	set phi1 -1.5707963267948966
    } else { set phi1 1.5707963267948966 }
    ComputeAux_AEA_LEAC LEAC $data $datum $phi1 $phi2
    return
}

proc ProjLEACPoint {data latd longd datum} {

    return [ProjAEAPoint $data $latd $longd $datum]
}

proc ProjLEACInvert {data x y} {

    return [ProjAEAInvert $data $x $y]
}

# ancillary procedures for AEA and LEAC projections

proc ComputeAux_AEA_LEAC {proj data datum phi1 phi2} {
    global $data MESS

    if { abs($phi1+$phi2) < 1e-10 } {
	GMMessage $MESS(badProjargs)
	ProjParams change $proj $data
	ProjAEAComputeAux $data $datum
	return
    }
    set d [EllipsdData $datum]
    set ${data}(m_a) [lindex $d 0] ; set f [lindex $d 1]
    set es [expr $f*(2-$f)]
    set e [set ${data}(m_e) [expr sqrt($es)]]
    set oes [set ${data}(m_0) [expr 1-$es]]
    set n [set ${data}(m_1) [expr sin($phi1)]]
    set cosp [expr cos($phi1)]
    set m1 [expr $cosp/sqrt(1.0-$es*$n*$n)]
    set ml1 [SmallQ $n $e $oes]
    if { abs($phi1-$phi2) > 1e-10 } {
	# secant
	set s2 [expr sin($phi2)]
	set m2 [expr cos($phi2)/sqrt(1.0-$es*$s2*$s2)]
	set ml2 [SmallQ $s2 $e $oes]
	set n [set ${data}(m_1) [expr ($m1*$m1-$m2*$m2)/($ml2-$ml1)]]
    }
    set ${data}(m_2) [expr 1.0-0.5*$oes*log((1.0-$e)/(1.0+$e))/$e]
    set c [set ${data}(m_3) [expr $m1*$m1+$n*$ml1]]
    set dd [set ${data}(m_4) [expr 1.0/$n]]
    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
    set ${data}(m_5) [expr $dd*sqrt($c-$n*[SmallQ sin($phi0) $e $oes])]
    return
}

proc LatAngle {qs te toes} {
    # compute latitude angle phi-1

    set phi [Aasin [expr 0.5*$qs]]
    if { $te < 1e-7 } { return $phi }
    for { set i 15 } { $i > 0 } { incr i -1 } {
	set sinpi [expr sin($phi)] ; set cospi [expr cos($phi)]
	set con [expr $te*$sinpi]
	set com [expr 1.0-$con*$con]
	set dphi [expr 0.5*$com*$com/$cospi* \
		       ($qs/$toes-$sinpi/$com+0.5/$te* \
		        log((1.0-$con)/(1.0+$con)))]
	set phi [expr $phi+$dphi]
	if { abs($dphi) < 1e-10 } { return $phi }
    }
    return 1e20
}

proc SmallQ {sinphi e oes} {
    # compute small q

    if { $e < 1e-7 } { return [expr $sinphi+$sinphi] }
    set con [expr $e*$sinphi]
    return [expr $oes*($sinphi/(1.0-$con*$con)- \
	               (0.5/$e)*log((1.0-$con)/(1.0+$con)))]
}

## Mercator projection: a special case of Lambert Conic Conformal
#   with the equator as standard parallel

# single standard parallel

proc ProjMerc1ComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data 

    set d [EllipsdData $datum]
    set a [lindex $d 0] ; set f [lindex $d 1]
    set ${data}(m_a) [expr 1.0*$a*[set ${data}(k0)]]
    set es [expr $f*(2-$f)]
    set ${data}(m_e) [expr sqrt($es)]
    set e4 [expr $es*$es] ; set e6 [expr $e4*$es]
    set ${data}(m_1) [expr $es* \
	    (0.5+0.20833333333333333333*$es+0.08333333333333333333*$e4+ \
             0.03611111111111111111*$e6)]
    set ${data}(m_2) [expr $es* \
	    (0.14583333333333333333*$es+0.12083333333333333333*$e4+ \
	     0.07039930555555555555*$e6)]
    set ${data}(m_3) [expr $es* \
	    (0.05833333333333333333*$e4+0.07232142857142857142*$e6)]
    set ${data}(m_4) [expr 0.02653149801587301587*$es*$e6]
    return
}

proc ProjMerc1Point {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set long0 [set ${data}(long0)]
    set longd [NormalLongCentred $longd $long0]
    set phi [expr $latd*0.01745329251994329576]
    set se [expr sin($phi)*[set ${data}(m_e)]]
    set mm [expr tan((1.5707963267948966+$phi)/2.0)* \
     	         pow((1.0-$se)/(1.0+$se), [set ${data}(m_e)]/2.0)]
    set x [expr [set ${data}(m_a)]*($longd-$long0)*0.01745329251994329576]
    set y [expr [set ${data}(m_a)]*log($mm)]
    return [list $x $y]
}

proc ProjMerc1Invert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    set ksi [expr 1.5707963267948966- \
	          2*atan(exp(-$y/[set ${data}(m_a)]))]
    set latd [expr ($ksi+[set ${data}(m_1)]*sin($ksi+$ksi)+ \
	            [set ${data}(m_2)]*sin(4*$ksi)+ \
	            [set ${data}(m_3)]*sin(6*$ksi)+ \
	            [set ${data}(m_4)]*sin(8*$ksi))* \
		   57.29577951308232087684]
    set longd [expr $x/[set ${data}(m_a)]*57.29577951308232087684+ \
	            [set ${data}(long0)]]
    return [list $latd [NormalLong $longd]]
}

# two standard parallels (symmetrical with respect to the equator)

proc ProjMerc2ComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data 

    set d [EllipsdData $datum]
    set a [lindex $d 0] ; set f [lindex $d 1]
    set es [expr $f*(2-$f)]
    set ${data}(m_e) [expr sqrt($es)]
    set phi1 [expr abs([set ${data}(lat1)])*0.01745329251994329576]
    set s [expr sin($phi1)]
    set k0 [expr cos($phi1)/sqrt(1-$es*$s*$s)]
    set ${data}(m_a) [expr 1.0*$a*$k0]
    set e4 [expr $es*$es] ; set e6 [expr $e4*$es]
    set ${data}(m_1) [expr $es* \
	    (0.5+0.20833333333333333333*$es+0.08333333333333333333*$e4+ \
             0.03611111111111111111*$e6)]
    set ${data}(m_2) [expr $es* \
	    (0.14583333333333333333*$es+0.12083333333333333333*$e4+ \
	     0.07039930555555555555*$e6)]
    set ${data}(m_3) [expr $es* \
	    (0.05833333333333333333*$e4+0.07232142857142857142*$e6)]
    set ${data}(m_4) [expr 0.02653149801587301587*$es*$e6]
    return
}

proc ProjMerc2Point {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    return [ProjMerc1Point $data $latd $longd $datum]
}

proc ProjMerc2Invert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    return [ProjMerc1Invert $data $x $y]
}

## Cassini-Soldner projection

proc ProjCSComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data 

    set d [EllipsdData $datum]
    set a [set ${data}(m_a) [lindex $d 0]] ; set f [lindex $d 1]
    set es [set ${data}(m_e) [expr $f*(2-$f)]]
    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
    set e4 [expr $es*$es] ; set e6 [expr $es*$e4]
    set m1 [set ${data}(m_1) [expr 1-0.25*$es-0.046875*$e4-0.01953125*$e6]]
    set m2 [set ${data}(m_2) [expr 0.375*$es+0.09375*$e4+0.0439453125*$e6]]
    set m3 [set ${data}(m_3) [expr 0.05859375*$e4+0.0439453125*$e6]]
    set m4 [set ${data}(m_4) [expr 0.01139322916666666667*$e6]]
    set ${data}(m_0) [expr $a*($m1*$phi0-$m2*sin($phi0+$phi0)+ \
	                          $m3*sin(4*$phi0)-$m4*sin(6*$phi0))]
    set se [expr sqrt(1-$es)]
    set e1 [expr (1.0-$se)/(1+$se)] ; set e12 [expr $e1*$e1]
    set ${data}(m_5) [expr $e1*(1.5-0.84375*$e12)]
    set ${data}(m_6) [expr $e12*(1.3125-1.71875*$e12)]
    set ${data}(m_7) [expr 1.57291666666666666667*$e12*$e1]
    set ${data}(m_8) [expr 2.142578125*$e12*$e12]
    return
}

proc ProjCSPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set long0 [set ${data}(long0)]
    set longd [NormalLongCentred $longd $long0]
    set a [set ${data}(m_a)] ; set es [set ${data}(m_e)]
    set phi [expr $latd*0.01745329251994329576]
    set cphi [expr cos($phi)] ; set sphi [expr sin($phi)]
    set A [expr ($longd-$long0)*0.01745329251994329576*$cphi]
    set A2 [expr $A*$A] ; set A4 [expr $A2*$A2]
    set t [expr tan($phi)] ; set T [expr $t*$t]
    set C [expr $es*$cphi*$cphi/(1.0-$es)]
    set v [expr $a/sqrt(1-$es*$sphi*$sphi)]
    set M [expr $a*([set ${data}(m_1)]*$phi- \
	            [set ${data}(m_2)]*sin($phi+$phi)+ \
	            [set ${data}(m_3)]*sin(4*$phi)- \
                    [set ${data}(m_4)]*sin(6*$phi))]
    set x [expr $v*$A*(1-$T*$A2/6.0-(8-$T+8*$C)*$T*$A4/120.0)]
    set y [expr $M-[set ${data}(m_0)]+$v*$t*(0.5*$A2+(5-$T+6*$C)*$A4/24.0)]
    return [list $x $y]
}

proc ProjCSInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data

    set a [set ${data}(m_a)] ; set es [set ${data}(m_e)]
    set miu1 [expr ([set ${data}(m_0)]+$y)/($a*[set ${data}(m_1)])]
    set phi1 [expr $miu1+[set ${data}(m_5)]*sin($miu1+$miu1)+ \
	           [set ${data}(m_6)]*sin(4*$miu1)+ \
	           [set ${data}(m_7)]*sin(6*$miu1)+ \
	           [set ${data}(m_8)]*sin(8*$miu1)]
    set sphi1 [expr sin($phi1)]
    set t [expr 1-$es*$sphi1*$sphi1]
    set v1 [expr $a/sqrt($t)] ; set rho1 [expr $v1*(1-$es)/$t]
    set t [expr tan($phi1)] ; set T1 [expr $t*$t]
    set D [expr $x/$v1] ; set D2 [expr $D*$D]
    set tt [expr (1+3*$T1)*$D2]
    set lat [expr ($phi1-$v1*$t/$rho1*$D2*(0.5-$tt/24.0))* \
	          57.29577951308232087684]
    set long [expr $D*(1-$T1*$D2/3.0+$tt*$D2/15.0)/cos($phi1)* \
	          57.29577951308232087684+[set ${data}(long0)]]
    return [list $lat [NormalLong $long]]
}

## American polyconic projection
# adapted from PROJ4.0

proc ProjAPOLYComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data

    set d [EllipsdData $datum]
    set ${data}(m_a) [lindex $d 0]
    set f [lindex $d 1]
    set es [set ${data}(m_5) [expr $f*(2-$f)]]
    MeridDistParams $data $es
    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
    set ${data}(m_6) \
	    [MeridionalDist $data $phi0 [expr sin($phi0)] [expr cos($phi0)]]
    return
}

proc ProjAPOLYPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set a [set ${data}(m_a)]
    set phi [expr $latd*0.01745329251994329576]
    set lam [expr $longd*0.01745329251994329576]
    if { abs($phi) <= 1e-10 } {
	set x [expr $lam*$a] ; set y [expr -$a*[set ${data}(m_6)]]
    } else {
	set sp [expr sin($phi)]
	if { abs([set cp [expr cos($phi)]]) > 1e-10 } {
	    set ms [expr $cp/sqrt(1-[set ${data}(m_5)]*$sp*$sp)/$sp]
	} else { set ms 0 }
	set t [expr $lam*$sp]
	set x [expr $a*$ms*sin($t)]
	set y [expr $a*(([MeridionalDist $data $phi $sp $cp]- \
		[set ${data}(m_6)])+$ms*(1-cos($t)))]
    }
    return [list $x $y]
}

proc ProjAPOLYInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data MESS TXT

    set a [set ${data}(m_a)]
    set x [expr 1.0*$x/$a] ; set y [expr 1.0*$y/$a]
    set y [expr $y+[set ${data}(m_6)]]
    if { abs($y) <= 1e-10 } {
	set phi 0 ; set lam $x
    } else {
	set r [expr $x*$x+$y*$y] ; set phi $y
	set es [set ${data}(m_5)]
	for { set i 0 } { $i < 20 } { incr i } {
	    set sp [expr sin($phi)] ; set cp [expr cos($phi)]
	    set s2ph [expr $sp*$cp]
	    if { abs($cp) < 1e-12 } {
		GMMessage [format $MESS(badinvproj) $TXT(PRJAPOLY)]
		set i 0
		break
	    }
	    set mlp [expr sqrt(1-$es*$sp*$sp)]
	    set c [expr $sp*$mlp/$cp]
	    set ml [MeridionalDist $data $phi $sp $cp]
	    set mlb [expr $ml*$ml+$r]
	    set mlp [expr (1-$es)/($mlp*$mlp*$mlp)]
	    set dphi [expr ($ml+$ml+$c*$mlb-2.0*$y*($c*$ml+1))/ \
		    ($es*$s2ph*($mlb-2*$y*$ml)/$c+ \
		     2*($y-$ml)*($c*$mlp-1.0/$s2ph)-$mlp-$mlp)]
	    set phi [expr $phi+$dphi]
	    if { abs($dphi) <= 1e-12 } { break }
	}
	if { $i == 0 } {
	    GMMessage [format $MESS(badinvproj) $TXT(PRJAPOLY)]
	}
	set c [expr sin($phi)]
	set lam [expr [Aasin [expr $x*tan($phi)*sqrt(1-$es*$c*$c)]]/sin($phi)]
    }
    set longd [expr $lam*57.29577951308232087684]
    set latd [expr $phi*57.29577951308232087684]
    return [list $latd $longd]
}

## auxiliary procs for polyconic projection
# meridional distance for ellipsoid and inverse

proc MeridDistParams {data es} {
    # compute parameters for meridional distance and inverse functions
    #  saving them at indices m_0, ..., m_4 of array $data
    global $data

    set t [expr $es*(0.046875+$es*(0.01953125+$es*0.01068115234375))]
    set ${data}(m_0) [expr 1-$es*(0.25+$t)]
    set ${data}(m_1) [expr $es*(0.75-$t)]
    set t [expr $es*$es]
    set ${data}(m_2) [expr $t*(0.46875-$es* \
	    (0.01302083333333333333+$es*0.00712076822916666666))]
    set t [expr $es*$t]
    set ${data}(m_3) \
	    [expr $t*(0.36458333333333333333-$es*0.00569661458333333333)]
    set ${data}(m_4) [expr $t*$es*0.3076171875]
    return
}

proc MeridionalDist {data phi sphi cphi} {
    # compute meridional distance for given latitude
    # assume parameters for transformation are in array $data at indices
    #  m_0, ..., m_4
    global $data

    set cphi [expr $cphi*$sphi]
    set sphi [expr $sphi*$sphi]
    return [expr [set ${data}(m_0)]*$phi-$cphi* \
	    ([set ${data}(m_1)]+$sphi* \
	     ([set ${data}(m_2)]+$sphi* \
	      ([set ${data}(m_3)]+$sphi*[set ${data}(m_4)])))]
}

proc LatFromMeridDist {data d} {
    # compute latitude from given meridional distance
    # assume parameters for transformation and e^2 are in array $data at
    #  indices m_0, ..., m_4 and m_5
    global $data MESS

    set es [set ${data}(m_5)]
    set k [expr 1.0/(1-$es)]    
    set phi $d
    for { set i 0 } { $i < 10 } { incr i } {
	set s [expr sin($phi)]
	set t [expr 1.0-$es*$s*$s]
	set dphi [expr [MeridionalDist $data $phi $s [expr cos($phi)]]* \
		($t*sqrt($t))*$k]
	set phi [expr $phi-$dphi]
	if { abs($dphi) < 1e-11 } { return $phi }
    }
    GMMessage $MESS(badinvmdist)
    return $phi
}

## Stereographic projection
#

proc ProjStereogrComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data

    set d [EllipsdData $datum]
    set ${data}(m_a) [lindex $d 0] ; set f [lindex $d 1]
    set e [set ${data}(m_e) [expr sqrt($f*(2-$f))]]
    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
    set k0 [set ${data}(k0)]
    set t [expr abs($phi0)]
    if { abs($t-1.5707963267948966) < 1e-10 } {
	if { $phi0 < 0 } {
	    set mode s_pole
	} else { set mode n_pole }
	set ${data}(m_0) [expr 2*$k0/sqrt(pow(1+$e,1+$e)*pow(1-$e,1-$e))]
    } elseif { $t > 1e-10 } {
	set mode obliq
	set t [expr sin($phi0)]
	set X [expr 2*atan([SSFunc_ $phi0 $t $e])-1.5707963267948966]
	set t [expr $t*$e]
	set ${data}(m_0) [expr 2*$k0*cos($phi0)/sqrt(1.0-$t*$t)]
	set ${data}(m_1) [expr sin($X)]
	set ${data}(m_2) [expr cos($X)]
    } else {
	set mode equat
	set ${data}(m_0) [expr 2*$k0]
	set ${data}(m_1) 0
	set ${data}(m_2) 1.0
    }
    set ${data}(m_3) $mode
    return
}

proc ProjStereogrPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    set long0 [set ${data}(long0)]
    set longd [NormalLongCentred $longd $long0]
    set e [set ${data}(m_e)]
    set akm1 [set ${data}(m_0)]
    set phi [expr $latd*0.01745329251994329576]
    set lam [expr ($longd-$long0)*0.01745329251994329576]
    set coslam [expr cos($lam)] ; set sinlam [expr sin($lam)]
    set sinphi [expr sin($phi)]
    switch [set ${data}(m_3)] {
	obliq -
	equat {
	    set X [expr 2*atan([SSFunc_ $phi $sinphi $e])-1.5707963267948966]
	    set sinX [expr sin($X)] ; set cosX [expr cos($X)]
	    set sinX1 [set ${data}(m_1)]
	    set cosX1 [set ${data}(m_2)]
	    set A [expr $akm1/($cosX1*(1+$sinX1*$sinX+$cosX1*$cosX*$coslam))]
	    set y [expr $A*($cosX1*$sinX-$sinX1*$cosX*$coslam)]
	    set x [expr $A*$cosX]
	}
	n_pole {
	    set x [expr $akm1*[ExpRedLat $phi $sinphi $e]]
	    set y [expr -$x*$coslam]
	}
	s_pole {
	    set x [expr $akm1*[ExpRedLat [expr -$phi] [expr -$sinphi] $e]]
	    set y [expr $x*$coslam]
	}
    }
    set a [set ${data}(m_a)]
    set x [expr $x*$sinlam*$a] ; set y [expr $y*$a]
    return [list $x $y]
}

proc ProjStereogrInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data MESS TXT

    set e [set ${data}(m_e)] ; set a [set ${data}(m_a)]
    set akm1 [set ${data}(m_0)]
    set long0 [set ${data}(long0)] ; set datum [set ${data}(datum)]
    set x [expr $x/$a] ; set y [expr $y/$a]
    if { abs($x)<1e-6 && abs($y)<1e-6 } {
	return [list [set ${data}(lat0)] $long0]
    }
    set rho [Hypot $x $y]
    set spole 0
    switch [set m [set ${data}(m_3)]] {
	obliq -
	equat {
	    set sinX1 [set ${data}(m_1)]
	    set cosX1 [set ${data}(m_2)]
	    set tp [expr 2*atan2($rho*$cosX1,$akm1)]
	    set cosphi [expr cos($tp)] ; set sinphi [expr sin($tp)]
	    if { $rho == 0.0 } {
		set phi_l [Aasin [expr $cosphi*$sinX1]]
	    } else {
		set phi_l [Aasin [expr $cosphi*$sinX1+($y*$sinphi*$cosX1/$rho)]]
	    }
	    set tp [expr tan(0.5*(1.5707963267948966+$phi_l))]
	    set x [expr $x*$sinphi]
	    set y [expr $rho*$cosX1*$cosphi-$y*$sinX1*$sinphi]
	    set halfpi 1.5707963267948966
	    set halfe [expr 0.5*$e]
	}
	n_pole {
	    set y [expr -$y]
	    set tp [expr -1.0*$rho/$akm1]
	    set phi_l [expr 1.5707963267948966-2*atan($tp)]
	    set halfpi -1.5707963267948966
	    set halfe [expr -0.5*$e]
	}
	s_pole {
	    set spole 1
	    set tp [expr -1.0*$rho/$akm1]
	    set phi_l [expr 1.5707963267948966-2*atan($tp)]
	    set halfpi -1.5707963267948966
	    set halfe [expr -0.5*$e]
	}
    }
    for { set i 8 } { [incr i -1] } { set phi_l $phi } {
	set sinphi [expr $e*sin($phi_l)]
	set phi [expr 2*atan($tp*pow((1+$sinphi)/(1-$sinphi),$halfe))-$halfpi]
	if { abs($phi_l-$phi) < 1e-10 } {
	    if { $spole } { set phi [expr -$phi] }
	    if { $x == 0 && $y == 0 } {
		set longd 0
	    } else { set longd [expr atan2($x,$y)*57.29577951308232087684] }
	    set latd [expr $phi*57.29577951308232087684]
	    set longd [expr $longd+$long0]
	    return [list $latd $longd]
	}
    }
    ## return 0 180 on error...
    GMMessage [format $MESS(badinvproj) $TXT(PRJStereogr)]
    return [list 0 180]
}

# Schreiber double projection: used for the RDG (Netherlands) grid
#  this is a special case of the oblique Stereographic projection
#

proc ProjSchreiberComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data

    set phi0 [set ${data}(m_0) 0.91029672689323934677]
    set lam0 [set ${data}(m_1) 0.09403203751960005358]
    set e [set ${data}(m_2) 0.08169683122252750299]
    set k0 0.9999079
    # 2*k*R
    set ${data}(m_3) [expr 12765289.142*$k0]
    set n 1.00047585668 ; set m 0.003773953832
    set tau0 [expr log([SSFunc_ $phi0 [expr sin($phi0)] $e])]
    set B0 [expr 2*atan(pow(2.7182818284590452354,$n*$tau0+$m))- \
	    1.5707963267948966]
    set ${data}(m_4) [expr sin($B0)]
    set ${data}(m_5) [expr cos($B0)]
    return
}

proc ProjSchreiberPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    # just for the sake of those not knowing that this projection
    #  was designed for use in the Netherlands...
    set long0 5.38763888888888888564
    set longd [NormalLongCentred $longd $long0]
    set n 1.00047585668 ; set m 0.003773953832
    set e [set ${data}(m_2)]
    set phi [expr $latd*0.01745329251994329576]
    set lam [expr $n*($longd*0.01745329251994329576-[set ${data}(m_1)])]
    set coslam [expr cos($lam)] ; set sinlam [expr sin($lam)]
    set sinphi [expr sin($phi)]
    set tau [expr log([SSFunc_ $phi $sinphi $e])]
    set B [expr 2*atan(pow(2.7182818284590452354,$n*$tau+$m))- \
	    1.5707963267948966]
    set sinB [expr sin($B)] ; set cosB [expr cos($B)]
    set sinB0 [set ${data}(m_4)]
    set cosB0 [set ${data}(m_5)]
    set tkR [set ${data}(m_3)]
    set A [expr $tkR/(1+$sinB*$sinB0+$cosB*$cosB0*$coslam)]
    set x [expr $A*$sinlam*$cosB]
    set y [expr $A*($sinB*$cosB0-$cosB*$sinB0*$coslam)]
    return [list $x $y]
}

proc ProjSchreiberInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data MESS TXT

    set datum [set ${data}(datum)]
    set long0 [expr [set ${data}(m_1)]*57.29577951308232087684]
    if { abs($x)<1e-6 && abs($y)<1e-6 } {
	set lat0 [expr [set ${data}(m_0)]*57.29577951308232087684]
	return [list $lat0 $long0]
    }
    set n 1.00047585668 ; set m 0.003773953832
    set sinB0 [set ${data}(m_4)]
    set cosB0 [set ${data}(m_5)]
    set tkR [set ${data}(m_3)]
    set rho [Hypot $x $y]
    set Psi [expr 2*atan2($rho,$tkR)]
    set sinPsi [expr sin($Psi)]
    set B [Aasin [expr cos($Psi)*$sinB0+$y*$sinPsi*$cosB0/$rho]]
    set lam [Aasin [expr $x*$sinPsi/$rho/cos($B)]]
    set longd [expr $lam/$n*57.29577951308232087684+$long0]
    set tau [expr (log(tan(0.5*(1.5707963267948966+$B)))-$m)/$n]
    set etau [expr pow(2.7182818284590452354,$tau)]
    set phi [expr 2*atan($etau)-1.5707963267948966]
    set e [set ${data}(m_2)] ; set he [expr 0.5*$e]
    for { set i 0 } { $i < 50 } { incr i } {
	set esp [expr $e*sin($phi)]
	set phi_ [expr 2*atan($etau*pow((1+$esp)/(1-$esp),$he))- \
		1.5707963267948966]
	if { abs($phi_-$phi) < 1e-10 } {
	    set latd [expr $phi*57.29577951308232087684]
	    return [list $latd $longd]
	}
	set phi $phi_
    }
    ## return 0 180 on error...
    GMMessage [format $MESS(badinvproj) $TXT(PRJSchreiber)]
    return [list 0 180]
}

## ancillary procs for Stereographic and Schreiber projections

proc SSFunc_ {phit sinphi eccen} {

    set sinphi [expr $sinphi*$eccen]
    return [expr tan(0.5*(1.5707963267948966+$phit))* \
	    pow((1-$sinphi)/(1+$sinphi),0.5*$eccen)]
}

# Swiss Oblique Mercator projection

proc ProjSOMComputeAux {data datum} {
    # compute auxiliary parameters from main parameters
    #  $data is name of global array for the parameters
    global $data

    set d [EllipsdData $datum]
    set a [lindex $d 0] ; set f [lindex $d 1]
    set es [expr $f*(2-$f)]
    set e [set ${data}(m_e) [expr sqrt($es)]]
    set he [set ${data}(m_0) [expr 0.5*$e]]
    set oes [expr 1-$es]
    set roes [set ${data}(m_6) [expr 1.0/$oes]]
    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
    set cp [expr cos($phi0)] ; set cp [expr $cp*$cp]
    set c [set ${data}(m_1) [expr sqrt(1+$es*$cp*$cp*$roes)]]
    set sp [expr sin($phi0)]
    set sp0 [set ${data}(m_3) [expr $sp/$c]]
    set p0 [Aasin $sp0]
    set ${data}(m_2) [expr cos($p0)]
    set sp [expr $sp*$e]
    set ${data}(m_4) [expr log(tan(0.78539816339744833+0.5*$p0))- \
	    $c*(log(tan(0.78539816339744833+0.5*$phi0))- \
	        $he*log((1.0+$sp)/(1.0-$sp)))]
    set ${data}(m_5) [expr $a*[set ${data}(k0)]*sqrt($oes)/(1.0-$sp*$sp)]
    return
}

proc ProjSOMPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    global $data

    set prdatum [set ${data}(datum)]
    if { $datum != $prdatum } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum $prdatum] { break }
    }
    # just for the sake of those not knowing that this projection
    #  was designed for use in Switzerland...
    set long0 [set ${data}(long0)]
    set longd [NormalLongCentred $longd $long0]
    set phi [expr $latd*0.01745329251994329576]
    set sp [expr [set ${data}(m_e)]*sin($phi)]
    set c [set ${data}(m_1)]
    set phip [expr 2*atan(exp($c*(log(tan(0.78539816339744833+0.5*$phi))- \
	    [set ${data}(m_0)]*log((1.0+$sp)/(1.0-$sp)))+ \
	    [set ${data}(m_4)]))-1.5707963267948966]
    set lam [expr $c*($longd-$long0)*0.01745329251994329576]
    set cp [expr cos($phip)]
    set phipp [Aasin [expr [set ${data}(m_2)]*sin($phip)- \
	    [set ${data}(m_3)]*$cp*cos($lam)]]
    set akR [set ${data}(m_5)]
    set x [expr $akR*[Aasin [expr $cp*sin($lam)/cos($phipp)]]]
    set y [expr $akR*log(tan(0.78539816339744833+0.5*$phipp))]
    return [list $x $y]
}

proc ProjSOMInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    global $data MESS TXT

    set datum [set ${data}(datum)]
    set akR [set ${data}(m_5)]
    set phipp [expr 2*(atan(exp(1.0*$y/$akR))-0.78539816339744833)]
    set lampp [expr 1.0*$x/$akR]
    set cp [expr cos($phipp)]
    set phip [Aasin [expr [set ${data}(m_2)]*sin($phipp)+ \
	    [set ${data}(m_3)]*$cp*cos($lampp)]]
    set lamp [Aasin [expr $cp*sin($lampp)/cos($phip)]]
    set e [set ${data}(m_e)]
    set he [set ${data}(m_0)] ; set c [set ${data}(m_1)]
    set roes [set ${data}(m_6)]
    set con [expr ([set ${data}(m_4)]- \
	    log(tan(0.78539816339744833+0.5*$phip)))/$c]
    set i 10
    while { [incr i -1] } {
	set esp [expr $e*sin($phip)]
	set delp [expr ($con+log(tan(0.78539816339744833+0.5*$phip))- \
		$he*log((1.0+$esp)/(1.0-$esp)))*(1.0-$esp*$esp)*cos($phip)* \
		$roes]
	set phip [expr $phip-$delp]
	if { abs($delp) < 1e-10 } { break }
    }
    if { $i } {
	set latd [expr $phip*57.29577951308232087684]
	set longd [expr 57.29577951308232087684*$lamp/$c+ \
		[set ${data}(long0)]]
	return [list $latd $longd]
    }
    ## return 0 180 on error...
    GMMessage [format $MESS(badinvproj) $TXT(PRJSOM)]
    return [list 0 180]
}

# ancillary procedures for Swiss Oblique Mercator and others

proc Aasin {v} {
    # arcsin function
    global MESS

    set av [expr abs($v)]
    if { $av >= 1 } {
	if { $av > 1.00000000000001 } {
	    GMMessage [format $MESS(badargtofunc) asin]
	}
	if { $v < 0 } {
	    return -1.5707963267948966
	}
	return 1.5707963267948966
    }
    return [expr asin($v)]
}

# ancillary procedures for LCC1, LCC2, Stereographic and Schreiber projections

proc Hypot {x y} {
    # computation of sqrt($x*$x+$y*$y) avoiding overflows

    if { $x < 0 } {
	set x [expr -$x]
    } elseif { $x == 0 } {
	if { $y < 0 } { return [expr -$y] }
	return $y
    }
    if { $y < 0 } {
	set y [expr -$y]
    } elseif { $y == 0 } { return $x }
    if { $x < $y } {
	set x [expr 1.0*$x/$y]
	return [expr $y*sqrt(1+$x*$x)]
    }
    set y [expr 1.0*$y/$x]
    return [expr $x*sqrt(1+$y*$y)]
}

#
# SL contribution
#
# Implementation of the EOV (Hungarian National Projection)
#
# 2008-07-15
#
# Based on publications by Jozsef Varga Phd.
#   http://www.agt.bme.hu/staff_h/varga/Osszes/Dok3uj.htm
#   http://www.agt.bme.hu/staff_h/varga/gps/kezdoknek.html
#

proc ProjEOVPoint {data latd longd datum} {
    # compute planar Cartesian coordinates for given position
    # although this projection computes x for northing and
    #  y for easting, for compatibility with the mapping procedures
    #  this procedure returns the usual (easting, northing) pair,
    #  and the grid coordinates are computed by inversing it
    #  $data is not used

    # MF change: use the 3 parameter conversion to the HD72 if needed
    if { $datum != "Hungarian Datum 1972" } {
	foreach "latd longd" \
		[ToDatum $latd $longd $datum "Hungarian Datum 1972"] { break }
    }
    # MF changes after this point: for efficiency, original code commented out
    # HD72 Phi, Lambda to radians
    set Phi [expr $latd*0.01745329251994329576]
    set Lambda [expr $longd*0.01745329251994329576]
    #---

    # HD72 Phi, Lambda to spherical phi, lambda
    ## GRS67 semi-major / semi-minor axis
    # set a 6378160.0
    # set b 6356774.516
    ## EOV parameters
    set k 1.003110007693
    set n 1.000719704936
    # set Lambda_0 [expr (19.0 + 02.0/60.0 + 54.8584/3600.0) / 360.0 * 2.0*$pi]
    # set Lambda_0 0.33246029532469185632
    ## GRS67 eccentricity
    # set epsilon [expr sqrt(($a*$a-$b*$b)/($a*$a))]
    set epsilon 0.0818205680555
    # set A [expr pow(tan($pi/4.0+$Phi/2.0), $n)]
    set A [expr pow(tan(0.785398163397448+0.5*$Phi), $n)]
    # set B [expr pow((1.0-$epsilon*sin($Phi))/(1.0+$epsilon*sin($Phi)),
    #                 $n*$epsilon/2.0)]
    set sinPhi [expr sin($Phi)]
    set B [expr pow((1.0-$epsilon*$sinPhi)/(1.0+$epsilon*$sinPhi), \
		    0.5*$n*$epsilon)]
    # set phi [expr 2.0 * atan($k*$A*$B) - $pi/2.0]
    set phi [expr 2.0 * atan($k*$A*$B) - 1.570796326794897]
    # set lambda [expr $n*($Lambda-$Lambda_0)]
    set lambda [expr $n*($Lambda-0.33246029532469185632)]

    # Spherical phi, lambda to auxiliary phi, lambda
    ## EOV starting point
    # set phi_0 [expr (47.0 + 06.0/60.0 + 00.0000/3600.0) / 360.0 * 2.0*$pi]
    set phi_0 0.8220500776893292307101144
    # set phi_a [expr asin(sin($phi)*cos($phi_0)
    #                       - cos($phi)*sin($phi_0)*cos($lambda))]
    set cosphi [expr cos($phi)]
    set phi_a [expr asin(sin($phi)*cos($phi_0)- \
			     $cosphi*sin($phi_0)*cos($lambda))]
    # set lambda_a [expr asin((cos($phi)*sin($lambda))/(cos($phi_a)))]
    set lambda_a [expr asin(($cosphi*sin($lambda))/cos($phi_a))]

    # Auxiliary phi, lambda to cartesian x, y (x == northing  &&  y == easting)
    ## Radius of the "new" Gauss-sphere
    # set R 6379743.001
    ## Reduction (scale factor)
    # set m_0 0.99993
    set Rm_0 6379296.41898993
    # set x [expr $R * $m_0 * log(tan($pi/4+$phi_a/2))]
    # set y [expr $R * $m_0 * $lambda_a]
    set x [expr $Rm_0 * log(tan(0.78539816339744830961+0.5*$phi_a))]
    set y [expr $Rm_0 * $lambda_a]

    # MF change: the false easting and false northing are added by the
    #  grid procedures; return the "raw" easting and northing
    return [list $y $x]
}

proc ProjEOVInvert {data x y} {
    # return list with latitude and longitude for point at $x,$y in the terrain
    #  where $x is the easting and $y the northing (against the EOV convention)
    #  with no false easting/northing added
    #  $data is not used
    # return 0, 180 if computation of latitude diverges
    global MESS TXT

    # Return to the EOV convention ($y == easting  &&  $x == northing)
    set z $x ; set x $y ; set y $z

    # MF changes after this point: for efficiency, original code commented out
    # Cartesian y, x to auxiliary phi, lambda
    # set R 6379743.001
    ## Reduction (scale factor)
    # set m_0 0.99993
    set Rm_0 6379296.41898993
    ## Pi :-P
    # set pi 3.1415926535897932384626434
    # set phi_a [expr 2*atan(exp($x/($R*$m_0)))-$pi/2]
    set phi_a [expr 2*atan(exp($x/$Rm_0))-1.57079632679489661923]
    # set lambda_a [expr $y/($R*$m_0)]
    set lambda_a [expr $y/($Rm_0)]

    # Auxiliary phi, lambda to spherical phi, lambda
    ## EOV starting point
    # set phi_0 [expr (47.0 + 06.0/60.0 + 00.0000/3600.0) / 360.0 * 2.0*$pi]
    # set phi_0 0.82205007768932923029
    set cosphi_a [expr cos($phi_a)]
    # set phi [expr asin( sin($phi_a)*cos($phi_0) +
    #                     cos($phi_a)*sin($phi_0)*cos($lambda_a) )]
    set phi [expr asin( sin($phi_a)*0.680720868959 + \
			    $cosphi_a*0.732542898787*cos($lambda_a) )]
    # set lambda [expr asin( (cos($phi_a)*sin($lambda_a)) / cos($phi) )]
    set lambda [expr asin( ($cosphi_a*sin($lambda_a)) / cos($phi) )]
    
    # Spherical phi, lambda to HD72 Phi, Lambda
    ## GRS67 semi-major / semi-minor axis
    # set a 6378160.0
    # set b 6356774.516
    ## EOV parameters
    set k 1.003110007693
    set n 1.000719704936
    # set Lambda_0 [expr (19.0 + 02.0/60.0 + 54.8584/3600.0) / 360.0 * 2.0*$pi]
    # set Lambda_0 0.33246029532469185632
    ## GRS67 eccentricity
    # set epsilon [expr sqrt(($a*$a-$b*$b)/($a*$a))]
    set epsilon 0.0818205680555
    ## Iteration to calculate Phi
    set Phi $phi
    set i 0
    # set A [expr tan($pi/4+$phi/2)]
    set A [expr tan(0.785398163397448+0.5*$phi)]
    set invn [expr 1/$n]
    while 1 {
        set B [expr $k *pow( (1-$epsilon*sin($Phi)) / (1+$epsilon*sin($Phi)), \
				 $n*$epsilon/2)]
        # set Phi_new [expr 2 * atan(pow($A/$B, 1/$n)) - $pi/2]
	set Phi_new [expr 2 * atan(pow(1.0*$A/$B, $invn)) - 1.570796326794897]
	if { abs($Phi-$Phi_new) <= 1.5E-9 } {
            break
	}
	if { [incr i] > 20} {
            GMMessage [format $MESS(badinvproj) $TXT(PRJEOV)]
            return [list 0 180]
	}
	set Phi $Phi_new
    }
    ## Calculate Lambda
    # set Lambda [expr $Lambda_0 + $lambda/$n]
    set Lambda [expr 0.33246029532469185632 + $lambda/$n]

    # MF change: results will be for the HD72
    # Convert the result from radian to DD.DDDDDDDD...
    set la [expr $Phi*57.29577951308232087684]
    set lo [expr $Lambda*57.29577951308232087684]
    #---
    return [list $la $lo]
}




