#
#  gpsman --- GPS Manager: a manager for GPS receiver data
#
#  Copyright (c) 2007 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: projections.tcl
#  Last change:  16 April 2007
#

## some information and algorithms used here adopted from
#  - libproj and PROJ4.3.3 - cartographic projection software
#      by Gerald I. Evenden (gie@kai.er.usgs.gov)
#  - Guidance Note Number 7, European Petroleum Survey Group
#      Revised November 1999
#  - Formulas and constants for the calculation of the Swiss conformal
#      cylindrical projection and for the transformation between coordinate
#      systems (http://www.swisstopo.ch/pub/data/geo/refsyse.pdf),
#      September 2001
#  - Het stelsel van de Rijksdriehoeksmeting en het European Terrestrial
#    Reference System 1989 (containing formulas and constants for the
#    Schreiber projection and the RD grid, used in The Netherlands;
#    partial English translation kindly provided by Rob Buitenhuis,
#    rob@buitenhs.demon.nl), http://www.rdnap.nl
#  - "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)
#  - IGN site (http://www.ign.fr), August 2003: French projections and
#    grids
#  - The California Spatial Library (http://www.gis.ca.gov/albers.epl),
#    October 2003: Teale-Albers projection and grid

##
## some algorithms and formulae adopted from:
#     - GKK grid conversion by Andreas Lange (Andreas.C.Lange@GMX.de)
#       in GPStrans 0.36
#     - Swedish grid conversion by Anders Lennartsson
#       (anders.lennartsson@sto.foa.se) in GPStrans 0.36
##
## information kindly provided by
#     - Dan Jacobson (jidanni@yahoo.com.tw) on the Taiwan Grid projection
#     - Alessandro Palmas (alpalmas@tin.it) on the French and Austrian grids
##

## "projection" here means the conversion from geodetic position to
# planar Cartesian coordinates im meters in the terrain
#  these coordinates can then be converted to map coordinates by the
# transformations defined in maptransf.tcl
#
# projection procs (direct and inverse) must work with longitude values in
#  the -179,180 range although they should use internally the equivalent
#  value in the -359,360 range that is closer to the central or false
#  longitude (if there is one):
#
#  no changes in the map datum ($Datum) will occur after initialization
##

## when adding new projection procedures these variables must be updated
#  as well as PRJ$name and parameters (except reserved ones) in files lang*.tcl
# see also the comments preceding each procedure for the UTM projection
##

set MAPKNOWNPROJS {UTM TM BMN BNG BWI CMP CTR GKK ITM KKJY KKJP SEG TWG
                   LCC1 LCC2 IcG Merc1 Merc2 CS APOLY Stereogr Schreiber RDG
                   SOM LV03 Lamb93 LambNTFe LambNTF AEA LEAC TAlbers}

     # indices of parameters array used for each main projection
     #  index "datum" is also used, but is dealt with in a different way
     #  indices starting with "_" are used for other values whose meanings
     #  are defined in the corresponding proc Proj${proj}ComputeAux
array set MAPPROJDATA {
    UTM   {UTMzone m_0}
    TM    {lat0 long0 k0}
    BMN   {bmnlong0 m_0}
    CTR   {ctrlong0 m_0}
    GKK   {gkklong0 m_0}
    KKJP  {kkjplong0 m_0}
    TWG   {twglong0 m_0}
    LCC1  {lat0 long0 k0 m_e m_a lcc_n lcc_rho0F lcc_sn}
    LCC2  {lat1 lat2 latF longF m_e m_a lcc_n lcc_rho0F lcc_sn}
    LambNTF {NTFzone m_e m_a lcc_n lcc_rho0F lcc_sn}
    AEA   {lat1 lat2 lat0 long0 m_e m_a m_0 m_1 m_2 m_3 m_4 m_5}
    LEAC  {lat1 polasp lat0 long0 m_e m_a m_0 m_1 m_2 m_3 m_4 m_5}
    Merc1 {long0 k0 m_e m_a m_1 m_2 m_3 m_4}
    Merc2 {lat1 long0 m_e m_a m_1 m_2 m_3 m_4}
    CS    {lat0 long0 m_e m_a m_0 m_1 m_2 m_3 m_4 m_5 m_6 m_7 m_8}
    APOLY {lat0 m_a m_0 m_1 m_2 m_3 m_4 m_5 m_6}
    Stereogr   {lat0 long0 k0 m_a m_e m_0 m_1 m_2 m_3}
    Schreiber  {m_0 m_1 m_2 m_3 m_4 m_5}
    SOM   {lat0 long0 k0 m_e m_0 m_1 m_2 m_3 m_4 m_5 m_6}
}

     # new types here imply changing proc BadParam in check.tcl
     # not in use, but defined: lat=MIN,MAX  long=MIN,MAX
array set MAPPROJDTYPE {
    lat0      lat        lat1      lat         lat2      lat
    latF      lat
    long0     long       longF     long
    k0        float>0
    bmnlong0  list=10.333333333,13.333333333,16.333333333
    ctrlong0  list=9,15
    gkklong0  list=0,3,6,9,12,15
    kkjplong0 list=21,24,27,30
    twglong0  list=115,117,119,121,123,125
    UTMzone   zone=UTM
    NTFzone   zone=LambNTF
    polasp    list:north,south
    lcc_n     reserved    lcc_rho0F reserved    lcc_sn    reserved
    m_e       reserved    m_a       reserved
    m_0       reserved    m_1       reserved
    m_2       reserved    m_3       reserved
    m_4       reserved    m_5       reserved
    m_6       reserved    m_7       reserved
    m_8       reserved
}

     # projections needing auxiliary parameters computed from the main ones
     #  a specific proc is needed to compute the auxiliary parameters
     #  (see, e.g., proc ProjLCC1ComputeAux)
set MAPPROJAUX {UTM LCC1 LCC2 Merc1 Merc2 CS APOLY Stereogr Schreiber SOM
                AEA LEAC}
     # auxiliary parameters depend on datum?
array set MAPPROJDEPDATUM {
    UTM  0      LCC1 1   LCC2  1      Merc1 1    Merc2 1    CS    1   APOLY 1
    Stereogr 1  SOM  1   Schreiber 0  LambNTF 1  AEA   1    LEAC  1
}

     # main projections admitting particular cases
set MAINPROJS {TM LCC1 LCC2 Merc1 Merc2 CS APOLY Stereogr SOM AEA LEAC}

     # parameters for special cases of Transverse Mercator:
     #  prefix of long0 parameter name, minimum and maximum values for
     #  central longitude, longitude range in each zone, scale factor
     
array set TMPARAM {
    BMN  {bmn 10.333333333 16.333333333 3 1}
    CTR  {ctr 9 15 6 0.9996}
    GKK  {gkk 0 15 3 1}
    KKJP {kkjp 21 30 3 1}
    TWG  {twg 115 125 2 0.9999}
}

     # parameters for the zones of Lambert NTF projection (LCC2)
     #  longitude of false origin, list of zones aligned with pairs
     #  of extreme zone latitudes and longitudes, for each zone
     #  false easting and false northing

array set LNTFPARAMS {
    longF 2.33720833333333333333
    zones {I II III IV}
    extrlats {48.15 51.3 45.45 48.15 42.3 45.45 41 42.5}
    extrlongs {-6 10 -6 10 -6 10 8 10}
    fen,I {600000 200000}
    fen,II {600000 200000}
    fen,III {600000 200000}
    fen,IV {234.358 185861.369}
}

    # pre-computed parameters for each zone of Lambert NTF projection (LCC2)

array set LNTFzI {
    NTFzone I
    lat1  50.3959116667
    lat2  48.5985227778
    longF  2.33720833333
    latF  49.5
    datum  {NTF (Nouvelle Triangulation de France)}
    lcc_rho0F  5457616.6655
    lcc_sn  1
    m_a  11603796.7976
    m_e  0.0824834000442
    lcc_n  0.760405965438
}

array set LNTFzII {
    NTFzone II
    lat1  45.8989188889
    lat2  47.6960144444
    longF  2.33720833333
    latF  46.8
    datum  {NTF (Nouvelle Triangulation de France)}
    lcc_rho0F  5999695.76018
    lcc_sn  1
    m_a  11745793.2222
    m_e  0.0824834000442
    lcc_n  0.72896862668
}

array set LNTFzIII {
    NTFzone III
    lat1  43.1992913889
    lat2  44.9960938889
    longF  2.33720833333
    latF  44.1
    datum  {NTF (Nouvelle Triangulation de France)}
    lcc_rho0F  6591905.05462
    lcc_sn  1
    m_a  11947992.3451
    m_e  0.0824834000442
    lcc_n  0.695912797771
}

