#
# This file is part of:
#
#  gpsman --- GPS Manager: a manager for GPS receiver data
#
# Copyright (c) 1998-2013 Miguel Filgueiras migfilg@t-online.de
#
#    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 3 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: compute.tcl
#  Last change:  6 October 2013
#
# Includes contributions by
#  - Matt Martin (matt.martin _AT_ ieee.org) marked "MGM contribution"
#  - Valere Robin (valere.robin _AT_ wanadoo.fr) marked "VR contribution"
#

## Some formulae kindly supplied by
#     Luisa Bastos, Universidade do Porto
#     Gil Goncalves, Universidade de Coimbra
#  Computation of area of spherical polygon adapted from sph_poly.c in
#   "Graphics Gems IV", edited by Paul Heckbert, Academic Press, 1994.
#  Formula for ellipsoid radius from
#   "Ellipsoidal Area Computations of Large Terrestrial Objects"
#   by Hrvoje Lukatela
#   http://www.geodyssey.com/papers/ggelare.html
#  Algorithm for finding a point at given distance and bearing from another
#   taken from the program "forward" available from
#   ftp://www.ngs.noaa.gov/pub/pcsoft/for_inv.3d/
#   It uses the modified Rainsford's Method with Helmert's elliptical terms,
#   and is effective in any azimuth and at any distance short of antipodal.
##

### positions and coordinates

proc Coord {pformt coord negh} {
    # convert coordinate $coord in format $pformt, with negative heading $negh,
    #  to signed degrees
    # $pformt in {DMS, DMM, DDD, GRA}

    set coord [string trim $coord]
    set sign 1
    set h [string index "$coord" 0]
    if { ! [regexp {[0-9]} $h] } {
	if { $h == $negh || $h == "-" } { set sign -1 }
	set coord [string range "$coord" 1 end]
    }
    switch $pformt {
	DMS {
	    scan "$coord" "%d %d %f" d m s
	    return [expr $sign*($d+$m/60.0+$s/3600.0)]
	}
	DMM {
	    scan "$coord" "%d %f" d m
	    return [expr $sign*($d+$m/60.0)]
	}
	DDD {
	    return [expr $sign*$coord]
	}
	GRA {
	    return [expr $sign*0.9*$coord]
	}
    }
}

proc FormatPosition {latd longd datum reqpfmt pfdatum args} {
    # produce formatted position from $latd and $longd in decimal
    #  degrees for $datum
    #  $reqpfmt is the required position format
    #  $pfdatum is the datum to use with this format
    #     if $pfdatum=="", $datum is used unless $reqpfmt needs a fixed datum
    #     otherwise $pfdatum must be compatible with $reqpfmt; if it is
    #     not a call to BUG will be made
    #  $args is either "" or is an alternative format to be used with $datum
    #     when conversion to $reqpfmt fails (out of grid scope)
    # result is a list with: position representation, position format
    #  used and datum used
    # if $args=="" and conversion to $reqpfmt yields an out of grid scope
    #  result ("--" in the zone or x-field), this result is used in
    #  building the position representation
    # a position representation is a list whose first two elements
    #  are lat and long in degrees, and whose rest is dependent on format
    #  type:
    #       ==latlong: lat and long in external format
    #       ==utm: zones east and north and x y coordinates
    #       ==grid: zone name and x y coordinates
    #       ==nzgrid: x y coordinates
    #       ==mh: Maidenhead-locator (6 characters)
    global POSTYPE

    set la $latd ; set lo $longd
    if { $pfdatum == "" } {
	set udatum $datum
    } else {
	set udatum $pfdatum
	if { $pfdatum != $datum } {
	    foreach "la lo" [ToDatum $la $lo $datum $pfdatum] {	break }
	}
    }
    switch [set ptype $POSTYPE($reqpfmt)] {
	latlong {
	    set p [FormatLatLong $la $lo $reqpfmt]
	}
	utm {
	    set p [linsert [DegreesToUTM $la $lo $udatum] 0 $la $lo]
	}
	grid -  nzgrid {
	    if { [set gd [GridDatum $reqpfmt $udatum]] != $udatum } {
		if { $pfdatum != "" } {
		    BUG calling FormatPosition with wrong pfdatum
		}
		set udatum $gd
		foreach "la lo" [ToDatum $latd $longd $datum $gd] { break }
	    }
	    if { $ptype == "grid" } {
		set p [DegreesTo$reqpfmt $la $lo $udatum]
	    } else { set p [DegreesToNZGrid $reqpfmt $la $lo $udatum] }
	    if { [lindex $p 0] == "--" && $args != "" } {
		# out of scope
		return [FormatPosition $latd $longd $datum $args ""]
	    }
	    set p [linsert $p 0 $la $lo]
	}
	mh {
	    set p [list $la $lo [DegreesToMHLoc $la $lo]]
	}
    }
    return [list $p $reqpfmt $udatum]
}

proc FormatLatLong {latd longd pfmt} {
    # produce formatted position from $latd and $longd in decimal
    #  degrees in latitude/longitude position format
    # result is the position representation: a list whose first two elements
    #  are lat and long in degrees, followed by lat and long in external format
    global COUTFMT

    if { $latd < 0 } {
	# VR contribution: 0-$latd instead of -$latd; same for $longd
	set lax [expr 0-$latd] ; set hlat S
    } else { set lax $latd ; set hlat N }
    if { $longd < 0 } {
	set lox [expr 0-$longd] ; set hlng W
    } else { set lox $longd ; set hlng E }
    if { $pfmt == "GRA" } {
	set lax [format $COUTFMT(GRA) [expr $latd/0.9]]
	set lox [format $COUTFMT(GRA) [expr $longd/0.9]]
	return [list $latd $longd $lax $lox]
    }
    return [list $latd $longd "$hlat[ExtDegrees $pfmt $lax]" \
		"$hlng[ExtDegrees $pfmt $lox]"]
}