array set LNTFzIV {
    NTFzone IV
    lat1  41.5603877778
    lat2  42.0000593542
    longF  2.33720833333
    latF  42.165
    datum  {NTF (Nouvelle Triangulation de France)}
    lcc_rho0F  7106637.87866
    lcc_sn  1
    m_a  12178813.9766
    m_e  0.0824834000442
    lcc_n  0.666276778805
}

     # particular cases of main projections:
     #  CMP  - used in Portuguese Military Maps
     #  BNG  - used for the British National Grid
     #  BWI  - used for the British West Indies Grid
     #  IcG  - used for the Iceland Grid
     #  ITM  - used for the Irish TM Grid
     #  KKJY - used for the Uniform Finnish Grid
     #  Lamb93   - used in the new French geodetic system RGF93
     #  LambNTFe - used for the French grid NTF II étendue
     #  LV03 - used for the Swiss grid in the LV03 frame
     #  RDG  - used for the Netherlands Grid
     #  SEG  - used for the Swedish Grid
     #  TAlbers - used for the Teale Albers grid (California, USA)
array set MAPPARTPROJ {
    CMP  TM
    BNG  TM
    BWI  TM
    IcG  LCC2
    ITM  TM
    Lamb93   LCC2
    LambNTFe LCC2
    LV03 SOM
    KKJY TM
    RDG  Schreiber
    SEG  TM
    TAlbers AEA
}

     # parameters aligned with MAPPROJDATA above
array set MAPPARTPDATA {
    CMP  {39.66666666666667 -8.13190611111111 1}
    BNG  {49.0 -2.0 0.9996012717}
    BWI  {0 -62.0 0.9995}
    IcG  {64.75 64.25 65 -19}
    ITM  {53.5 -8.0 1.000035}
    KKJY {0 27.0 1}
    Lamb93  {44.0 49.0 46.5 3.0}
    LambNTFe {45.8989188888888888889 47.69601444444444444444
                 46.8 2.33720833333333333333}
    LV03 {46.9524055556 7.43958333333 1}
    RDG  {}
    SEG  {0 15.808277777778 1}
    TAlbers {34.0 40.5 0 -120.0}
}

     # projections and known grids
array set GRIDFOR {
    AEA  UTM/UPS
    APOLY UTM/UPS
    BMN  BMN
    BNG  BNG
    BWI  BWI
    CMP  CMP
    CS   UTM/UPS
    CTR  CTR
    GKK  GKK
    IcG  IcG
    ITM  ITM
    KKJY KKJY
    KKJP KKJP
    Lamb93   Lamb93
    LambNTF  LambNTF
    LambNTFe LambNTFe
    LCC1 UTM/UPS
    LCC2 UTM/UPS
    LEAC UTM/UPS
    LV03 LV03
    Merc1 UTM/UPS
    Merc2 UTM/UPS
    RDG  RDG
    SEG  SEG
    Schreiber RDG
    SOM  UTM/UPS
    Stereogr UTM/UPS
    TAlbers  TAlbers
    TM   UTM/UPS
    TWG  TWG
    UTM  UTM/UPS
}

     # grids with no zones
array set GRIDDEF {
    BWI {BWI 400000 0 KM {0 1e6 0 1e6} {0 10 -66 -55}}
    CMP {CMP 200000 300000 KM {64e3 368e3 0 580e3} {36.95 42.17 -9.66 -6.08}
         Lisboa}
    IcG {IcG 500000 500000 KM {0 1e6 0 1e6} {63.5 65 -30 -10} {WGS 84}}
    Lamb93 {Lamb93 700000 6600000 KM {0 1e10 0 1e10}
            {41 52 -6 11} {WGS 84}}
    LambNTFe {LambNTFe 600000 2200000 KM {0 1e10 0 1e10}
 	      {41 52 -6 11} {NTF (Nouvelle Triangulation de France)}}
    LV03 {LV03 600000 200000 KM {0 1e6 0 1e6} {45 48 6 11} CH-1903}
    RDG {RDG 155000 463000 KM {0 29e4 29e4 63e4} {50.3 53.45 3 7.45}
         {Rijks Driehoeksmeting}}
    SEG {SEG 1500000 0 KM {12e5 19e5 61e5 77e5} {54 70 10 26}}
    TAlbers {TAlbers 0 -4000000 KM {-1.5e6 1.5e6 -1.5e6 1.5e6}
	     {25 45 -130 -110} {NAD27 CONUS}}
}

     # zones of grids (regular expressions to recognize valid zones)
     #  UTM is added below after use of this array in colecting the grid
     #  names
array set GRIDZN {
    BMN  {^M((28)|(31)|(34))$}
    BNG  {^([A-H]|[J-Z])([A-H]|[J-Z])$}
    CTR  {^[1-2]$}
    GKK  {^[0-5]$}
    ITM  {^[A-H]|[J-Z]$}
    KKJY {^27E$}
    KKJP {^[1-4]$}
    LambNTF {^I(I|II|V)?$}
    TWG  {^[1-6]$}
}

     # position formats; $GRIDS may change dynamically
     # type for grids set below and in proc SetUpNZGrid
     #  x,y indices set only for numeric formats
array set POSTYPE {
    DDD  latlong
    DDD,xyixs  {1 0}
    DMS  latlong
    DMM  latlong
    GRA  latlong
    GRA,xyixs  {3 2}
    UTM/UPS  utm
    UTM/UPS,xyixs {4 5}
    MH   mh
}

     # non-numeric position formats
set NNUMPFORMATS {DMM DMS MH}
set PFORMATS ""
foreach f [set NONGRIDFMTS {DMS DMM DDD GRA UTM/UPS MH}] {
    lappend PFORMATS $TXT($f)
    # but UTM/UPS is set to 1 below
    set ZGRID($f) 0
}
set GRIDS [lsort [concat [array names GRIDZN] [array names GRIDDEF]]]
set PFORMATS [concat $PFORMATS $GRIDS]

     # must be here, after computing $GRIDS
     #  there is a redundancy here: UTM for the projection and UTM/UPS for
     #  the position format...
array set GRIDZN {
    UTM  {^([0-5]?[0-9]|60)[A-HJ-NP-Z]$}
    UTM/UPS  {^([0-5]?[0-9]|60)[A-HJ-NP-Z]$}
}

     # initialize projection and grid information
set ZGRID(UTM/UPS) 1
set TXT(PRJUTM/UPS) $TXT(UTM/UPS)

foreach g [array names TMPARAM] {
    # central longitude parameter
    set TXT([lindex $TMPARAM($g) 0]long0) $TXT(long0)
}

foreach g $GRIDS {
    # all grids named as they should appear in TXT
    set TXT($g) $g ; set INVTXT($g) $g
    # grid has more than 1 zone (see SetUpNZGrid for those that do not)
    set ZGRID($g) 1
    set POSTYPE($g) grid ; set POSTYPE($g,xyixs) {3 4}
    # at present, all grids with >1 zone with fixed datum must be given
    #  after this foreach
    set GRD${g}(datum) ""
}

    # all grids with >1 zone with fixed datum must be given here!
set GRDLambNTF(datum) "NTF (Nouvelle Triangulation de France)"

## grids

proc SetUpNZGrids {} {
    # initialize grids with no zones defined in array GRIDDEF
    global GRIDDEF

    foreach g [array names GRIDDEF] {
	SetUPNZGrid $g $GRIDDEF($g)
    }
    return
}

proc SetUPNZGrid {gr params} {
    # initialize grid with no zones
    #  $params is list of values as in array GRIDDEF
    global ZGRID GRIDZN MAPPROJAUX MAPPROJDEPDATUM MAPPARTPROJ MESS \
	    DSCALEFOR GRD$gr POSTYPE

    set ZGRID($gr) 0 ; set GRIDZN($gr) {}
    set POSTYPE($gr) nzgrid ; set POSTYPE($gr,xyixs) {2 3}
    if { [llength $params] < 6 } { BUG bad params for SetUPNZGrid }
    foreach "proj fe fn unit xybounds llbounds datum" $params {}
    if { $datum == "" } {
	set fixdatum 0
	if { [lsearch -exact $MAPPROJAUX $proj] != -1 && \
		$MAPPROJDEPDATUM($proj) } {
	    GMMessage [format $MESS(gridneedsdatum) $TXT($proj)]
	    set datum "WGS 84" ; set fixdatum 1
	}
    } else { set fixdatum 1 }
    ProjInitPart $proj GRD$gr $datum
    set GRD${gr}(_proj) $MAPPARTPROJ($proj)
    foreach f "fe fn" {
	set GRD${gr}(_$f) [set $f]
    }
    set GRD${gr}(_unit) $unit
    # scale factor to get metre from user unit
    set GRD${gr}(_uscale) \
	    [expr 1000.0*$DSCALEFOR(sub,$unit)/$DSCALEFOR($unit)]
    foreach f "xmin xmax ymin ymax" b $xybounds {
	set GRD${gr}(_$f) $b
    }
    foreach f "lamin lamax lomin lomax" b $llbounds {
	set GRD${gr}(_$f) $b
    }
    if { $fixdatum } { set GRD${gr}(datum) $datum }
    return
}