proc ExtDegrees {pformt degs} {
    # from signed degrees $degs to external representation of format $pformt
    # $pformt in {DMS, DMM, DDD, DMSsimpl}
    #  DMSsimpl is similar to DMS, but for values less than 1 degree the
    #  representation is either MM'SS.S" or SS.S"
    global COUTFMT

    switch -glob $pformt {
	DMS* {
	      if { $degs < 0 } {
		  set degs [expr -$degs] ; set sign -1
	      } else { set sign 1 }
	      set d [expr int($degs)]
	      set degs [expr ($degs-$d)*60]
	      set d [expr $sign*$d]
	      set m [expr int($degs)]
	      set s [expr ($degs-$m)*60]
	      if { $s > 59.95 } { set s 0 ; incr m }
	      if { $m > 59 } { set m 0 ; incr d }
	      set sf $COUTFMT(sec)
	      if { $d > 0 || $pformt == "DMS" } {
		  return [format "%d %02d $sf" $d $m $s]
	      }
	      if { $m > 0 } {
		  return [format "%02d'$sf" $m $s]
	      }
	      return [format "$sf\"" $s]
	}
	DMM {
	      if { $degs < 0 } {
		  set degs [expr -$degs] ; set sign -1
	      } else { set sign 1 }
	      set d [expr int($degs)]
	      set degs [expr ($degs-$d)*60]
	      set d [expr $sign*$d]
	      if { $degs > 59.995 } { set degs 0 ; incr d }
	      return [format "%d $COUTFMT(min)" $d $degs]
	}
	DDD {
	      return [format $COUTFMT(deg) $degs]
	}
    }
}

proc ChangeCoordSign {pos} {

    if { $pos != "N" } { return W }
    return S
}

proc ConvertPos {pfrmt zone x y datum npfrmt} {
    # create position representation in format $npfrmt from position
    #  $x, $y (as decimal numbers) in format $pfrmt
    #         if $pfrmt in {DDD, GRA} take $x as long and $y as lat
    #  $zone is either a list with ze and zn for UTM/UPS, or the zone
    #    name for other grids (and is meaningless for other formats
    # values of $x, $y are checked, as well as suitability of $datum
    #  and of validity of new format
    # return -1 on error
    global POSTYPE MESS TXT COUTFMT

    switch [set ptype $POSTYPE($pfrmt)] {
	latlong {
	    # must be either DDD or GRA
	    if { $pfrmt == "GRA" } {
		set x [expr $x*0.9] ; set y [expr $y*0.9]
	    }
	    if { $y < 0 } {
		set la [expr -$y] ; set hlat S
	    } else {
		set la $y ; set hlat N
	    }
	    if { $x < 0 } {
		set lo [expr -$x] ; set hlng W
	    } else {
		set lo $x ; set hlng E
	    }
	    if { ! [CheckCoord GMMessage $pfrmt $la N 90] || \
		    ! [CheckCoord GMMessage $pfrmt $lo E 180] } {
		return -1
	    }
	    if { $npfrmt == "DDD" } {
		return [list $y $x "${hlat}[format $COUTFMT(deg) $la]" \
			           "${hlng}[format $COUTFMT(deg) $lo]"]
	    }
	    set la $y ; set lo $x
	}
	utm {
	    foreach "ze zn" $zone {}
	    foreach "la lo" [UTMToDegrees $ze $zn $x $y $datum] { break }
	    if { $npfrmt == $pfrmt } {
		return [list $la $lo $ze $zn $x $y]
	    }
	}
	grid -  nzgrid {
	    if { [BadDatumFor $pfrmt $datum GMMessage] != 0 } {
		return -1
	    }
	    if { [set p [GridToDegrees $pfrmt $zone $x $y $datum]] == 0 } {
		GMMessage $MESS(outofgrid)
		return -1
	    }
	    if { $npfrmt == $pfrmt } {
		if { $ptype == "grid" } {
		    return [linsert $p end $zone $x $y]
		} else { return [linsert $p end $x $y] }
	    }
	}
	default {
	    BUG calling ConvertPos with wrong format
	}
    }
    # $datum is suitable for $npfrmt
    set p [lindex [FormatPosition $la $lo $datum $npfrmt $datum] 0]
    if { [lindex $p 2] == "--" } {
	GMMessage $MESS(outofgrid)
	return -1
    }
    return $p
}

### distances and geographic bearings
# formulae from Spherical Trignometry
# the Law of Cosines, used below to compute distances and bearings, is known
#  to be unsuitable for computations because of excessive round-off errors;
#  alternative procedures are supplied in file acccomp.tcl and are loaded
#  if $ACCFORMULAE (accurate formulae option) is set

proc ComputeDist {p1 p2 datum} {
    # distance (km) between positions $p1 and $p2 with same datum
            # formulae kindly supplied by Luisa Bastos, Universidade do Porto
            #   and Gil Goncalves, Universidade de Coimbra

    set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1]
    set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1]
    if { $lad1==$lad2 && $lod1==$lod2 } { return 0 }
    set la1 [expr $lad1*0.01745329251994329576]
    set lo1 [expr $lod1*0.01745329251994329576]
    set la2 [expr $lad2*0.01745329251994329576]
    set lo2 [expr $lod2*0.01745329251994329576]
    set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)]
    if { $x >= 1 } { return 0 }
    return [expr 1e-3*[lindex [EllipsdData $datum] 0]*acos($x)]
}