proc DegreesToNZGrid {grid latd longd datum} {
    # convert lat/long in signed degrees to coordinates in no-zone grid
    # return "-- 0" on error
    global GRD$grid

    if { $latd < [set GRD${grid}(_lamin)] || \
	    $latd > [set GRD${grid}(_lamax)] || \
	    $longd < [set GRD${grid}(_lomin)] || \
	    $longd > [set GRD${grid}(_lomax)] } {
	return "-- 0"
    }
    set proj [set GRD${grid}(_proj)]
    if { [set GRD${grid}(datum)] == "" } {
	set GRD${grid}(datum) $datum
	set nodatum 1
    } else { set nodatum 0 }
    set p [Proj${proj}Point GRD$grid $latd $longd $datum]
    if { $nodatum } { set GRD${grid}(datum) "" }
    set s [set GRD${grid}(_uscale)]
    set x [expr int(round(1.0*[lindex $p 0]/$s+[set GRD${grid}(_fe)]))]
    set y [expr int(round(1.0*[lindex $p 1]/$s+[set GRD${grid}(_fn)]))]
    if { $x < [set GRD${grid}(_xmin)] || \
	    $x > [set GRD${grid}(_xmax)] || \
	    $y < [set GRD${grid}(_ymin)] || \
	    $y > [set GRD${grid}(_ymax)] } {
	return "-- 0"
    }
    return [list $x $y]
}

proc GridToDegrees {grid zone x y datum} {
    # convert grid coords to lat/long in signed degrees
    # return 0 on error, -1 if $datum is not that of $grid,
    #  otherwise list with lat and long
    global ZGRID GRD$grid

    if { $ZGRID($grid) } { return [${grid}ToDegrees $zone $x $y $datum] }
    # no-zone grid
    if { $x < [set GRD${grid}(_xmin)] || \
	    $x > [set GRD${grid}(_xmax)] || \
	    $y < [set GRD${grid}(_ymin)] || \
	    $y > [set GRD${grid}(_ymax)] } {
	return 0
    }
    set s [set GRD${grid}(_uscale)]
    set x [expr $s*($x-[set GRD${grid}(_fe)])]
    set y [expr $s*($y-[set GRD${grid}(_fn)])]
    if { [set gdatum [set GRD${grid}(datum)]] == "" } {
	set GRD${grid}(datum) $datum
	set nodatum 1
    } elseif { $gdatum != $datum } {
	return -1
    } else { set nodatum 0 }
    set proj [set GRD${grid}(_proj)]
    set p [Proj${proj}Invert GRD$grid $x $y]
    if { $nodatum } { set GRD${grid}(datum) "" }
    set latd [lindex $p 0] ; set longd [lindex $p 1]
    if { $latd < [set GRD${grid}(_lamin)] || \
	    $latd > [set GRD${grid}(_lamax)] || \
	    $longd < [set GRD${grid}(_lomin)] || \
	    $longd > [set GRD${grid}(_lomax)] } {
	return 0
    }
    return $p
}

proc GridDatum {grid args} {
    # return datum needed by grid, or first element of $args if none
    global GRD$grid

    if { [set d [set GRD${grid}(datum)]] == "" } {
	return [lindex $args 0]
    }
    return $d
}

## grids having zones
# UTM grid (see also posncomp.tcl)

proc CompUTMOnZone {lat long longz} {
    # UTM/UPS position representation of position with lat/long signed degrees
    #  on zone containing meridian with longitude $longz (signed degrees)
    global UTMlat0 UTMk0 Datum

    if { $lat>=-80 && $lat<=84 } {
	set long0 [expr 6*int(floor($longz/6.0))+3]
	set cs [ConvToTM $lat $long $UTMlat0 $long0 $UTMk0 $Datum]
	set x [expr int(5e5+[lindex $cs 0])] ; set y [lindex $cs 1]
	if { $lat < 0 } { set y [expr int(1e7+$y)] }
	return [list $x $y]
    }
    return [ConvToUPS $lat $long $Datum]
}

# BMN, CTR, GKK, basic KKJ, and Taiwan grid coordinates

proc DegreesToBMN {lat long datum} {
    # convert from lat/long in signed degrees to Austrian BMN grid coords
    # zone codes: M28, M31, M34

    if { $long < 8.833333333 || $long > 17.833333333 || $lat < 0 } {
	return "-- 0 0"
    }
    set z [expr int(($long-8.833333333)/3)]
    set fe [expr 150000+300000*$z]
    set long0 [expr 10.833333333+3*$z]
    set cs [ConvToTM $lat $long 0 $long0 1 $datum]
    if { [set x [expr round([lindex $cs 0]+$fe)]] < 0 || \
	    [set y [expr round([lindex $cs 1])]] < 0 } {
	return "-- 0 0"
    }
    return [list M[expr $z*3+28] $x $y]
}

proc BMNToDegrees {zone x y datum} {
    # convert from Austrian BMN grid coords to lat/long in signed degrees
    global GRIDZN

    if { $x < 0 || $y < 0 || $x > 6e6 || $y > 1e7 || \
	    ! [regexp $GRIDZN(BMN) $zone] } {
	return 0
    }
    set z [expr ([string range $zone 1 2]-28)/3]
    set fe [expr 150000+300000*$z]
    set long0 [expr 10.833333333+3*$z]
    return [ConvFromTM [expr $x-$fe] $y 0 $long0 1 $datum]
}

proc DegreesToCTR {lat long datum} {
    # convert from lat/long in signed degrees to Italian CTR grid coords
    # zone codes: 1-2

    if { $long < 6 || $long > 18 || $lat < 0 } {
	return "-- 0 0"
    }
    if { [set z [expr int(($long-6)/6)+1]] == 1 } {
	set fe 150e4
    } else { set fe 252e4 }
    set long0 [expr 3+6*$z]
    set cs [ConvToTM $lat $long 0 $long0 0.9996 $datum]
    if { [set x [expr round([lindex $cs 0]+$fe)]] < 0 || \
	    [set y [expr round([lindex $cs 1])]] < 0 } {
	return "-- 0 0"
    }
    return [list $z $x $y]
}

proc CTRToDegrees {zone x y datum} {
    # convert from Italian CTR grid coords to lat/long in signed degrees
    global GRIDZN

    if { $x < 0 || $y < 0 || $x > 6e6 || $y > 1e7 || \
	    ! [regexp $GRIDZN(CTR) $zone] } {
	return 0
    }
    if { $z == 1 } {
	set fe 150e4
    } else { set fe 252e4 }
    set long0 [expr 3+6*$z]
    return [ConvFromTM [expr $x-$fe] $y 0 $long0 0.9996 $datum]
}

proc DegreesToGKK {lat long datum} {
    # convert from lat/long in signed degrees to German Krueger grid coords
    # zone codes: 0-5

    if { $long < -1.5 || $long > 16.5 || $lat < 0 } {
	return "-- 0 0"
    }
    set z [expr int(($long+1.5)/3)]
    set long0 [expr 3.0*$z]
    set cs [ConvToTM $lat $long 0 $long0 1.0 $datum]
    if { [set x [expr round([lindex $cs 0]+5e5+1e6*$z)]] < 0 || \
	    [set y [expr round([lindex $cs 1])]] < 0 } {
	return "-- 0 0"
    }
    return [list $z $x $y]
}

proc GKKToDegrees {zone x y datum} {
    # convert from German Krueger grid coords to lat/long in signed degrees
    global GRIDZN

    if { $x < 0 || $y < 0 || $x > 6e6 || $y > 1e7 || \
	    ! [regexp $GRIDZN(GKK) $zone] } {
	return 0
    }
    set long0 [expr 3.0*$zone]
    set x [expr $x-5e5-1e6*$zone]
    return [ConvFromTM $x $y 0 $long0 1.0 $datum]
}

proc DegreesToKKJP {lat long datum} {
    # convert from lat/long in signed degrees to basic Finnish grid coords
    # zone codes: 1-4

    if { $long < 19.5 || $long > 31.5 || $lat < 0 } {
	return "-- 0 0"
    }
    set z [expr int(($long-16.5)/3)]
    set long0 [expr 3.0*$z+18]
    set cs [ConvToTM $lat $long 0 $long0 1.0 $datum]
    if { [set x [expr round([lindex $cs 0]+5e5+1e6*$z)]] < 0 || \
	    [set y [expr round([lindex $cs 1])]] < 0 } {
	return "-- 0 0"
    }
    return [list $z $x $y]
}

proc KKJPToDegrees {zone x y datum} {
    # convert from basic Finnish grid coords to lat/long in signed degrees
    global GRIDZN

    if { $x < 0 || $y < 0 || $x > 5e6 || $y > 1e7 || \
	    ! [regexp $GRIDZN(KKJP) $zone] } {
	return 0
    }
    set long0 [expr 3.0*$zone+18]
    set x [expr $x-5e5-1e6*$zone]
    return [ConvFromTM $x $y 0 $long0 1.0 $datum]
}

proc DegreesToTWG {lat long datum} {
    # convert from lat/long in signed degrees to Taiwan grid coords
    # zone codes: 1-6

    if { $long < 114 || $long > 126 || $lat < 0 } {
	return "-- 0 0"
    }
    set z [expr int(($long-114)/2)+1]
    set long0 [expr 2.0*$z+113]
    set cs [ConvToTM $lat $long 0 $long0 0.9999 $datum]
    if { [set x [expr round([lindex $cs 0]+25e4)]] < 0 || \
	    [set y [expr round([lindex $cs 1])]] < 0 } {
	return "-- 0 0"
    }
    return [list $z $x $y]
}