proc ComputeBear {p1 p2 datum} {
    # bearing from position $p1 to $p2 with same datum
    #  $datum not used here but kept for compatibility with replacement
    #    proc in acccomp.tcl
            # formulae kindly supplied by Luisa Bastos, Universidade do Porto
            #   and Gil Goncalves, Universidade de Coimbra

    set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1]
    set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1]
    if { $lad1==$lad2 && $lod1==$lod2 } { return 0 }
    set la1 [expr $lad1*0.01745329251994329576]
    set lo1 [expr $lod1*0.01745329251994329576]
    set la2 [expr $lad2*0.01745329251994329576]
    set lo2 [expr $lod2*0.01745329251994329576]
    set da [expr $la2-$la1] ; set do [expr $lo2-$lo1]
    if { [expr abs($da)] < 1e-20 } {
 	if { [expr abs($do)] < 1e-20 } {
 	    set b 0
 	} elseif { $do < 0 } {
 	    set b 270
 	} else { set b 90 }
    } elseif { [expr abs($do)] < 1e-20 } {
 	if { $da < 0 } {
 	    set b 180
 	} else { set b 0 }
    } else {
 	set b [expr round(atan2(sin($do), \
 		                tan($la2)*cos($la1)-sin($la1)*cos($do)) \
 			  *57.29577951308232087684)]
	if { $b < 0 } {
	    if { $do < 0 } { incr b 360 } else { incr b 180 }
	} elseif { $do < 0 } { incr b 180 }
    }
    return $b
}    

proc ComputeDistBear {p1 p2 datum} {
    # distance between and bearing from positions $p1 and $p2 with same datum
            # formulae kindly supplied by Luisa Bastos, Universidade do Porto
            #   and Gil Goncalves, Universidade de Coimbra

    set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1]
    set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1]
    if { $lad1==$lad2 && $lod1==$lod2 } { return "0 0" }
    set la1 [expr $lad1*0.01745329251994329576]
    set lo1 [expr $lod1*0.01745329251994329576]
    set la2 [expr $lad2*0.01745329251994329576]
    set lo2 [expr $lod2*0.01745329251994329576]
    set a [lindex [EllipsdData $datum] 0]
           # distance (km)
    set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)]
    if { $x < 1 } {
	set d [expr 1e-3*$a*acos($x)]
    } else { set d 0 }
           # bearing
    set da [expr $la2-$la1] ; set do [expr $lo2-$lo1]
    if { [expr abs($da)] < 1e-20 } {
 	if { [expr abs($do)] < 1e-20 } {
 	    set b 0
 	} elseif { $do < 0 } {
 	    set b 270
 	} else { set b 90 }
    } elseif { [expr abs($do)] < 1e-20 } {
 	if { $da < 0 } {
 	    set b 180
 	} else { set b 0 }
    } else {
 	set b [expr round(atan2(sin($do), \
 		                tan($la2)*cos($la1)-sin($la1)*cos($do)) \
 			  *57.29577951308232087684)]
	if { $b < 0 } {
	    if { $do < 0 } { incr b 360 } else { incr b 180 }
	} elseif { $do < 0 } { incr b 180 }
    }
    return [list $d $b]
}

proc SetDatumData {datum} {
    # this proc sets datum parameters as global variables for repeated use
    #  in conversions; see procs ComputeDistFD, ComputeDistBearFD
    global DatumA DatumF FDDatum

    if { $FDDatum != $datum } {
	set dt [EllipsdData $datum]
	set DatumA [lindex $dt 0] ; set DatumF [lindex $dt 1]
	set FDDatum $datum
    }
    return
}

proc ComputeDistFD {p1 p2} {
    # compute distance (km) between positions $p1 and $p2 assuming datum
    #  parameters where set by calling SetDatumData
    global DatumA
            # formulae kindly supplied by Luisa Bastos, Universidade do Porto
            #   and Gil Goncalves, Universidade de Coimbra

    set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1]
    set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1]
    if { $lad1==$lad2 && $lod1==$lod2 } { return 0 }
    set la1 [expr $lad1*0.01745329251994329576]
    set lo1 [expr $lod1*0.01745329251994329576]
    set la2 [expr $lad2*0.01745329251994329576]
    set lo2 [expr $lod2*0.01745329251994329576]
    set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)]
    if { $x >= 1 } { return 0 }
    return [expr 1e-3*$DatumA*acos($x)]
}

proc ComputeDistBearFD {p1 p2} {
    # compute distance (km) between and bearing from positions $p1 and $p2
    #  assuming datum parameters where set by calling SetDatumData
    global DatumA
            # formulae kindly supplied by Luisa Bastos, Universidade do Porto
            #   and Gil Goncalves, Universidade de Coimbra

    set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1]
    set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1]
    if { $lad1==$lad2 && $lod1==$lod2 } { return "0 0" }
    set la1 [expr $lad1*0.01745329251994329576]
    set lo1 [expr $lod1*0.01745329251994329576]
    set la2 [expr $lad2*0.01745329251994329576]
    set lo2 [expr $lod2*0.01745329251994329576]
           # distance (km)
    set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)]
    if { $x < 1 } {
	set d [expr 1e-3*$DatumA*acos($x)]
    } else { set d 0 }
           # bearing
    set da [expr $la2-$la1] ; set do [expr $lo2-$lo1]
    if { [expr abs($da)] < 1e-20 } {
 	if { [expr abs($do)] < 1e-20 } {
 	    set b 0
 	} elseif { $do < 0 } {
 	    set b 270
 	} else { set b 90 }
    } elseif { [expr abs($do)] < 1e-20 } {
 	if { $da < 0 } {
 	    set b 180
 	} else { set b 0 }
    } else {
 	set b [expr round(atan2(sin($do), \
 		                tan($la2)*cos($la1)-sin($la1)*cos($do)) \
 			  *57.29577951308232087684)]
	if { $b < 0 } {
	    if { $do < 0 } { incr b 360 } else { incr b 180 }
	} elseif { $do < 0 } { incr b 180 }
    }
    return [list $d $b]
}