proc TWGToDegrees {zone x y datum} {
    # convert from Taiwan grid coords to lat/long in signed degrees
    global GRIDZN

    if { $x < 0 || $y < 0 || $x > 5e6 || $y > 1e7 || \
	    ! [regexp $GRIDZN(TWG) $zone] } {
	return 0
    }
    set long0 [expr 2.0*$zone+113]
    set x [expr $x-25e4]
    return [ConvFromTM $x $y 0 $long0 0.9999 $datum]
}

# Lambert NTF grid coordinates

proc DegreesToLambNTF {lat long datum} {
    # convert from lat/long in signed degrees to French Lambert NTF grid coords
    # zone codes: I, II, III, IV
    global LNTFPARAMS

    set projdatum "NTF (Nouvelle Triangulation de France)"
    if { $datum !=  $projdatum } {
	foreach "lat long" \
	    [ToDatum $lat $long $datum $projdatum] { break }
    }
    set zone "--"
    foreach z $LNTFPARAMS(zones) "lamin lamax" $LNTFPARAMS(extrlats) \
	    "lomin lomax" $LNTFPARAMS(extrlongs) {
	if { $lamin <= $lat && $lat <= $lamax && $lomin <= $long && \
		 $long <= $lomax } {
	    set zone $z ; break
	}
    }
    if { $zone == "--" } { return "-- 0 0" }
    foreach "x y" [ProjLCC2Point LNTFz$zone $lat $long $projdatum] {}
    foreach "fe fn" $LNTFPARAMS(fen,$zone) {}
    return [list $zone [expr round($x+$fe)] [expr round($y+$fn)]]
}

proc LambNTFToDegrees {zone x y datum} {
    # convert from French NTF grid coords to lat/long in signed degrees
    global GRIDZN LNTFPARAMS

    if { $x < 0 || $y < 0 || $x > 1e6 || $y > 1e6 || \
	     ! [regexp $GRIDZN(LambNTF) $zone] } {
	return 0
    }
    foreach "fe fn" $LNTFPARAMS(fen,$zone) {}
    set p [ProjLCC2Invert LNTFz$zone [expr $x-$fe] [expr $y-$fn]]
    set projdatum "NTF (Nouvelle Triangulation de France)"
    if { $datum !=  $projdatum } {
	set p [ToDatum [lindex $p 0] [lindex $p 1] $projdatum $datum]
    }
    return $p
}

## Maidenhead locators

set MHCHARS ABCDEFGHIJKLMNOPQRSTUVWX

proc DegreesToMHLoc {latd longd} {
    # convert from lat/long in signed degrees to Maidenhead locator
    #  (6 characters)
    global MHCHARS

    set longd [expr $longd+180] ; set latd [expr $latd+90]
    set k [expr int(floor($longd/20))] ; set c1 [string index $MHCHARS $k]
    set longd [expr $longd-$k*20]
    set k [expr int(floor($latd/10))]
    set c2 [string index $MHCHARS $k]
    set latd [expr $latd-$k*10]
    set c3 [expr int(floor($longd/2))]
    set longd [expr $longd-$c3*2]
    set c4 [expr int(floor($latd))]
    set latd [expr $latd-$c4]
    set c5 [string index $MHCHARS [expr int(floor(($longd*60)/5))]]
    set c6 [string index $MHCHARS [expr int(floor(($latd*60)/2.5))]]
    return ${c1}${c2}${c3}${c4}${c5}${c6}
}

proc MHLocToDegrees {mhl} {
    # convert from Maidenhead locator (6 characters) to degrees
    global MHCHARS

    foreach i "0 1 2 3 4 5" {
	set c$i [string index $mhl $i]
    }
    set longd [expr 20*[string first $c0 $MHCHARS]+2*$c2+ \
	    5*[string first $c4 $MHCHARS]/60.0-180]
    set latd [expr 10*[string first $c1 $MHCHARS]+$c3+ \
	    2.5*[string first $c5 $MHCHARS]/60.0-90]
    return [list $latd $longd]
}

## user-defined projections and grids

proc DefineProjection {} {
    # set-up a user-defined case of a main projection
    global MAINPROJS MAPKNOWNPROJS MAPPROJDATA MAPPROJDTYPE MAPPARTPDATA \
	    MAPPARTPROJ GRIDFOR UPName UPAbbrev UPData UProjs TXT MESS

    catch {unset UPData} ; catch {unset UPName} ; catch {unset UPAbbrev}
    set nps ""
    foreach mp $MAINPROJS {
	lappend nps $TXT(PRJ$mp)
    }
    while 1 {
	set mp [GMChooseFrom single $TXT(baseproj) 20 $nps $MAINPROJS \
		"UPName UPAbbrev" [list "=$TXT(name)" "=$TXT(abbrev)"]]
	if { $mp == "" } { return }
	if { [ProjNamesOk] } { break }
    }
    foreach e $MAPPROJDATA($mp) {
	if { $MAPPROJDTYPE($e) != "reserved" } {
	    set UPData($e) ""
	}
    }
    set UPData(_grid) $GRIDFOR($mp)
    set UPData(_abbr) $UPAbbrev
    set UPData(datum) ""
    if { ! [ProjParams define $mp UPData] } { return }
    lappend UProjs $UPAbbrev ; lappend MAPKNOWNPROJS $UPAbbrev
    set MAPPARTPROJ($UPAbbrev) $mp
    set TXT(PRJ$UPAbbrev) $UPName
    set vals ""
    foreach e $MAPPROJDATA($mp) {
	if { $MAPPROJDTYPE($e) != "reserved" } {
	    lappend vals $UPData($e)
	}
    }
    set MAPPARTPDATA($UPAbbrev) $vals
    if { [set GRIDFOR($UPAbbrev) $UPData(_grid)] == $UPAbbrev } {
	SetUPNZGrid $UPAbbrev $UPData(_gparams)
	AddUserGrid $UPAbbrev
    }
    SaveUserProjsGrids
    return
}

proc ProjNamesOk {} {
    # check that projection name and abbreviation in $UPName and $UPAbbrev
    #  are acceptable
    global UPName UPAbbrev MAPKNOWNPROJS TXT MESS

    if { $UPName == "" || $UPAbbrev == "" } {
	GMMessage $MESS(projnameabbr) ; return 0
    } elseif { [regexp {[ ]} $UPAbbrev] } {
	GMMessage $MESS(abbrevhasspaces) ; return 0
    } elseif { [lsearch -exact $MAPKNOWNPROJS $UPAbbrev] != -1 } {
	GMMessage $MESS(abbrevinuse) ; return 0
    } else {
	foreach p $MAPKNOWNPROJS {
	    if { $TXT(PRJ$p) == $UPName } {
		GMMessage $MESS(nameinuse) ; return 0
	    }
	}
    }
    return 1
}

proc OpenUserProjection {proj} {
    # edit user-defined projection
    global MAPKNOWNPROJS MAPPROJDATA MAPPROJDTYPE MAPPARTPROJ MAPPARTPDATA \
	    GRIDFOR ZGRID MapProjection UPData UProjs TXT MESS MapPFormat \
	    UPName UPAbbrev

    catch {unset UPData}
    set mp $MAPPARTPROJ($proj)
    foreach e $MAPPROJDATA($mp) v $MAPPARTPDATA($proj) {
	if { $MAPPROJDTYPE($e) != "reserved" } {
	    set UPData($e) $v
	}
    }
    set gr $GRIDFOR($proj)
    if { $ZGRID($gr) } {
	set UPData(datum) ""
    } else {
	global GRD$gr
	set UPData(datum) [set GRD${gr}(datum)]
    }
    set UPData(_grid) $gr
    set UPData(_abbr) $proj
    if { $MapProjection == $proj } {
	GMMessage $MESS(projinuse)
	ProjParams edit $mp UPData
	return
    }
    set grinuse 0
    if { $gr == $proj } {
	set selfgrid 1 ; set lp ""
	foreach p [array names GRIDFOR] {
	    if { $GRIDFOR($p) == $proj && $p != $proj } {
		lappend lp $p ; set grinuse 1
	    }
	}
	if { $MapPFormat == $gr } {
	    lappend lp $TXT(nameMap) ; set grinuse 1
	}
	if { $grinuse && ! \
		[GMConfirm [format $MESS(gridinuse) $lp]] } { return }
    } else { set selfgrid 0 }
    if { [set op [ProjParams edit $mp UPData]] == 0 } {
	return
    }
    if { $op == 2 } {
	# create new projection
	set proj $UPAbbrev
	set projname $UPName
	set grinuse 0
    }
    set nogrid [string compare $UPData(_grid) $proj]
    if { $grinuse && ($op == -1 || $nogrid) } {
	GMMessage [format $MESS(gridinusenochg) $lp]
	return
    }
    switch -- $op {
	1 {
	    # save new parameters
	    set vals ""
	    foreach e $MAPPROJDATA($mp) {
		if { $MAPPROJDTYPE($e) != "reserved" } {
		    lappend vals $UPData($e)
		}
	    }
	    set MAPPARTPDATA($proj) $vals
	    set GRIDFOR($proj) $UPData(_grid)
	    if { ! $nogrid } {
		SetUPNZGrid $proj $UPData(_gparams)
		if { ! $selfgrid } {
		    AddUserGrid $proj
		}
	    } elseif { $selfgrid } {
		ForgetUserGrid $proj
	    }
	}
	2 {
	    # create new
	    lappend UProjs $UPAbbrev ; lappend MAPKNOWNPROJS $UPAbbrev
	    set MAPPARTPROJ($UPAbbrev) $mp
	    set TXT(PRJ$UPAbbrev) $UPName
	    set vals ""
	    foreach e $MAPPROJDATA($mp) {
		if { $MAPPROJDTYPE($e) != "reserved" } {
		    lappend vals $UPData($e)
		}
	    }
	    set MAPPARTPDATA($UPAbbrev) $vals
	    if { [set GRIDFOR($UPAbbrev) $UPData(_grid)] == $UPAbbrev } {
		SetUPNZGrid $UPAbbrev $UPData(_gparams)
		AddUserGrid $UPAbbrev
	    }
	    SaveUserProjsGrids
	}
	-1 {
	    # forget projection
	    set UProjs [Delete $UProjs $proj]
	    set MAPKNOWNPROJS [Delete $MAPKNOWNPROJS $proj]
	    if { $selfgrid } {
		ForgetUserGrid $proj
	    }
	    unset MAPPARTPROJ($proj)
	    unset TXT(PRJ$proj)
	    unset MAPPARTPDATA($proj)
	    unset GRIDFOR($proj)
	}
    }
    SaveUserProjsGrids
    return
}

proc AddUserGrid {grid} {
    # add user grid to list of known grids
    global GRIDS PFORMATS TXT INVTXT POSTYPE

    set TXT($grid) $grid ; set INVTXT($grid) $grid
    set GRIDS [lsort [linsert $GRIDS 0 $grid]]
    set PFORMATS [concat [list $TXT(DMS) $TXT(DMM) $TXT(DDD) $TXT(UTM/UPS)] \
	    $GRIDS]
    set POSTYPE($grid) nzgrid ; set POSTYPE($grid,xyixs) {2 3}
    return
}

proc ForgetUserGrid {grid} {
    # delet user grid from list of known grids
    global GRIDS PFORMATS TXT INVTXT GRD$grid

    unset TXT($grid) ; unset INVTXT($grid)
    set GRIDS [Delete $GRIDS $grid]
    set PFORMATS [Delete $PFORMATS $grid]
    unset GRD$grid
    return
}

proc SaveUserProjsGrids {} {
    # save user-defined projections and grids
    global UFile UProjs MAPPARTPROJ MAPPARTPDATA GRIDFOR ZGRID TXT INVTXT \
	    MESS

    if { [catch {set f [open $UFile(proj,grid) w]}] } {
	GMMessage $MESS(cantwrtprgr)
	return
    }
    puts $f ""
    puts $f "# $MESS(written) GPSMan [NowTZ]"
    puts $f "# $MESS(editrisk)"
    puts $f ""
    puts $f "set UProjs \{$UProjs\}"
    puts $f ""
    foreach p $UProjs {
	puts $f "set MAPPARTPROJ($p) $MAPPARTPROJ($p)"
	puts $f "set MAPPARTPDATA($p) \{$MAPPARTPDATA($p)\}"
	puts $f "set GRIDFOR($p) $GRIDFOR($p)"
	puts $f "set TXT(PRJ$p) \"$TXT(PRJ$p)\""
	if { $GRIDFOR($p) == $p } {
	    global GRD$p

	    puts $f "global GRD$p"
	    puts $f "array set GRD$p \{"
	    foreach a [array names GRD$p] {
		puts $f "    $a  \"[set GRD${p}($a)]\""
	    }
	    puts $f "\}"
	    puts $f "AddUserGrid $p"
	    puts $f "set ZGRID($p) 0"
	}
	puts $f ""
    }
    if { $UProjs != "" } {
	puts $f "set MAPKNOWNPROJS \[concat \$MAPKNOWNPROJS \$UProjs\]"
	puts $f ""
    }
    close $f
    return
}

proc SetupUserProjs {} {
    # process information on user-defined projections and grids
    global UFile UProjs MAPKNOWNPROJS MAPPARTPROJ MAPPARTPDATA GRIDFOR \
	    ZGRID TXT INVTXT MESS

    source $UFile(proj,grid)
    # check for multiple entries in $MAPKNOWNPROJS
    set l $MAPKNOWNPROJS ; set MAPKNOWNPROJS ""
    while { $l != "" } {
	set p [lindex $l 0] ; set l [lreplace $l 0 0]
	if { [set i [lsearch -exact $l $p]] != -1 } {
	    GMMessage [format $MESS(redefproj) $TXT(PRJ$p)]
	    set l [lreplace $l $i $i]
	    while { [set i [lsearch -exact $l $p]] != -1 } {
		set l [lreplace $l $i $i]
	    }
	}
	lappend MAPKNOWNPROJS $p
    }
    return
}

## changing, defining and editing parameters