proc CompWPDistBear {wp1 wp2} {
    # distance (km) between and bearing from WPs with names $wp1, $wp2
    global WPPFrmt WPPosn WPDatum

    set i1 [IndexNamed WP $wp1]
    set i2 [IndexNamed WP $wp2]
    if { $i1<0 || $i2<0 } { return "--- ---" }
    if { $i1 == $i2 } { return "0 0" }
    return [CompDistBearDatums $WPPosn($i1) $WPDatum($i1) \
	                          $WPPosn($i2) $WPDatum($i2)]
}

proc CompDistBearDatums {p1 d1 p2 d2} {
    # distance (km) between and bearing from positions $p1 and $p2
    #  with datums $d1 and $d2
    if { $d1 != $d2 } {
	set p2 [ToDatum [lindex $p2 0] [lindex $p2 1] $d2 $d1]
    }
    return [ComputeDistBear $p1 $p2 $d1]
}

## point at given distance and bearing from another point

proc CoordsAtDistBear {pos dist bear datum} {
    # compute coordinates of point at given distance (km) and bearing
    #  (0-360 degrees) from another point
    # algorithm taken from the program "forward" available from
    #  ftp://www.ngs.noaa.gov/pub/pcsoft/for_inv.3d/
    #  It uses the modified Rainsford's Method with Helmert's elliptical
    #  terms, and is effective in any azimuth and at any distance short of
    #  antipodal.
    #  $dist is constrained here to be less than pi*semi-major axis
    global MESS

    foreach "lat0 long0" $pos { break }
    foreach "a f" [EllipsdData $datum] { break }
    if { [set dist [expr $dist*1000]] >= 3.1415926*$a } {
	GMMessage $MESS(distlarge)
	return [list $lat $long]
    }
    set phi [expr $lat0*0.01745329251994329576]
    set lam [expr $long0*0.01745329251994329576]
    set azm [expr $bear*0.01745329251994329576]
    set r [expr 1-$f]
    set tu [expr $r*sin($phi)/cos($phi)]
    set sf [expr sin($azm)] ; set cf [expr cos($azm)]
    set baz 0
    if { [set cf [expr cos($azm)]] != 0 } {
	set baz [expr atan2($tu,$cf)*2]
    }
    set cu [expr 1.0/sqrt($tu*$tu+1)]
    set su [expr $tu*$cu]
    set sa [expr $cu*$sf]
    set c2a [expr -$sa*$sa+1]
    set x [expr sqrt((1.0/$r/$r-1)*$c2a+1)+1]
    set x [expr ($x-2.0)/$x]
    set c [expr 1-$x]
    set c [expr ($x*$x/4.0+1.0)/$c]
    set d [expr (0.375*$x*$x-1)*$x]
    set y [set tu [expr 1.0*$dist/$r/$a/$c]]
    while 1 {
	set sy [expr sin($y)] ; set cy [expr cos($y)]
	set cz [expr cos($baz+$y)]
	set e [expr $cz*$cz*2-1]
	set c $y
	set x [expr $e*$cy]
	set y [expr ((($sy*$sy*4-3)*($e+$e-1)*$cz*$d/6.0+$x)*$d/4.0-$cz)* \
		    $sy*$d+$tu]
	if { abs($y-$c) <= 5e-14 } { break }
    }
    set baz [expr $cu*$cy*$cf-$su*$sy]
    set c [expr $r*sqrt($sa*$sa+$baz*$baz)]
    set d [expr $su*$cy+$cu*$sy*$cf]
    set phi2 [expr atan2($d,$c)]
    set c [expr $cu*$cy-$su*$sy*$cf]
    set x [expr atan2($sy*$sf,$c)]
    set c [expr ((-3*$c2a+4)*$f+4)*$c2a*$f/16.0]
    set d [expr (($e*$cy*$c+$cz)*$sy*$c+$y)*$sa]
    set lam2 [expr $lam+$x-(1-$c)*$d*$f]
    # set baz [expr atan2($sa,$baz)+3.14159265358979323846]

    set lat1 [expr $phi2*57.29577951308232087684]
    set long1 [expr $lam2*57.29577951308232087684]
    return [list $lat1 $long1]
}

### area

proc ProjectedArea {wpixs} {
    # compute area of polygon whose boundary is the polyline formed by the
    #  WPs in given list
    # polygon cannot be self-intersecting (no test on this!)
    # return value in square km
    global WPPosn WPDatum ASKPROJPARAMS

    set ix [lindex $wpixs 0]
    set datum $WPDatum($ix) ; set p $WPPosn($ix)
    set ps [list [list [set lat [lindex $p 0]] [lindex $p 1] $datum]]
    foreach ix [lreplace $wpixs 0 0] {
	set dt $WPDatum($ix) ; set p $WPPosn($ix)
	if { $dt != $datum } {
	    set p [ToDatum [lindex $p 0] [lindex $p 1] $dt $datum]
	}
	lappend ps [list [lindex $p 0] [lindex $p 1] $datum]
    }
    if { $lat>=-80 && $lat<=84 } {
	set proj TM
    } else {
	# in fact this will use UPS
	set proj UTM
    }
    set ask $ASKPROJPARAMS
    set ASKPROJPARAMS 0
    set xy [ProjInit $proj AProj $datum $ps]
    set ASKPROJPARAMS $ask
    set xmin [set x0 [lindex $xy 0]] ; set ymin [set y0 [lindex $xy 1]]
    set xs "" ; set ys ""
    foreach p [lreplace $ps 0 0] {
	set xy [eval Proj${proj}Point AProj $p]
	set x [lindex $xy 0] ; set y [lindex $xy 1]
	lappend xs $x ; lappend ys $y
	if { $x < $xmin } { set xmin $x }
	if { $y < $ymin } { set ymin $y }
    }
    set sum 0
    set x0 [set x [expr $x0-$xmin]] ; set y0 [set y [expr $y0-$ymin]]
    foreach x1 $xs y1 $ys {
	set x1 [expr $x1-$xmin] ; set y1 [expr $y1-$ymin]
	set sum [expr $sum+$x*$y1-$x1*$y]
	set x $x1 ; set y $y1
    }
    set sum [expr $sum+$x*$y0-$x0*$y]
    return [expr abs(0.5e-6*$sum)]
}