proc ProjParams {how proj data} {
    # open dialog for defining or confirming the projection parameters
    #  $how in {set, change, define, edit}
    #  $proj in $MAPKNOWNPROJS is
    #    set, change: projection whose parameters must be set or may be changed
    #    define, edit: main projection, a particular case of which is
    #     being defined or edited
    #  $data is name of global array with projection parameters
    # if $how==change may turn ASKPROJPARAMS to false in current session
    # return 1 if there were changes or the parameters were defined;
    #        0 if there were no changes or the operation was cancelled;
    #        -1 if the projection is to be forgotten, when editing;
    #        2 if a new projection is to be created, when editing, with
    #          name and abbreviation given by $UPName and $UPAbbrev
    global $data MAPPROJDATA MAPPROJDTYPE GRIDS TXT MESS COLOUR ASKPROJPARAMS \
	    EPOSX EPOSY ProjRes Param DLUNIT DISTUNIT Datum MAPPROJAUX \
	    MAPPROJDEPDATUM MAPPARTPROJ UPName UPAbbrev

    catch {unset Param}
    toplevel .proj
    wm protocol .proj WM_DELETE_WINDOW { set ProjRes ok }
    wm title .proj "$TXT(projection)/GPS Manager"
    wm transient .proj
    wm geometry .proj +$EPOSX+$EPOSY
    bind .proj <Key-Return> { set ProjRes ok }

    frame .proj.fr -relief flat -borderwidth 5 -bg $COLOUR(selbg)
    label .proj.fr.title -text "???" -relief sunken
    label .proj.fr.text -text $MESS(projchg)
    set frs .proj.fr.frsel
    frame $frs -relief flat -borderwidth 0
    set n 0
    foreach p $MAPPROJDATA($proj) {
	set rs 1
	switch -glob $MAPPROJDTYPE($p) {
	    lat -  long -  lat=* -
	    long=* {
		regsub {=.+} $MAPPROJDTYPE($p) "" cw
		set Param($p) DDD
		set fr $frs.f$p
		frame $fr -relief flat -borderwidth 0
		menubutton $fr.pfmt -text $TXT(DDD) -relief raised \
			-width 6 -direction below -menu $fr.pfmt.m
		menu $fr.pfmt.m -tearoff 0
		foreach f "DMS DMM DDD" {
		    $fr.pfmt.m add command -label $TXT($f) \
			    -command "ParamPFormt $fr $cw $f $p"
		}
		entry $fr.val -width 12
		$fr.val insert 0 [set ${data}($p)]
		grid configure $fr.pfmt -column 0 -row 0 -sticky nesw
		grid configure $fr.val -column 1 -row 0 -sticky nesw
		grid configure $fr -column 1 -row $n -sticky nesw
	    }
	    list=* -  list:* {
		set val [set ${data}($p)]
		set fr $frs.f$p
		frame $fr -relief flat -borderwidth 0
		set type $MAPPROJDTYPE($p)
		set l [split [string range $type 5 end] ,]
		if { [string index $type 4] == ":" } {
		    set lt ""
		    foreach v $l {
			lappend lt $TXT($v)
		    }
		} else { set lt $l }
		set r 0 ; set c 0 ; set vc 0 ; set bi ""
		foreach v $l tv $lt {
		    radiobutton $fr.r$vc -text $tv -variable Param($p) \
			    -value $v -anchor w -selectcolor $COLOUR(check)
		    grid configure $fr.r$vc -column $c -row $r -sticky nesw
		    if { $val == $v } { set bi $vc }
		    if { [incr c] > 2 } { set c 0 ; incr r ; incr rs }
		    incr vc
		}
		if { $bi != "" } { $fr.r$bi invoke }
		grid configure $fr -column 1 -row $n -sticky nesw
	    }
	    reserved {
		continue
	    }
	    default {
		entry $frs.e$p -width 12
		$frs.e$p insert 0 [set ${data}($p)]
		grid configure $frs.e$p -column 1 -row $n \
			-sticky nesw
	    }
	}
	label $frs.l$p -text $TXT($p)
	grid configure $frs.l$p -column 0 -row $n -rowspan $rs \
		-sticky w
	incr n
    }
    if { $n == 0 } {
	destroy .proj
	return 0
    }
    if { $how != "change" && $how != "set" } {
	set grid 1
	set oldgr [set ${data}(_grid)]
	set abbr [set ${data}(_abbr)]
	if { $how == "edit" && $oldgr == $abbr } {
	    set selfgrid 1 ; set t UTM/UPS
	} else { set selfgrid 0 ; set t $oldgr }
	label $frs.lgr -text "$TXT(grid):"
	menubutton $frs.m_guse -text $TXT($t) -relief raised \
		-direction below -menu $frs.m_guse.m
	radiobutton $frs.r_guse -text $TXT(use) -variable Param(_guse) \
		-value 1 -anchor w -selectcolor $COLOUR(check) \
		-command "ParamGrid $abbr $data"
	menu $frs.m_guse.m -tearoff 0
	foreach g [linsert $GRIDS 0 UTM/UPS] {
	    if { $g != $abbr } {
		$frs.m_guse.m add command -label $TXT($g) \
			-command "set ${data}(_grid) {$g} ;\
			          $frs.m_guse configure -text {$TXT($g)}"
	    }
	}
	radiobutton $frs.r_gdef -text $TXT(create) -variable Param(_guse) \
		-value 0 -anchor w -selectcolor $COLOUR(check) \
		-command "ParamGrid $abbr $data"
	label $frs.l_gdef -text [set ${data}(_abbr)] -width 8
	frame $frs.zg -relief flat -borderwidth 0
	label $frs.zg.lu -text $TXT(unit)
	set Param(_gunit) $DISTUNIT
	menubutton $frs.zg.mu -text $DLUNIT($DISTUNIT,subdist) -relief raised \
		-direction below -menu $frs.zg.mu.m
	menu $frs.zg.mu.m -tearoff 0
	foreach u "KM NAUTMILE STATMILE" {
	    set nm $DLUNIT($u,subdist)
	    $frs.zg.mu.m add command -label $nm \
		    -command "$frs.zg.mu configure -text {$nm} ; \
		              set Param(_gunit) $u"
	}
	foreach k "fe fn" t "easting northing" {
	    label $frs.zg.l$k -text $TXT(f$t)
	    entry $frs.zg.e$k -width 12
	}
	frame $frs.zg.zb -relief flat -borderwidth 0
	label $frs.zg.zb.lbd -text "$TXT(bounds):"
	label $frs.zg.zb.lbdn -text $TXT(min)
	label $frs.zg.zb.lbdx -text $TXT(max)
	set zbn 0
	grid configure $frs.zg.zb.lbd -column 0 -row $zbn -sticky w
	grid configure $frs.zg.zb.lbdn -column 1 -row $zbn -sticky snew
	grid configure $frs.zg.zb.lbdx -column 2 -row $zbn -sticky snew
	incr zbn
	foreach k "x y la lo" t "easting northing lat long" {
	    label $frs.zg.zb.l$k -text $TXT($t)
	    entry $frs.zg.zb.e${k}n -width 12
	    entry $frs.zg.zb.e${k}x -width 12
	    grid configure $frs.zg.zb.l$k -column 0 -row $zbn -sticky w
	    grid configure $frs.zg.zb.e${k}n -column 1 -row $zbn -sticky snew
	    grid configure $frs.zg.zb.e${k}x -column 2 -row $zbn -sticky snew
	    incr zbn
	}
	if { [set datum [set ${data}(datum)]] == "" } {
	    set datum $Datum ; set nodatum 1
	} else { set nodatum 0 }
	menubutton $frs.zg.datum -text $datum -relief raised \
		-direction below -menu $frs.zg.datum.m
	menu $frs.zg.datum.m -tearoff 0
	FillDatumMenu $frs.zg.datum.m ParamGNewDatum

	if { [lsearch -exact $MAPPROJAUX $proj] != -1 && \
		$MAPPROJDEPDATUM($proj) } {
	    label $frs.zg.cfd -text $TXT(fixeddatum)
	    set Param(_gfdatum) 1
	} else {
	    checkbutton $frs.zg.cfd -text $TXT(fixeddatum) \
		    -variable Param(_gfdatum) -onvalue 1 -offvalue 0 \
		    -selectcolor $COLOUR(check) -command ParamGDatum
	    if { $nodatum } {
		$frs.zg.cfd deselect
	    } else { $frs.zg.cfd invoke }
	}

	if { $selfgrid } {
	    ParamGridSet $oldgr
	} else {
	    # defaults: use pre-defined grid
	    $frs.r_guse invoke
	}

	set zgn 0
	grid configure $frs.zg.lu -column 0 -row $zgn -sticky w
	grid configure $frs.zg.mu -column 1 -row $zgn -sticky snew
	incr zgn
	grid configure $frs.zg.lfe -column 0 -row $zgn -sticky w
	grid configure $frs.zg.efe -column 1 -row $zgn -sticky snew
	incr zgn
	grid configure $frs.zg.lfn -column 0 -row $zgn -sticky w
	grid configure $frs.zg.efn -column 1 -row $zgn -sticky snew
	incr zgn
	grid configure $frs.zg.zb -column 0 -row $zgn -columnspan 2 -sticky e \
		-ipadx 10
	incr zgn
	grid configure $frs.zg.cfd -column 0 -row $zgn -sticky w
	grid configure $frs.zg.datum -column 1 -row $zgn -sticky snew

	grid configure $frs.lgr -column 0 -row $n -columnspan 2 \
		-sticky w -ipady 4
	incr n
	grid configure $frs.r_guse -column 0 -row $n -sticky w
	grid configure $frs.m_guse -column 1 -row $n -sticky snew
	incr n
	grid configure $frs.r_gdef -column 0 -row $n -sticky w
	grid configure $frs.l_gdef -column 1 -row $n -sticky w
	incr n
	grid configure $frs.zg -column 0 -row $n -columnspan 2 -sticky e \
		-ipadx 10
    } else {
	set grid 0 ; set oldgr "" ; set selfgrid 0
    }
    frame .proj.fr.frbs -relief flat -borderwidth 0
    # text of next button is reconfigured below
    button .proj.fr.frbs.ok -text $TXT(ok) \
	    -command { .proj.fr.frbs.ok configure -state normal
		       set ProjRes ok }
    button .proj.fr.frbs.rv -text $TXT(revert) \
	    -command ".proj.fr.frbs.rv configure -state normal ; \
	              ProjParamsRevert $proj $data $grid \"$oldgr\" $selfgrid"
    button .proj.fr.frbs.cnc -text $TXT(cancel) -command {set ProjRes cancel}
    switch $how {
	set {
	    pack .proj.fr.frbs.ok .proj.fr.frbs.rv .proj.fr.frbs.cnc \
		    -side left -padx 5
	    pack .proj.fr.title .proj.fr.text $frs .proj.fr.frbs \
		    -side top -pady 5
	}
	change {
	    checkbutton .proj.fr.stop -text $TXT(dontaskagain) \
		    -variable ASKPROJPARAMS -onvalue 0 -offvalue 1 \
		    -selectcolor $COLOUR(check)
	    pack .proj.fr.frbs.ok .proj.fr.frbs.rv -side left -padx 5
	    pack .proj.fr.title .proj.fr.text $frs .proj.fr.stop \
		    .proj.fr.frbs -side top -pady 5
	}
	define {
	    .proj.fr.frbs.ok configure -text $TXT(create)
	    pack .proj.fr.frbs.ok .proj.fr.frbs.cnc -side left -padx 5
	    pack .proj.fr.title .proj.fr.text $frs .proj.fr.frbs \
		    -side top -pady 5
	}
	edit {
	    .proj.fr.frbs.ok configure -text $TXT(change)
	    button .proj.fr.frbs.fgt -text $TXT(forget) \
		    -command ".proj.fr.frbs.fgt configure -state normal ; \
		              set ProjRes forget"
	    button .proj.fr.frbs.crt -text $TXT(create) \
		    -command ".proj.fr.frbs.crt configure -state normal ; \
		              set ProjRes new"
	    pack .proj.fr.frbs.ok .proj.fr.frbs.rv .proj.fr.frbs.crt \
		     .proj.fr.frbs.fgt .proj.fr.frbs.cnc -side left -padx 2
	    pack .proj.fr.title .proj.fr.text $frs .proj.fr.frbs \
		    -side top -pady 5
	}
    }
    pack .proj.fr -side top
    update idletasks
    # cannot use RaiseWindow because of menus
    set gs [grab current]
    grab .proj
    while 1 {
	tkwait variable ProjRes
	switch $ProjRes {
	    cancel {
		set r 0 ; break
	    }
	    forget {
		if { ! [GMConfirm \
			[format $MESS(askforget) $TXT(projection)]] } {
		    set ProjRes 1
		    continue
		}
		set r -1 ; break
	    }
	}
	set l "" ; set newdata same
	foreach p $MAPPROJDATA($proj) {
	    switch -glob [set type $MAPPROJDTYPE($p)] {
		lat -  long -  lat=* -
		long=* {
		    regsub {=.+} $type "" cw
		    set nv [ParamGetCoord $p $cw [$frs.f$p.val get]]
		    if { $nv == "nil" } { set newdata error ; break }
		    set nv [lindex $nv 0]
		    set pw $frs.f$p.val
		}
		list=* -  list:* {
		    # there can be no wrong values in list; no need for $pw
		    set nv $Param($p)
		}
		reserved { continue }
		default {
		    set nv [$frs.e$p get]
		    set pw $frs.e$p
		}
	    }
	    if { $nv != [set ${data}($p)] } {
		if { [BadParam $TXT($p) $type $nv] } {
		    focus $pw
		    set newdata error
		    break
		}
		lappend l $p $nv
		set newdata diff
	    }
	}
	if { $newdata != "error" } {
	    if { $ProjRes == "new" } {
		while 1 {
		    if { [GMChooseParams $TXT(projection) "UPName UPAbbrev" \
			    [list "=$TXT(name)" "=$TXT(abbrev)"]] } {
			if { [ProjNamesOk] } { set canc 0 ; break }
		    } else {
			set canc 1 ; break
		    }
		}
		if { $canc } { continue }
		set abbr $UPAbbrev
	    }
	    if { $newdata == "diff" } {
		array set $data $l
	    }
	    if { $grid } {
		global GRD$oldgr

		set newgrid same
		if { $Param(_guse) } {
		    if { $oldgr != [set ${data}(_grid)] } {
			set newgrid diff
		    }
		} else {
		    # check parameters for new/self grid
		    set unit $Param(_gunit)
		    set ${data}(_grid) $abbr
		    if { $selfgrid && \
			     $unit != [set GRD${oldgr}(_unit)] } {
			set sameunit 0
		    } else { set sameunit 1 }
		    set diff 0
		    foreach e "efe efn zb.exn zb.exx zb.eyn zb.eyx \
			    zb.elan zb.elax zb.elon zb.elox" \
			    t "float float float float float float \
			    lat lat long long" \
			    v "fe fn xmin xmax ymin ymax lamin lamax \
			    lomin lomax" \
			    n "feasting fnorthing bounds bounds bounds bounds \
			    bounds bounds bounds bounds" {
			set $v [$frs.zg.$e get]
			if { $selfgrid && $sameunit && \
				 [set $v] == [set GRD${oldgr}(_$v)] } {
			    continue
			}
			if { [BadParam $TXT($n) $t [set $v]] } {
			    focus $frs.zg.$e
			    set newgrid error
			    break
			}
			set diff 1
		    }
		    if { $newgrid == "error" } { continue }
		    if { $diff } {
			if { $xmin >= $xmax || $ymin >= $ymax || \
				$lamin >= $lamax || $lomin >= $lomax } {
			    GMMessage [format $MESS(badparam) $TXT(bounds)]
			    continue
			}
			set newgrid diff
		    }
		    if { $Param(_gfdatum) } {
			set datum [$frs.zg.datum cget -text]
			set newgrid diff
		    } elseif { ! $nodatum } {
			set datum ""
			set newgrid diff
		    } else { set datum "" }
		}
		if { $newgrid == "diff" } {
		    set newdata diff
		} else {
		    # make sure grid parameters are set
		    foreach v "fe fn unit xmin xmax ymin ymax \
                               lamin lamax lomin lomax" {
			set $v [set GRD${oldgr}(_$v)]
		    }
		}
		set ${data}(_gparams) \
		    [list $abbr $fe $fn $unit \
			 [list $xmin $xmax $ymin $ymax] \
			 [list $lamin $lamax $lomin $lomax] \
			 $datum]
	    }
	    if { $ProjRes == "new" } {
		set r 2
	    } elseif { $newdata == "diff" } {
		set r 1
	    } else { set r 0 }
	    break
	}
    }
    DestroyRGrabs .proj $gs
    update idletasks
    return $r
}