proc SphericalArea {wpixs} {
    # compute area of polygon whose boundary is the polyline formed by the
    #  WPs in given list
    # computation is based on a spherical approximation
    # return negative value if there are precision errors
    #  otherwise value in square km
    global WPPosn WPDatum

    set ix [lindex $wpixs 0]
    set datum $WPDatum($ix) ; set p [set p0 $WPPosn($ix)]
    set phi0 [set phi [expr [lindex $p 0]*0.01745329251994329576]]
    set lam0 [set lam [expr [lindex $p 1]*0.01745329251994329576]]
    set c0 [expr cos($phi0)]
    # computation of mean radius of the ellipsoid at a point
    #  as described in Hrvoje Lukatela, "Ellipsoidal Area Computations of
    #  Large Terrestrial Objects", http://www.geodyssey.com/papers/ggelare.html
    set dt [EllipsdData $datum]
    set a [lindex $dt 0] ; set b [expr $a*(1-[lindex $dt 1])]
    set s [expr sin($phi)]
    set a2 [expr $a*$a]
    # radius in km
    set r [expr 1e-3*$a2*$b/($a2*$c0*$c0+$b*$b*$s*$s)]
    # computation of area adapted from sph_poly.c in
    #  "Graphics Gems IV", edited by Paul Heckbert, Academic Press, 1994.
    set srad 0
    SetDatumData $datum
    foreach ixn [lreplace $wpixs 0 0] {
	set pos $WPPosn($ixn)
	if { $WPDatum($ixn) != $datum } {
	    set pos [ToDatum [lindex $pos 0] [lindex $pos 1] \
		    $WPDatum($ixn) $datum]
	}
	set phi1 [expr [lindex $pos 0]*0.01745329251994329576]
	set lam1 [expr [lindex $pos 1]*0.01745329251994329576]
	set c1 [expr cos($phi1)]
	if { $lam0 != $lam1 } {
	    set HavA [expr ((1-cos($phi1-$phi0))/2.0)+ \
		           ((1-cos($lam1-$lam0))/2.0)*$c0*$c1]
	    set A [expr 2*asin(sqrt($HavA))]
	    set B [expr 1.5707963267948966192313-$phi1]
	    set C [expr 1.5707963267948966192313-$phi0]
	    set S [expr 0.5*($A+$B+$C)]
            set T [expr tan($S/2.0)*tan(($S-$A)/2.0)* \
		    tan(($S-$B)/2.0)*tan(($S-$C)/2.0)]
	    if { abs($T) < 1e-8 } { return -1 }
	    set E [expr abs(4*atan(sqrt(abs($T))))]
	    if { $lam1 < $lam0 } { set E [expr -$E] }
	    set srad [expr $srad+$E]
	}
	set phi0 $phi1 ; set lam0 $lam1 ; set c0 $c1
    }
    set lam1 $lam
    if { $lam0 != $lam1 } {
	set phi1 $phi ; set c1 [expr cos($lam1)]
	set HavA [expr ((1-cos($phi1-$phi0))/2.0)+ \
		           ((1-cos($lam1-$lam0))/2.0)*$c0*$c1]
	set A [expr 2*asin(sqrt($HavA))]
	set B [expr 1.5707963267948966192313-$phi1]
	set C [expr 1.5707963267948966192313-$phi0]
	set S [expr 0.5*($A+$B+$C)]
	set T [expr tan($S/2.0)*tan(($S-$A)/2.0)* \
		    tan(($S-$B)/2.0)*tan(($S-$C)/2.0)]
	if { abs($T) < 1e-8 } { return -1 }
        set E [expr abs(4*atan(sqrt(abs($T))))]
	if { $lam1 < $lam0 } { set E [expr -$E] }
	set srad [expr $srad+$E]
    }
    return [expr abs($r*$r*$srad)]
}

### datums

proc ChangeTPsDatum {tps datum1 datum2} {
    # convert position of TR points on list $tps from $datum1 to $datum2
    # change only the first 4 elements (position) of each TP representation
    # this depends on DataIndex (see setup.tcl, Storage), not used here
    #  for speed

    if { $datum1 == $datum2 } { return $tps }
    set l ""
    foreach tp $tps {
	set np [lindex [FormatPosition [lindex $tp 0] [lindex $tp 1] \
			    $datum1 DMS $datum2] 0]
	lappend l [concat $np [lrange $tp 4 end]]
    }
    return $l
}

proc ChangeLPsDatum {lps odatum datum pformt} {
    # change datum of LN list of points
    # change only the first element (position) of each LP representation
    # this depends on DataIndex (see setup.tcl, Storage), not used here
    #  for speed
    # return -1 if the current position format requires a fixed datum
    #  different from $datum
    global MESS POSTYPE

    if { $odatum == $datum } { return $lps }
    if { [set gd [BadDatumFor $pformt $datum Ignore]] != 0 } {
	GMMessage [format $MESS(cantchggriddatum) $gd]
	return -1
    }
    set l ""
    foreach lp $lps {
	foreach "latd longd" [lindex $lp 0] { break }
	lappend l [lreplace $lp 0 0 \
		       [lindex [FormatPosition $latd $longd $odatum \
				    $pformt $datum] 0]]
    }
    return $l
}

### altitude

proc UserAltitude {altlst} {
    # return altitude in user units from internal representation
    global AltUnit ALSCALE

    if { $altlst == "" } { return "" }
    if { $AltUnit == "M" } { return [lindex $altlst 0] }
    if { $AltUnit == [lindex $altlst 2] } { return [lindex $altlst 1] }
    return [expr [lindex $altlst 0]/$ALSCALE]
}

proc AltitudeList {a} {
    # check value for altitude and return internal representation as list
    # return nil on error
    global AltUnit ALSCALE

    if { $a == "" } { return "" }
    if { ! [CheckSignedFloat Ignore $a] } { return nil }
    if { $AltUnit == "M" } { return $a }
    return [list [expr $a*$ALSCALE] $a $AltUnit]
}

### dates

proc Today {dformt} {
    # build representation of current date under format $dformt

    set dt [clock format [clock seconds] -format "%Y %m %d %H %M %S"]
    scan $dt %d%0d%0d%0d%0d%0d y m d h mn s
    return [FormatDate $dformt $y $m $d $h $mn $s]
}

proc TodayUTC {dformt} {
    # VR contribution
    # build representation of current date under format $dformt

    set dt [clock format [clock seconds] -format "%Y %m %d %H %M %S"]
    scan $dt %d%0d%0d%0d%0d%0d y m d h mn s
    return [eval FormatDate $dformt [LocalTimeAndUTC $y $m $d $h $mn $s Local]]
}

proc FormatDay {dformt y m d} {
    # build representation of date (without time of day) under format $dformt
    # (see proc FormatDate for date with time of day)
    #  $dformt in
    #   {YYYYMMDD, YYYY-MM-DD, MMDDYYYY, DDMMMYYYY, YYYY/MM/DD, ISO8601}
    # changes in proc FormatDate are likely to affect this proc!

    set cdate [FormatDate $dformt $y $m $d 0 0 0]
    if { $dformt == "YYYY/MM/DD" } {
	return [string range $cdate 0 9]
    }
    return [lindex $cdate 0]
}

proc FormatDate {dformt y m d h mn s} {
    # build representation of date (including time of day) under format $dformt
    # (see proc FormatDay for date without time of day)
    #  $dformt either in
    #   {YYYYMMDD, YYYY-MM-DD, MMDDYYYY, DDMMMYYYY, YYYY/MM/DD, ISO8601}
    #   or DDMMYY in which case the date will have no time of day
    # changes here must be reflected in proc FormatDay and in proc ScanDate
    # when adding new formats here, array DATEW must be updated in main.tcl
    global MONTHNAMES

    # make sure there are no leading 0's
    foreach v "y m d h mn s" {
	scan [set $v] %0d $v
    }
    set h [format "%02d:%02d:%02d" $h $mn $s]
    switch $dformt {
	DDMMYY {
	    # MGM contribution
	    return [format "%02d%02d%02d" $d $m [expr $y % 100]]
	}
	ISO8601 {
	    # VR contribution
	    return [format "%4d-%02d-%02dT%sZ" $y $m $d $h]
	}
    	YYYYMMDD { return [format "%4d.%02d.%02d %s" $y $m $d $h] }
	YYYY-MM-DD { return [format "%4d-%02d-%02d %s" $y $m $d $h] }
	MMDDYYYY { return [format "%02d/%02d/%4d %s" $m $d $y $h] }
	DDMMMYYYY {
	    incr m -1 ; set m [lindex $MONTHNAMES $m]
	    return [format "%02d-%s-%4d %s" $d $m $y $h]
	}
	YYYY/MM/DD { return [format "%4d/%02d/%02d-%s" $y $m $d $h] }
    }
}

proc FormatTime {secs args} {
    # build represention of seconds in hours:minutes:seconds format
    #  $secs may be a float or an integer
    # if $args!="" it gives the number of decimal places for the seconds field
    #  that defaults to 0 with $secs rounded to integer

    if { $args != "" && $args > 0 } {
	set fss "%02.${args}f" ; set fss0 "%8.${args}f"
    } else {
	set secs [expr round($secs)]
	set fss "%02d" ; set fss0 "%8d"
    }
    set s0 [expr int($secs)]
    set s [expr $s0%60] ; set x [expr ($s0-$s)/60]
    set s [expr $s+$secs-$s0]
    set mn [expr $x%60] ; set h [expr ($x-$mn)/60]
    if { $h > 0 } { return [format "%2d:%02d:$fss" $h $mn $s] }
    if { $mn > 0 } { return [format "%5d:$fss" $mn $s] }
    return [format $fss0 $s]
}

proc TimeToSecs {time} {
    # convert represention of seconds in hours:minutes:seconds format
    #  to seconds (hours and minutes fields may be missing)
    # return integer or float according to seconds field

    set n [llength [set l [split $time ":"]]]
    set s [lindex $l end]
    if { [string first 0 $s] == 0 } { set s [string range $s 1 end] }
    if { $n == 1 } { return $s }
    set ll ""
    foreach f $l {
	scan $f %0d f
	lappend ll $f
    }
    scan [lindex $l 0] %0d f
    if { $n == 2 } {
	return [expr $s+60*$f]
    }
    scan [lindex $l 1] %0d g
    return [expr $s+60*$g+3600*$f]
}

proc NowTZ {} {
    # current date under default format with time zone appended
    global DateFormat
    
    return "[Today $DateFormat] ([clock format 0 -format %Z])"
}

proc Now {} {
    # current date under default format
    global DateFormat
    
    return [Today $DateFormat]
}

array set DAYSOF {
     1 31   2 28   3 31  4 30  5 31  6 30  7 31  8 31  9 30
    10 31  11 30  12 31
}