proc ParamGridSet {gr} {
    # show parameters for grid $gr
    global GRD$gr Param DLUNIT

    set frs .proj.fr.frsel
    $frs.r_gdef invoke
    set Param(_gunit) [set GRD${gr}(_unit)]
    $frs.zg.mu configure -text $DLUNIT($Param(_gunit),subdist)
    foreach k "fe fn" {
	$frs.zg.e$k delete 0 end
	$frs.zg.e$k insert 0 [set GRD${gr}(_$k)]
    }
    foreach k "x y la lo" {
	foreach kk "n x" p "min max" {
	    $frs.zg.zb.e${k}$kk delete 0 end
	    $frs.zg.zb.e${k}$kk insert 0 [set GRD${gr}(_${k}$p)]
	}
    }
    return
}

proc ParamGDatum {} {
    # fixed datum for grid has been de-/selected
    global Param

    if { $Param(_gfdatum) } {
	.proj.fr.frsel.zg.datum configure -state normal
    } else {
	.proj.fr.frsel.zg.datum configure -state disabled
    }
    return
}

proc ParamGNewDatum {datum args} {
    # new datum has been chosen for grid

    .proj.fr.frsel.zg.datum configure -text $datum
    return
}

proc ParamGrid {gr data} {
    # choice of grid type has changed
    global Param $data INVTXT

    if { $Param(_guse) } {
	.proj.fr.frsel.m_guse configure -state normal
	set ${data}(_grid) $INVTXT([.proj.fr.frsel.m_guse cget -text])
	foreach w [winfo children .proj.fr.frsel.zg] {
	    catch {$w configure -state disabled}
	}
	foreach w [winfo children .proj.fr.frsel.zg.zb] {
	    catch {$w configure -state disabled}
	}
    } else {
	.proj.fr.frsel.m_guse configure -state disabled
	set ${data}(_grid) $gr
	foreach w [winfo children .proj.fr.frsel.zg] {
	    catch {$w configure -state normal}
	}
	foreach w [winfo children .proj.fr.frsel.zg.zb] {
	    catch {$w configure -state normal}
	}
	if { [winfo class .proj.fr.frsel.zg.cfd] != "Label" && \
		! $Param(_gfdatum) } {
	    .proj.fr.frsel.zg.datum configure -state disabled	    
	}
    }
    return
}

proc ParamGetCoord {param whc val} {
    # convert coordinate $val to signed degrees and check it
    #  $whc in {lat, long}
    #  $param is parameter name
    # return "empty" if empty, "nil" on error, otherwise
    #  list with value in signed degrees, heading (as letter) and
    #  absolute value in degrees
    global Param

    if { $val == "" } { return empty }
    if { $whc == "lat" } {
	set negh S ; set posh N ; set wC Lat
    } else { set negh W ; set posh E ; set wC Long }
    if { ! [Check$wC GMMessage $val $Param($param)] } { return nil }
    set val [Coord $Param($param) $val $negh]
    if { $val < 0 } {
	return [list $val $negh [expr -$val]]
    }
    return [list $val $posh $val]
}

proc ParamPFormt {fr whc pfmt param} {
    # change format of lat/long in dialog on the projection parameters
    #  $whc in {lat, long}
    #  $param is parameter name
    global TXT Param

    if { $pfmt == $Param($param) } { return }
    set val [ParamGetCoord $param $whc [$fr.val get]]
    if { $val == "nil" } { return }
    if { $val != "empty" } {
	$fr.val delete 0 end
	$fr.val insert 0 [lindex $val 1][ExtDegrees $pfmt [lindex $val 2]]
    }
    $fr.pfmt configure -text $TXT($pfmt)
    set Param($param) $pfmt
    return
}

proc ProjParamsRevert {proj data grid oldgr selfgrid} {
    # revert to initial values in dialog on the projection parameters
    #  $selfgrid is set if editing a projection defining a grid
    global $data MAPPROJDATA MAPPROJDTYPE TXT Param DLUNIT DISTUNIT Datum

    foreach p $MAPPROJDATA($proj) {
	switch -glob $MAPPROJDTYPE($p) {
	    lat -  long -  lat=* -
	    long=* {
		set fr .proj.fr.frsel.f$p
		$fr.pfmt configure -text $TXT(DDD)
		$fr.val delete 0 end
		$fr.val insert 0 [set ${data}($p)]
	    }
	    list=* -  list:* {
		.proj.fr.frsel.f$p.r[set ${data}($p)] invoke
	    }
	    reserved { }
	    default {
		.proj.fr.frsel.e$p delete 0 end
		.proj.fr.frsel.e$p insert 0 [set ${data}($p)]
	    }
	}
    }
    if { $grid } {
	set ${data}(_grid) $oldgr
	if { [set datum [set ${data}(datum)]] == "" } {
	    set datum $Datum
	    catch {.proj.fr.frsel.zg.cfd deselect}
	} else {
	    catch {.proj.fr.frsel.zg.cfd select}
	}
	.proj.fr.frsel.zg.datum configure -text $datum
	foreach e "efe efn zb.exn zb.exx zb.eyn zb.eyx \
		zb.elan zb.elax zb.elon zb.elox" {
	    .proj.fr.frsel.zg.$e delete 0 end
	}
	if { $selfgrid } {
	    ParamGridSet $oldgr
	} else {
	    .proj.fr.frsel.m_guse configure -text $TXT($oldgr)
	    .proj.fr.frsel.r_guse invoke
	    set Param(_gunit) $DISTUNIT
	    .proj.fr.frsel.zg.mu configure -text $DLUNIT($DISTUNIT,subdist)
	}
    }
    update idletasks
    return
}