proc ScanDate {date} {
    # convert date string to a list with year, month (in 1..12), day,
    #  hours, minutes and seconds
    # if there is time zone information, return local time
    # return empty list on error
    # see proc FormatDate for possible date formats
    # changes in proc FormatDate are likely to affect this proc!

    set date [string trim $date]
    if { [regexp {^([0-9][0-9])([0-9][0-9])([0-9][0-9])$} $date \
	      x d m y] } {
	# DDMMYY, no time of day
	#  YY taken as 2000+YY if YY<70, and 1900+YY otherwise
	foreach v "d m y" {
	    scan [set $v] %02d $v
	}
	if { $y < 70 } {
	    incr y 2000
	} else { incr y 1900 }
	return [list $y $m $d 0 0 0]
    }
    set goodn -1
    # includes VR contribution: decimal comma, empty tail
    if { [regexp {^ *([-0-9]+)T([0-9:]+([.,][0-9]+)?)([^ ]*) *$} \
	     $date m ymd hms m tail] } {
	# coping with ISO8601 format as in:
	#   2000-01-01T12:00:00Z
	#   2000-01-01T12:00:00.00000Z
	#   2004-06-07T13:51:27.2515650-07:00
	#   2009-09-23T07:02:46.773
	#   2009-09-23T07:02:46,773
	# must replace decimal comma by period
	regsub {,} $hms "." hms
	set n 1
	if { [scan $hms %02d:%02d:%0f h mn s] == 3 && \
		 [scan $ymd %4d-%02d-%02d y m d] == 3 } {
	    if { [set s [expr round($s)]] > 59 } {
		incr s -60
		if { [incr mn] > 59 } {
		    incr mn -60
		    if { [incr h] > 23 } {
			incr h -24
			foreach "y m d" [NextDay $y $m $d] {}
		    }
		}
	    }
	    if { [scan $tail %03d:%02d hz mz] == 2 } {
		set goodn 1
		# compute UTC date
		if { $hz > 0 } { set mz [expr -$mz] }
		set hz [expr -$hz]
		if { [incr mn $mz] > 59 } {
		    incr h ; incr mn -60
		} elseif { $mn < 0 } {
		    incr h -1 ; incr mn 60
		}
		if { [incr h $hz] > 23 } {
		    incr h -24
		    foreach "y m d" [NextDay $y $m $d] {}
		} elseif { $h < 0 } {
		    incr h 24
		    foreach "y m d" [PreviousDay $y $m $d] {}
		}
		# convert to local date
		foreach "y m d h mn s" \
		    [LocalTimeAndUTC $y $m $d $h $mn $s UTC] {}
	    } elseif { $tail == "Z" } {
		set goodn 1
		# convert to local date
		foreach "y m d h mn s" \
		    [LocalTimeAndUTC $y $m $d $h $mn $s UTC] {}
	    } elseif { $tail == "" } {
		# local time : should convert to UTC 
		set goodn 1 
	    }
	}
    } else {
	set l [split $date -] ; set n [llength $l]
	if { $n == 1 } {
	    # no "-"s
	    set l [split $date /] ; set n [llength $l]
	    if { $n == 1 } {
		# neither "-"s nor "/"s
		set l [split $date .] ; set n [llength $l]
		if { $n == 3 } {
		    # YYYY.MM.DD HH:MM:SS
		    set y [lindex $l 0]
		    set m [string trimleft [lindex $l 1] 0]
		    set n [scan [lindex $l 2] "%0d %0d:%0d:%0d" d h mn s]
		    set goodn 4
		}
	    } elseif { $n == 3 } {
		# MM/DD/YYYY HH:MM:SS
		set m [string trimleft [lindex $l 0] 0]
		set d [string trimleft [lindex $l 1] 0]
		set n [scan [lindex $l 2] "%0d %0d:%0d:%0d" y h mn s]
		set goodn 4
	    }
	} elseif { $n == 3 } {
	    set mt [lindex $l 1]
	    if { [set m [Month $mt]] } {
		# DD-MMM-YYYY HH:MM:SS
		set d [string trimleft [lindex $l 0] 0]
		set n [scan [lindex $l 2] "%d %0d:%0d:%0d" y h mn s]
		set goodn 4
	    } elseif { [scan $mt %0d m] } {
		# YYYY-MM-DD HH:MM:SS
		set y [lindex $l 0]
		set n [scan [lindex $l 2] "%0d %0d:%0d:%0d" d h mn s]
		set goodn 4
	    }
	} elseif { $n == 2 } {
	    # YYYY/MM/DD-HH:MM:SS
	    set n [scan [lindex $l 0] "%0d/%0d/%0d" y m d]
	    incr n [scan [lindex $l 1] "%0d:%0d:%0d" h mn s]
	    set goodn 6
	}
    }
    if { $n == $goodn && [CheckDateEls $y $m $d $h $mn $s] } {
	return [list $y $m $d $h $mn $s]
    }
    return ""
}

proc DateToSecs {y m d h mn s} {
    # convert date to seconds ellapsed since beginning of $YEAR0, not
    #  necessarily a leap year
    global YEAR0

    return [DateToSecsFrom $y $m $d $h $mn $s $YEAR0]
}

proc DateToSecsFrom {y m d h mn s year0} {
    # convert date to seconds ellapsed since beginning of $year0, not
    #  necessarily a leap year
    global DAYSOF

    # make sure there are no leading 0's
    foreach v "y m d h mn s" {
	scan [set $v] %0d $v
    }
    # days in current month and in whole years except February 29ths
    set days [expr 365*($y-$year0)+$d-1]
    set yy [expr $year0+$year0%4]
    while { $yy < $y } {
	if { $yy%100!=0 || $yy%400==0 } { incr days }
	incr yy 4
    }
    # add days in whole months of current year
    if { $m>2 && $y%4==0 && ($y%100!=0 || $y%400==0) } { incr days }
    incr m -1
    while { $m > 0 } { incr days $DAYSOF($m) ; incr m -1 }
    return [expr $days*86400+$h*3600+$mn*60+$s]
}

proc DateFromSecs {secs} {
    # build date from seconds ellapsed since beginning of $YEAR0, not
    #  necessarily a leap year, using default format
    global DateFormat

    return [eval FormatDate $DateFormat [DateIntsFromSecs $secs]]
}