## dialog for selecting datum and projection

proc ChooseDatumProjection {args} {
    # dialog for selection of datum and projection
    #  $args is "" or contains a list of latd,longd,datum to be projected
    # return 0 on failure, or projection internal name on success, in which
    #  case datum and projection parameters are stored in global array CDPData
    global Datum MapProjection CDPDatum CDPProj CDPProjT CDPData ProjRes \
	    MAPPROJDATA MAPPARTPROJ MAPPARTPDATA MAPPROJAUX EPOSX EPOSY \
	    TXT COLOUR DATUMWIDTH

    set CDPDatum $Datum ; set CDPProj $MapProjection
    set CDPProjT $TXT(PRJ$MapProjection)
    toplevel .cdp
    wm protocol .cdp WM_DELETE_WINDOW { set ProjRes ok }
    wm title .cdp "$TXT(projection)/GPS Manager"
    wm transient .cdp
    wm geometry .cdp +$EPOSX+$EPOSY
    bind .cdp <Key-Return> { set ProjRes ok }

    frame .cdp.fr -relief flat -borderwidth 5 -bg $COLOUR(selbg)
    label .cdp.fr.title -text "???" -relief sunken

    set frs .cdp.fr.fbs1
    frame $frs  -relief flat -borderwidth 0
    menubutton $frs.datum -textvariable CDPDatum -menu $frs.datum.m \
	    -width $DATUMWIDTH
    menu $frs.datum.m \
	    -postcommand "FillDatumMenu $frs.datum.m CDPChangeDatum"
    menubutton $frs.cdp -textvariable CDPProjT -width 20 -menu $frs.cdp.m
    menu $frs.cdp.m -postcommand "FillProjsMenu $frs.cdp.m \
		                     MAPKNOWNPROJS CDPChangeProjection"

    set frb .cdp.fr.fbs2
    frame $frb  -relief flat -borderwidth 0
    button $frb.ok -text $TXT(ok) -command { set ProjRes ok }
    button $frb.cnc -text $TXT(cancel) -command { set ProjRes cancel }

    pack $frs.datum $frs.cdp -side left -padx 2
    pack $frb.ok $frb.cnc -side left -padx 2
    pack .cdp.fr.title $frs $frb -side top -pady 5
    pack .cdp.fr
    update idletasks
    # cannot use RaiseWindow because of menus
    # RaiseWindow .cdp
    Raise .cdp
    set gs [grab current]
    grab .cdp
    tkwait variable ProjRes
    DestroyRGrabs .cdp $gs
    if { $ProjRes == "cancel" } {
	return 0
    }
    # initialize projection parameters
    catch {unset CDPData}
    set CDPData(datum) $CDPDatum

    if { [catch {set mp $MAPPARTPROJ($CDPProj)}] } {
	set ps [lindex $args 0]
	if { $ps == "" } {
	    # main projection: ask parameters
	    foreach p $MAPPROJDATA($CDPProj) {
		set CDPData($p) ""
	    }
	    if { [ProjParams set $CDPProj CDPData] == 0 } { return 0 }
	} else {
	    ProjInit $CDPProj CDPData $CDPDatum $ps
	}
	set mp $CDPProj
    } else {
	foreach e $MAPPROJDATA($mp) v $MAPPARTPDATA($CDPProj) {
	    set CDPData($e) $v
	}
    }
    if { [lsearch -exact $MAPPROJAUX $mp] != -1 } {
	Proj${mp}ComputeAux CDPData $CDPDatum
    }
    set CDPData(main_proj) $mp
    return $CDPProj
}

proc CDPChangeDatum {datum args} {
    # record selected datum
    #  $args not used but is needed as this is called-back from a menu
    global CDPDatum

    set CDPDatum $datum
    return
}

proc CDPChangeProjection {proj} {
    # record selected projection
    global CDPProj CDPProjT TXT

    set CDPProj $proj ; set CDPProjT $TXT(PRJ$proj)
    return
}

## menus with projections and grids

proc FillProjsMenu {menu listname comm} {
    # create menu entries from projections in list given by $listname for
    #  calling command whose argument is each projection
    global $listname MAXMENUITEMS TXT

    if { [winfo exists $menu] } {
	$menu delete 0 end
    }
    set n 0 ; set m 0
    foreach pr [set $listname] {
	if { $n > $MAXMENUITEMS } {
	    $menu add cascade -label "$TXT(more) ..." -menu $menu.m$m
	    set menu $menu.m$m ; destroy $menu
	    menu $menu -tearoff 0
	    set n 0 ; incr m
	}
	$menu add command -label $TXT(PRJ$pr) -command "$comm $pr"
	incr n
    }
    return
}

## initialization of particular cases of main projections

proc MapInitPartProj {proj data datum ps} {
    # initialize particular projection for the map
    #  $data is name of global array for parameters
    # return projection of first position

    set mp [ProjInitPart $proj $data $datum]
    set p [lindex $ps 0]
    set latd [lindex $p 0] ; set longd [lindex $p 1] ; set pdatum [lindex $p 2]
    return [Proj${mp}Point $data $latd $longd $pdatum]
}

proc ProjInitPart {proj data datum} {
    # set fixed parameters for a particular case of a main projection
    #  $data is name of global array for parameters
    # return name of main projection
    global $data MAPPARTPROJ MAPPARTPDATA MAPPROJDATA MAPPROJAUX

    set mp $MAPPARTPROJ($proj)
    foreach e $MAPPROJDATA($mp) v $MAPPARTPDATA($proj) {
	set ${data}($e) $v
    }
    if { [lsearch -exact $MAPPROJAUX $mp] != -1 } {
	Proj${mp}ComputeAux $data $datum
    }
    set ${data}(datum) $datum
    return $mp
}

## setting up the map coordinates format and datum

proc ChangeMapPFormat {pfmt} {
    # change position format to use with map cursor coordinates
    # if it is a grid with a fixed datum, change $MapPFDatum
    global GRIDS ZGRID Datum MapPFormat MapPFDatum MapPFNeedsDatum MapEmpty

    set MapPFormat $pfmt
    set MapPFNeedsDatum 0
    if { [lsearch -exact $GRIDS $pfmt] != -1 && ! $ZGRID($pfmt) } {
	global GRD$pfmt
	if { [set datum [set GRD${pfmt}(datum)]] != "" } {
	    set MapPFDatum $datum
	    set MapPFNeedsDatum 1
	}	
    }
    return
}

proc ChangeMPFDatum {datum args} {
    # change datum of map cursor coordinates
    #  $args is not used but is needed as this is called-back from
    #  menus built by proc FillDatumMenu (geod.tcl)
    # do nothing if datum is not appropriate to current cursor position format
    global MapPFDatum MapPFNeedsDatum MapEmpty

    if { $MapPFDatum != $datum } {
	if { $MapPFNeedsDatum } { return }
	set MapPFDatum $datum
    }
    return
}

## setting up of a projection for the map

proc MapProjectionIs {proj} {
    # set global variables so that map projection is $proj
    #  $proj in $MAPKNOWNPROJS
    global MapProjection MapProjTitle MapProjInitProc GRIDFOR MAPPARTPROJ TXT \
	    CMDLINE MapPFNeedsDatum MapPFDatum

    set MapProjection $proj ; set MapProjTitle $TXT(PRJ$proj)
    if { ! $CMDLINE } {
	ChangeMapPFormat $GRIDFOR($proj)
	if { $MapPFNeedsDatum } { ChangeMapDatum $MapPFDatum }
    }
    # names of projection procs
    if { [catch {set mp $MAPPARTPROJ($proj)}] } {
	# main projection
	set MapProjInitProc ProjInit
    } else {
	# particular case of a main projection
	set MapProjInitProc MapInitPartProj
	set proj $mp
    }
    foreach a "Point Invert" {
	global MapProj${a}Proc

	set MapProj${a}Proc Proj${proj}$a
    }
    return
}

proc ChangeMapDatum {datum args} {
    # change map datum
    #  $args is not used but is needed as this is called-back from
    #  menus built by proc FillDatumMenu (geod.tcl)
    # to be called only when map is empty
    global Datum

    set Datum $datum
    return
}

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

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
	}
	default {
	    BUG "no initialization proc for projection"
	}
    }
    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" that only has 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 -  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-$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_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 x [expr $A*$cosX]
	    set y [expr $A*($cosX1*$sinX-$sinX1*$cosX*$coslam)]
	}
	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 ${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)]
	    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 -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  -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)]
}