proc DateIntsFromSecs {secs} {
    # compute date from seconds ellapsed since beginning of $YEAR0,
    #  not necessarily a leap year
    # return list with year, month, day, hour, minutes, seconds as integers
    global DAYSOF YEAR0

    set s [expr $secs%60] ; set x [expr ($secs-$s)/60]
    set mn [expr $x%60] ; set x [expr ($x-$mn)/60]
    set h [expr $x%24] ; set x [expr ($x-$h)/24]
    set y [expr int($x/365)] ; set yd [expr $y*365]
    set yy [expr $YEAR0+$YEAR0%4] ; incr y $YEAR0
    while { $yy < $y } {
	if { $yy%100!=0 || $yy%400==0 } { incr yd }
	incr yy 4
    }
    if { $yd > $x } {
	incr y -1
	if { $y%4==0 && ($y%100!=0 || $y%400==0) } { incr yd -1 }
	incr yd -365
    }
    incr x [expr 1-$yd] ; set m 1
    while { $x > $DAYSOF($m) } {
	if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } {
	    if { $x == 29 } { break }
	    incr x -1
	}
	incr x -$DAYSOF($m) ; incr m
    }
    return [list $y $m $x $h $mn $s]
}

proc DateFromSecsFmt {secs fmt} {
    # build date from seconds ellapsed since beginning of $YEAR0, not
    #  necessarily a leap year, using given format
    global DAYSOF YEAR0

    set s [expr $secs%60] ; set x [expr ($secs-$s)/60]
    set mn [expr $x%60] ; set x [expr ($x-$mn)/60]
    set h [expr $x%24] ; set x [expr ($x-$h)/24]
    set y [expr int($x/365)] ; set yd [expr $y*365]
    set yy [expr $YEAR0+$YEAR0%4] ; incr y $YEAR0
    while { $yy < $y } {
	if { $yy%100!=0 || $yy%400==0 } { incr yd }
	incr yy 4
    }
    if { $yd > $x } {
	incr y -1
	if { $y%4==0 && ($y%100!=0 || $y%400==0) } { incr yd -1 }
	incr yd -365
    }
    incr x [expr 1-$yd] ; set m 1
    while { $x > $DAYSOF($m) } {
	if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } {
	    if { $x == 29 } { break }
	    incr x -1
	}
	incr x -$DAYSOF($m) ; incr m
    }
    return [FormatDate $fmt $y $m $x $h $mn $s]
}

proc NextDay {y m d} {
    # return list with year, month and day for the day following the
    #  given date
    global DAYSOF

    if { [incr d] > $DAYSOF($m) } {
	if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } {
	    if { $d != 29 } {
		set d 1 ; set m 3
	    }
	} else {
	    set d 1
	    if { [incr m] == 13 } {
		incr y ; set m 1
	    }
	}
    }
    return [list $y $m $d]
}

proc PreviousDay {y m d} {
    # return list with year, month and day for the day before the
    #  given date
    global DAYSOF

    if { [incr d -1] == 0 } {
	if { [incr m -1] == 0 } {
	    set m 12 ; incr y -1
	}
	if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } {
	    set d 29
	} else {
	    set d $DAYSOF($m)
	}
    }
    return [list $y $m $d]
}

proc Month {name} {
    # from month name to 0 (error) or 1..12
    global ALLMONTH

    for {set m 1} { $m<13 } { incr m } {
	if { [lsearch -exact $ALLMONTH($m) $name] != -1 } { return $m }
    }
    return 0
}

proc LocalTimeAndUTC {y m d h mn s origin} {
    # convert from/to UTC to/from local time
    #  $origin =="UTC"  converts from UTC
    #          =="local"  converts to UTC
    # return list with same elements as in the args
    global TimeOffset DAYSOF

    if { $origin != "UTC" } {
	incr h $TimeOffset
    } else { set h [expr $h-$TimeOffset] }
    if { $h < 0 } {
	incr h 24
	if { [incr d -1] == 0 } {
	    if { [incr m -1] == 0 } {
		set m 12 ; set d 31 ; incr y -1
	    } else {
		if { $m == 2 && $y%4 == 0 && \
			($y%100 != 0 || $y%400 == 0) } {
		    set d 29
		} else {
		    set d $DAYSOF($m)
		}
	    }
	}
    } elseif { $h > 23 } {
	incr h -24
	if { [incr d] > $DAYSOF($m) } {
	    if { $m == 2 } {
		if { $y%4 != 0 || ($y%100 == 0 && $y%400 != 0) } {
		    set m 3 ; set d 1
		}
	    } elseif { $m == 12 } {
		set m 1 ; set d 1 ; incr y
	    } else {
		incr m ; set d 1
	    }
	}
    }
    return [list $y $m $d $h $mn $s]
}

proc CompatibleDates {y0 m0 d0 secs1} {
    # check whether two dates are compatible
    #  $y0, $m0, $d0 give one date (year, month name and day) and may
    #  be * (to match anything)
    #  $secs1 is the other date in seconds since $YEAR0
    # return 1 if compatible

    if { $m0 != "*" && [set m0 [Month $m0]] == 0 } {
	BUG calling CompatibleDates with bad month name
    }
    foreach "y1 m1 d1 h1 mn1 s1" [DateIntsFromSecs $secs1] {}
    if { [string match $d0-$m0-$y0 $d1-$m1-$y1] } { return 1 }
    return 0
}

### varia

proc VectorBearing {vx vy} {
    # compute bearing (angle from y-axis clockwise, in 0..359) of vector
    #  with given components
    # return "_" if components are too small

    if { abs($vy)+abs($vx) > 1e-20 } {
	set b [expr round(atan2($vx,$vy)*57.29577951308232087684)]
	if { $b < 0 } { incr b 360 }
    } else { set b "_" }
    return $b
}