1# 2# This file is part of: 3# 4# gpsman --- GPS Manager: a manager for GPS receiver data 5# 6# Copyright (c) 1998-2013 Miguel Filgueiras migfilg@t-online.de 7# 8# This program is free software; you can redistribute it and/or modify 9# it under the terms of the GNU General Public License as published by 10# the Free Software Foundation; either version 3 of the License, or 11# (at your option) any later version. 12# 13# This program is distributed in the hope that it will be useful, 14# but WITHOUT ANY WARRANTY; without even the implied warranty of 15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16# GNU General Public License for more details. 17# 18# You should have received a copy of the GNU General Public License 19# along with this program. 20# 21# File: compute.tcl 22# Last change: 6 October 2013 23# 24# Includes contributions by 25# - Matt Martin (matt.martin _AT_ ieee.org) marked "MGM contribution" 26# - Valere Robin (valere.robin _AT_ wanadoo.fr) marked "VR contribution" 27# 28 29## Some formulae kindly supplied by 30# Luisa Bastos, Universidade do Porto 31# Gil Goncalves, Universidade de Coimbra 32# Computation of area of spherical polygon adapted from sph_poly.c in 33# "Graphics Gems IV", edited by Paul Heckbert, Academic Press, 1994. 34# Formula for ellipsoid radius from 35# "Ellipsoidal Area Computations of Large Terrestrial Objects" 36# by Hrvoje Lukatela 37# http://www.geodyssey.com/papers/ggelare.html 38# Algorithm for finding a point at given distance and bearing from another 39# taken from the program "forward" available from 40# ftp://www.ngs.noaa.gov/pub/pcsoft/for_inv.3d/ 41# It uses the modified Rainsford's Method with Helmert's elliptical terms, 42# and is effective in any azimuth and at any distance short of antipodal. 43## 44 45### positions and coordinates 46 47proc Coord {pformt coord negh} { 48 # convert coordinate $coord in format $pformt, with negative heading $negh, 49 # to signed degrees 50 # $pformt in {DMS, DMM, DDD, GRA} 51 52 set coord [string trim $coord] 53 set sign 1 54 set h [string index "$coord" 0] 55 if { ! [regexp {[0-9]} $h] } { 56 if { $h == $negh || $h == "-" } { set sign -1 } 57 set coord [string range "$coord" 1 end] 58 } 59 switch $pformt { 60 DMS { 61 scan "$coord" "%d %d %f" d m s 62 return [expr $sign*($d+$m/60.0+$s/3600.0)] 63 } 64 DMM { 65 scan "$coord" "%d %f" d m 66 return [expr $sign*($d+$m/60.0)] 67 } 68 DDD { 69 return [expr $sign*$coord] 70 } 71 GRA { 72 return [expr $sign*0.9*$coord] 73 } 74 } 75} 76 77proc FormatPosition {latd longd datum reqpfmt pfdatum args} { 78 # produce formatted position from $latd and $longd in decimal 79 # degrees for $datum 80 # $reqpfmt is the required position format 81 # $pfdatum is the datum to use with this format 82 # if $pfdatum=="", $datum is used unless $reqpfmt needs a fixed datum 83 # otherwise $pfdatum must be compatible with $reqpfmt; if it is 84 # not a call to BUG will be made 85 # $args is either "" or is an alternative format to be used with $datum 86 # when conversion to $reqpfmt fails (out of grid scope) 87 # result is a list with: position representation, position format 88 # used and datum used 89 # if $args=="" and conversion to $reqpfmt yields an out of grid scope 90 # result ("--" in the zone or x-field), this result is used in 91 # building the position representation 92 # a position representation is a list whose first two elements 93 # are lat and long in degrees, and whose rest is dependent on format 94 # type: 95 # ==latlong: lat and long in external format 96 # ==utm: zones east and north and x y coordinates 97 # ==grid: zone name and x y coordinates 98 # ==nzgrid: x y coordinates 99 # ==mh: Maidenhead-locator (6 characters) 100 global POSTYPE 101 102 set la $latd ; set lo $longd 103 if { $pfdatum == "" } { 104 set udatum $datum 105 } else { 106 set udatum $pfdatum 107 if { $pfdatum != $datum } { 108 foreach "la lo" [ToDatum $la $lo $datum $pfdatum] { break } 109 } 110 } 111 switch [set ptype $POSTYPE($reqpfmt)] { 112 latlong { 113 set p [FormatLatLong $la $lo $reqpfmt] 114 } 115 utm { 116 set p [linsert [DegreesToUTM $la $lo $udatum] 0 $la $lo] 117 } 118 grid - nzgrid { 119 if { [set gd [GridDatum $reqpfmt $udatum]] != $udatum } { 120 if { $pfdatum != "" } { 121 BUG calling FormatPosition with wrong pfdatum 122 } 123 set udatum $gd 124 foreach "la lo" [ToDatum $latd $longd $datum $gd] { break } 125 } 126 if { $ptype == "grid" } { 127 set p [DegreesTo$reqpfmt $la $lo $udatum] 128 } else { set p [DegreesToNZGrid $reqpfmt $la $lo $udatum] } 129 if { [lindex $p 0] == "--" && $args != "" } { 130 # out of scope 131 return [FormatPosition $latd $longd $datum $args ""] 132 } 133 set p [linsert $p 0 $la $lo] 134 } 135 mh { 136 set p [list $la $lo [DegreesToMHLoc $la $lo]] 137 } 138 } 139 return [list $p $reqpfmt $udatum] 140} 141 142proc FormatLatLong {latd longd pfmt} { 143 # produce formatted position from $latd and $longd in decimal 144 # degrees in latitude/longitude position format 145 # result is the position representation: a list whose first two elements 146 # are lat and long in degrees, followed by lat and long in external format 147 global COUTFMT 148 149 if { $latd < 0 } { 150 # VR contribution: 0-$latd instead of -$latd; same for $longd 151 set lax [expr 0-$latd] ; set hlat S 152 } else { set lax $latd ; set hlat N } 153 if { $longd < 0 } { 154 set lox [expr 0-$longd] ; set hlng W 155 } else { set lox $longd ; set hlng E } 156 if { $pfmt == "GRA" } { 157 set lax [format $COUTFMT(GRA) [expr $latd/0.9]] 158 set lox [format $COUTFMT(GRA) [expr $longd/0.9]] 159 return [list $latd $longd $lax $lox] 160 } 161 return [list $latd $longd "$hlat[ExtDegrees $pfmt $lax]" \ 162 "$hlng[ExtDegrees $pfmt $lox]"] 163} 164 165proc ExtDegrees {pformt degs} { 166 # from signed degrees $degs to external representation of format $pformt 167 # $pformt in {DMS, DMM, DDD, DMSsimpl} 168 # DMSsimpl is similar to DMS, but for values less than 1 degree the 169 # representation is either MM'SS.S" or SS.S" 170 global COUTFMT 171 172 switch -glob $pformt { 173 DMS* { 174 if { $degs < 0 } { 175 set degs [expr -$degs] ; set sign -1 176 } else { set sign 1 } 177 set d [expr int($degs)] 178 set degs [expr ($degs-$d)*60] 179 set d [expr $sign*$d] 180 set m [expr int($degs)] 181 set s [expr ($degs-$m)*60] 182 if { $s > 59.95 } { set s 0 ; incr m } 183 if { $m > 59 } { set m 0 ; incr d } 184 set sf $COUTFMT(sec) 185 if { $d > 0 || $pformt == "DMS" } { 186 return [format "%d %02d $sf" $d $m $s] 187 } 188 if { $m > 0 } { 189 return [format "%02d'$sf" $m $s] 190 } 191 return [format "$sf\"" $s] 192 } 193 DMM { 194 if { $degs < 0 } { 195 set degs [expr -$degs] ; set sign -1 196 } else { set sign 1 } 197 set d [expr int($degs)] 198 set degs [expr ($degs-$d)*60] 199 set d [expr $sign*$d] 200 if { $degs > 59.995 } { set degs 0 ; incr d } 201 return [format "%d $COUTFMT(min)" $d $degs] 202 } 203 DDD { 204 return [format $COUTFMT(deg) $degs] 205 } 206 } 207} 208 209proc ChangeCoordSign {pos} { 210 211 if { $pos != "N" } { return W } 212 return S 213} 214 215proc ConvertPos {pfrmt zone x y datum npfrmt} { 216 # create position representation in format $npfrmt from position 217 # $x, $y (as decimal numbers) in format $pfrmt 218 # if $pfrmt in {DDD, GRA} take $x as long and $y as lat 219 # $zone is either a list with ze and zn for UTM/UPS, or the zone 220 # name for other grids (and is meaningless for other formats 221 # values of $x, $y are checked, as well as suitability of $datum 222 # and of validity of new format 223 # return -1 on error 224 global POSTYPE MESS TXT COUTFMT 225 226 switch [set ptype $POSTYPE($pfrmt)] { 227 latlong { 228 # must be either DDD or GRA 229 if { $pfrmt == "GRA" } { 230 set x [expr $x*0.9] ; set y [expr $y*0.9] 231 } 232 if { $y < 0 } { 233 set la [expr -$y] ; set hlat S 234 } else { 235 set la $y ; set hlat N 236 } 237 if { $x < 0 } { 238 set lo [expr -$x] ; set hlng W 239 } else { 240 set lo $x ; set hlng E 241 } 242 if { ! [CheckCoord GMMessage $pfrmt $la N 90] || \ 243 ! [CheckCoord GMMessage $pfrmt $lo E 180] } { 244 return -1 245 } 246 if { $npfrmt == "DDD" } { 247 return [list $y $x "${hlat}[format $COUTFMT(deg) $la]" \ 248 "${hlng}[format $COUTFMT(deg) $lo]"] 249 } 250 set la $y ; set lo $x 251 } 252 utm { 253 foreach "ze zn" $zone {} 254 foreach "la lo" [UTMToDegrees $ze $zn $x $y $datum] { break } 255 if { $npfrmt == $pfrmt } { 256 return [list $la $lo $ze $zn $x $y] 257 } 258 } 259 grid - nzgrid { 260 if { [BadDatumFor $pfrmt $datum GMMessage] != 0 } { 261 return -1 262 } 263 if { [set p [GridToDegrees $pfrmt $zone $x $y $datum]] == 0 } { 264 GMMessage $MESS(outofgrid) 265 return -1 266 } 267 if { $npfrmt == $pfrmt } { 268 if { $ptype == "grid" } { 269 return [linsert $p end $zone $x $y] 270 } else { return [linsert $p end $x $y] } 271 } 272 } 273 default { 274 BUG calling ConvertPos with wrong format 275 } 276 } 277 # $datum is suitable for $npfrmt 278 set p [lindex [FormatPosition $la $lo $datum $npfrmt $datum] 0] 279 if { [lindex $p 2] == "--" } { 280 GMMessage $MESS(outofgrid) 281 return -1 282 } 283 return $p 284} 285 286### distances and geographic bearings 287# formulae from Spherical Trignometry 288# the Law of Cosines, used below to compute distances and bearings, is known 289# to be unsuitable for computations because of excessive round-off errors; 290# alternative procedures are supplied in file acccomp.tcl and are loaded 291# if $ACCFORMULAE (accurate formulae option) is set 292 293proc ComputeDist {p1 p2 datum} { 294 # distance (km) between positions $p1 and $p2 with same datum 295 # formulae kindly supplied by Luisa Bastos, Universidade do Porto 296 # and Gil Goncalves, Universidade de Coimbra 297 298 set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1] 299 set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1] 300 if { $lad1==$lad2 && $lod1==$lod2 } { return 0 } 301 set la1 [expr $lad1*0.01745329251994329576] 302 set lo1 [expr $lod1*0.01745329251994329576] 303 set la2 [expr $lad2*0.01745329251994329576] 304 set lo2 [expr $lod2*0.01745329251994329576] 305 set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)] 306 if { $x >= 1 } { return 0 } 307 return [expr 1e-3*[lindex [EllipsdData $datum] 0]*acos($x)] 308} 309 310proc ComputeBear {p1 p2 datum} { 311 # bearing from position $p1 to $p2 with same datum 312 # $datum not used here but kept for compatibility with replacement 313 # proc in acccomp.tcl 314 # formulae kindly supplied by Luisa Bastos, Universidade do Porto 315 # and Gil Goncalves, Universidade de Coimbra 316 317 set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1] 318 set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1] 319 if { $lad1==$lad2 && $lod1==$lod2 } { return 0 } 320 set la1 [expr $lad1*0.01745329251994329576] 321 set lo1 [expr $lod1*0.01745329251994329576] 322 set la2 [expr $lad2*0.01745329251994329576] 323 set lo2 [expr $lod2*0.01745329251994329576] 324 set da [expr $la2-$la1] ; set do [expr $lo2-$lo1] 325 if { [expr abs($da)] < 1e-20 } { 326 if { [expr abs($do)] < 1e-20 } { 327 set b 0 328 } elseif { $do < 0 } { 329 set b 270 330 } else { set b 90 } 331 } elseif { [expr abs($do)] < 1e-20 } { 332 if { $da < 0 } { 333 set b 180 334 } else { set b 0 } 335 } else { 336 set b [expr round(atan2(sin($do), \ 337 tan($la2)*cos($la1)-sin($la1)*cos($do)) \ 338 *57.29577951308232087684)] 339 if { $b < 0 } { 340 if { $do < 0 } { incr b 360 } else { incr b 180 } 341 } elseif { $do < 0 } { incr b 180 } 342 } 343 return $b 344} 345 346proc ComputeDistBear {p1 p2 datum} { 347 # distance between and bearing from positions $p1 and $p2 with same datum 348 # formulae kindly supplied by Luisa Bastos, Universidade do Porto 349 # and Gil Goncalves, Universidade de Coimbra 350 351 set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1] 352 set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1] 353 if { $lad1==$lad2 && $lod1==$lod2 } { return "0 0" } 354 set la1 [expr $lad1*0.01745329251994329576] 355 set lo1 [expr $lod1*0.01745329251994329576] 356 set la2 [expr $lad2*0.01745329251994329576] 357 set lo2 [expr $lod2*0.01745329251994329576] 358 set a [lindex [EllipsdData $datum] 0] 359 # distance (km) 360 set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)] 361 if { $x < 1 } { 362 set d [expr 1e-3*$a*acos($x)] 363 } else { set d 0 } 364 # bearing 365 set da [expr $la2-$la1] ; set do [expr $lo2-$lo1] 366 if { [expr abs($da)] < 1e-20 } { 367 if { [expr abs($do)] < 1e-20 } { 368 set b 0 369 } elseif { $do < 0 } { 370 set b 270 371 } else { set b 90 } 372 } elseif { [expr abs($do)] < 1e-20 } { 373 if { $da < 0 } { 374 set b 180 375 } else { set b 0 } 376 } else { 377 set b [expr round(atan2(sin($do), \ 378 tan($la2)*cos($la1)-sin($la1)*cos($do)) \ 379 *57.29577951308232087684)] 380 if { $b < 0 } { 381 if { $do < 0 } { incr b 360 } else { incr b 180 } 382 } elseif { $do < 0 } { incr b 180 } 383 } 384 return [list $d $b] 385} 386 387proc SetDatumData {datum} { 388 # this proc sets datum parameters as global variables for repeated use 389 # in conversions; see procs ComputeDistFD, ComputeDistBearFD 390 global DatumA DatumF FDDatum 391 392 if { $FDDatum != $datum } { 393 set dt [EllipsdData $datum] 394 set DatumA [lindex $dt 0] ; set DatumF [lindex $dt 1] 395 set FDDatum $datum 396 } 397 return 398} 399 400proc ComputeDistFD {p1 p2} { 401 # compute distance (km) between positions $p1 and $p2 assuming datum 402 # parameters where set by calling SetDatumData 403 global DatumA 404 # formulae kindly supplied by Luisa Bastos, Universidade do Porto 405 # and Gil Goncalves, Universidade de Coimbra 406 407 set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1] 408 set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1] 409 if { $lad1==$lad2 && $lod1==$lod2 } { return 0 } 410 set la1 [expr $lad1*0.01745329251994329576] 411 set lo1 [expr $lod1*0.01745329251994329576] 412 set la2 [expr $lad2*0.01745329251994329576] 413 set lo2 [expr $lod2*0.01745329251994329576] 414 set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)] 415 if { $x >= 1 } { return 0 } 416 return [expr 1e-3*$DatumA*acos($x)] 417} 418 419proc ComputeDistBearFD {p1 p2} { 420 # compute distance (km) between and bearing from positions $p1 and $p2 421 # assuming datum parameters where set by calling SetDatumData 422 global DatumA 423 # formulae kindly supplied by Luisa Bastos, Universidade do Porto 424 # and Gil Goncalves, Universidade de Coimbra 425 426 set lad1 [lindex $p1 0] ; set lod1 [lindex $p1 1] 427 set lad2 [lindex $p2 0] ; set lod2 [lindex $p2 1] 428 if { $lad1==$lad2 && $lod1==$lod2 } { return "0 0" } 429 set la1 [expr $lad1*0.01745329251994329576] 430 set lo1 [expr $lod1*0.01745329251994329576] 431 set la2 [expr $lad2*0.01745329251994329576] 432 set lo2 [expr $lod2*0.01745329251994329576] 433 # distance (km) 434 set x [expr cos($lo1-$lo2)*cos($la1)*cos($la2)+sin($la1)*sin($la2)] 435 if { $x < 1 } { 436 set d [expr 1e-3*$DatumA*acos($x)] 437 } else { set d 0 } 438 # bearing 439 set da [expr $la2-$la1] ; set do [expr $lo2-$lo1] 440 if { [expr abs($da)] < 1e-20 } { 441 if { [expr abs($do)] < 1e-20 } { 442 set b 0 443 } elseif { $do < 0 } { 444 set b 270 445 } else { set b 90 } 446 } elseif { [expr abs($do)] < 1e-20 } { 447 if { $da < 0 } { 448 set b 180 449 } else { set b 0 } 450 } else { 451 set b [expr round(atan2(sin($do), \ 452 tan($la2)*cos($la1)-sin($la1)*cos($do)) \ 453 *57.29577951308232087684)] 454 if { $b < 0 } { 455 if { $do < 0 } { incr b 360 } else { incr b 180 } 456 } elseif { $do < 0 } { incr b 180 } 457 } 458 return [list $d $b] 459} 460 461proc CompWPDistBear {wp1 wp2} { 462 # distance (km) between and bearing from WPs with names $wp1, $wp2 463 global WPPFrmt WPPosn WPDatum 464 465 set i1 [IndexNamed WP $wp1] 466 set i2 [IndexNamed WP $wp2] 467 if { $i1<0 || $i2<0 } { return "--- ---" } 468 if { $i1 == $i2 } { return "0 0" } 469 return [CompDistBearDatums $WPPosn($i1) $WPDatum($i1) \ 470 $WPPosn($i2) $WPDatum($i2)] 471} 472 473proc CompDistBearDatums {p1 d1 p2 d2} { 474 # distance (km) between and bearing from positions $p1 and $p2 475 # with datums $d1 and $d2 476 if { $d1 != $d2 } { 477 set p2 [ToDatum [lindex $p2 0] [lindex $p2 1] $d2 $d1] 478 } 479 return [ComputeDistBear $p1 $p2 $d1] 480} 481 482## point at given distance and bearing from another point 483 484proc CoordsAtDistBear {pos dist bear datum} { 485 # compute coordinates of point at given distance (km) and bearing 486 # (0-360 degrees) from another point 487 # algorithm taken from the program "forward" available from 488 # ftp://www.ngs.noaa.gov/pub/pcsoft/for_inv.3d/ 489 # It uses the modified Rainsford's Method with Helmert's elliptical 490 # terms, and is effective in any azimuth and at any distance short of 491 # antipodal. 492 # $dist is constrained here to be less than pi*semi-major axis 493 global MESS 494 495 foreach "lat0 long0" $pos { break } 496 foreach "a f" [EllipsdData $datum] { break } 497 if { [set dist [expr $dist*1000]] >= 3.1415926*$a } { 498 GMMessage $MESS(distlarge) 499 return [list $lat $long] 500 } 501 set phi [expr $lat0*0.01745329251994329576] 502 set lam [expr $long0*0.01745329251994329576] 503 set azm [expr $bear*0.01745329251994329576] 504 set r [expr 1-$f] 505 set tu [expr $r*sin($phi)/cos($phi)] 506 set sf [expr sin($azm)] ; set cf [expr cos($azm)] 507 set baz 0 508 if { [set cf [expr cos($azm)]] != 0 } { 509 set baz [expr atan2($tu,$cf)*2] 510 } 511 set cu [expr 1.0/sqrt($tu*$tu+1)] 512 set su [expr $tu*$cu] 513 set sa [expr $cu*$sf] 514 set c2a [expr -$sa*$sa+1] 515 set x [expr sqrt((1.0/$r/$r-1)*$c2a+1)+1] 516 set x [expr ($x-2.0)/$x] 517 set c [expr 1-$x] 518 set c [expr ($x*$x/4.0+1.0)/$c] 519 set d [expr (0.375*$x*$x-1)*$x] 520 set y [set tu [expr 1.0*$dist/$r/$a/$c]] 521 while 1 { 522 set sy [expr sin($y)] ; set cy [expr cos($y)] 523 set cz [expr cos($baz+$y)] 524 set e [expr $cz*$cz*2-1] 525 set c $y 526 set x [expr $e*$cy] 527 set y [expr ((($sy*$sy*4-3)*($e+$e-1)*$cz*$d/6.0+$x)*$d/4.0-$cz)* \ 528 $sy*$d+$tu] 529 if { abs($y-$c) <= 5e-14 } { break } 530 } 531 set baz [expr $cu*$cy*$cf-$su*$sy] 532 set c [expr $r*sqrt($sa*$sa+$baz*$baz)] 533 set d [expr $su*$cy+$cu*$sy*$cf] 534 set phi2 [expr atan2($d,$c)] 535 set c [expr $cu*$cy-$su*$sy*$cf] 536 set x [expr atan2($sy*$sf,$c)] 537 set c [expr ((-3*$c2a+4)*$f+4)*$c2a*$f/16.0] 538 set d [expr (($e*$cy*$c+$cz)*$sy*$c+$y)*$sa] 539 set lam2 [expr $lam+$x-(1-$c)*$d*$f] 540 # set baz [expr atan2($sa,$baz)+3.14159265358979323846] 541 542 set lat1 [expr $phi2*57.29577951308232087684] 543 set long1 [expr $lam2*57.29577951308232087684] 544 return [list $lat1 $long1] 545} 546 547### area 548 549proc ProjectedArea {wpixs} { 550 # compute area of polygon whose boundary is the polyline formed by the 551 # WPs in given list 552 # polygon cannot be self-intersecting (no test on this!) 553 # return value in square km 554 global WPPosn WPDatum ASKPROJPARAMS 555 556 set ix [lindex $wpixs 0] 557 set datum $WPDatum($ix) ; set p $WPPosn($ix) 558 set ps [list [list [set lat [lindex $p 0]] [lindex $p 1] $datum]] 559 foreach ix [lreplace $wpixs 0 0] { 560 set dt $WPDatum($ix) ; set p $WPPosn($ix) 561 if { $dt != $datum } { 562 set p [ToDatum [lindex $p 0] [lindex $p 1] $dt $datum] 563 } 564 lappend ps [list [lindex $p 0] [lindex $p 1] $datum] 565 } 566 if { $lat>=-80 && $lat<=84 } { 567 set proj TM 568 } else { 569 # in fact this will use UPS 570 set proj UTM 571 } 572 set ask $ASKPROJPARAMS 573 set ASKPROJPARAMS 0 574 set xy [ProjInit $proj AProj $datum $ps] 575 set ASKPROJPARAMS $ask 576 set xmin [set x0 [lindex $xy 0]] ; set ymin [set y0 [lindex $xy 1]] 577 set xs "" ; set ys "" 578 foreach p [lreplace $ps 0 0] { 579 set xy [eval Proj${proj}Point AProj $p] 580 set x [lindex $xy 0] ; set y [lindex $xy 1] 581 lappend xs $x ; lappend ys $y 582 if { $x < $xmin } { set xmin $x } 583 if { $y < $ymin } { set ymin $y } 584 } 585 set sum 0 586 set x0 [set x [expr $x0-$xmin]] ; set y0 [set y [expr $y0-$ymin]] 587 foreach x1 $xs y1 $ys { 588 set x1 [expr $x1-$xmin] ; set y1 [expr $y1-$ymin] 589 set sum [expr $sum+$x*$y1-$x1*$y] 590 set x $x1 ; set y $y1 591 } 592 set sum [expr $sum+$x*$y0-$x0*$y] 593 return [expr abs(0.5e-6*$sum)] 594} 595 596proc SphericalArea {wpixs} { 597 # compute area of polygon whose boundary is the polyline formed by the 598 # WPs in given list 599 # computation is based on a spherical approximation 600 # return negative value if there are precision errors 601 # otherwise value in square km 602 global WPPosn WPDatum 603 604 set ix [lindex $wpixs 0] 605 set datum $WPDatum($ix) ; set p [set p0 $WPPosn($ix)] 606 set phi0 [set phi [expr [lindex $p 0]*0.01745329251994329576]] 607 set lam0 [set lam [expr [lindex $p 1]*0.01745329251994329576]] 608 set c0 [expr cos($phi0)] 609 # computation of mean radius of the ellipsoid at a point 610 # as described in Hrvoje Lukatela, "Ellipsoidal Area Computations of 611 # Large Terrestrial Objects", http://www.geodyssey.com/papers/ggelare.html 612 set dt [EllipsdData $datum] 613 set a [lindex $dt 0] ; set b [expr $a*(1-[lindex $dt 1])] 614 set s [expr sin($phi)] 615 set a2 [expr $a*$a] 616 # radius in km 617 set r [expr 1e-3*$a2*$b/($a2*$c0*$c0+$b*$b*$s*$s)] 618 # computation of area adapted from sph_poly.c in 619 # "Graphics Gems IV", edited by Paul Heckbert, Academic Press, 1994. 620 set srad 0 621 SetDatumData $datum 622 foreach ixn [lreplace $wpixs 0 0] { 623 set pos $WPPosn($ixn) 624 if { $WPDatum($ixn) != $datum } { 625 set pos [ToDatum [lindex $pos 0] [lindex $pos 1] \ 626 $WPDatum($ixn) $datum] 627 } 628 set phi1 [expr [lindex $pos 0]*0.01745329251994329576] 629 set lam1 [expr [lindex $pos 1]*0.01745329251994329576] 630 set c1 [expr cos($phi1)] 631 if { $lam0 != $lam1 } { 632 set HavA [expr ((1-cos($phi1-$phi0))/2.0)+ \ 633 ((1-cos($lam1-$lam0))/2.0)*$c0*$c1] 634 set A [expr 2*asin(sqrt($HavA))] 635 set B [expr 1.5707963267948966192313-$phi1] 636 set C [expr 1.5707963267948966192313-$phi0] 637 set S [expr 0.5*($A+$B+$C)] 638 set T [expr tan($S/2.0)*tan(($S-$A)/2.0)* \ 639 tan(($S-$B)/2.0)*tan(($S-$C)/2.0)] 640 if { abs($T) < 1e-8 } { return -1 } 641 set E [expr abs(4*atan(sqrt(abs($T))))] 642 if { $lam1 < $lam0 } { set E [expr -$E] } 643 set srad [expr $srad+$E] 644 } 645 set phi0 $phi1 ; set lam0 $lam1 ; set c0 $c1 646 } 647 set lam1 $lam 648 if { $lam0 != $lam1 } { 649 set phi1 $phi ; set c1 [expr cos($lam1)] 650 set HavA [expr ((1-cos($phi1-$phi0))/2.0)+ \ 651 ((1-cos($lam1-$lam0))/2.0)*$c0*$c1] 652 set A [expr 2*asin(sqrt($HavA))] 653 set B [expr 1.5707963267948966192313-$phi1] 654 set C [expr 1.5707963267948966192313-$phi0] 655 set S [expr 0.5*($A+$B+$C)] 656 set T [expr tan($S/2.0)*tan(($S-$A)/2.0)* \ 657 tan(($S-$B)/2.0)*tan(($S-$C)/2.0)] 658 if { abs($T) < 1e-8 } { return -1 } 659 set E [expr abs(4*atan(sqrt(abs($T))))] 660 if { $lam1 < $lam0 } { set E [expr -$E] } 661 set srad [expr $srad+$E] 662 } 663 return [expr abs($r*$r*$srad)] 664} 665 666### datums 667 668proc ChangeTPsDatum {tps datum1 datum2} { 669 # convert position of TR points on list $tps from $datum1 to $datum2 670 # change only the first 4 elements (position) of each TP representation 671 # this depends on DataIndex (see setup.tcl, Storage), not used here 672 # for speed 673 674 if { $datum1 == $datum2 } { return $tps } 675 set l "" 676 foreach tp $tps { 677 set np [lindex [FormatPosition [lindex $tp 0] [lindex $tp 1] \ 678 $datum1 DMS $datum2] 0] 679 lappend l [concat $np [lrange $tp 4 end]] 680 } 681 return $l 682} 683 684proc ChangeLPsDatum {lps odatum datum pformt} { 685 # change datum of LN list of points 686 # change only the first element (position) of each LP representation 687 # this depends on DataIndex (see setup.tcl, Storage), not used here 688 # for speed 689 # return -1 if the current position format requires a fixed datum 690 # different from $datum 691 global MESS POSTYPE 692 693 if { $odatum == $datum } { return $lps } 694 if { [set gd [BadDatumFor $pformt $datum Ignore]] != 0 } { 695 GMMessage [format $MESS(cantchggriddatum) $gd] 696 return -1 697 } 698 set l "" 699 foreach lp $lps { 700 foreach "latd longd" [lindex $lp 0] { break } 701 lappend l [lreplace $lp 0 0 \ 702 [lindex [FormatPosition $latd $longd $odatum \ 703 $pformt $datum] 0]] 704 } 705 return $l 706} 707 708### altitude 709 710proc UserAltitude {altlst} { 711 # return altitude in user units from internal representation 712 global AltUnit ALSCALE 713 714 if { $altlst == "" } { return "" } 715 if { $AltUnit == "M" } { return [lindex $altlst 0] } 716 if { $AltUnit == [lindex $altlst 2] } { return [lindex $altlst 1] } 717 return [expr [lindex $altlst 0]/$ALSCALE] 718} 719 720proc AltitudeList {a} { 721 # check value for altitude and return internal representation as list 722 # return nil on error 723 global AltUnit ALSCALE 724 725 if { $a == "" } { return "" } 726 if { ! [CheckSignedFloat Ignore $a] } { return nil } 727 if { $AltUnit == "M" } { return $a } 728 return [list [expr $a*$ALSCALE] $a $AltUnit] 729} 730 731### dates 732 733proc Today {dformt} { 734 # build representation of current date under format $dformt 735 736 set dt [clock format [clock seconds] -format "%Y %m %d %H %M %S"] 737 scan $dt %d%0d%0d%0d%0d%0d y m d h mn s 738 return [FormatDate $dformt $y $m $d $h $mn $s] 739} 740 741proc TodayUTC {dformt} { 742 # VR contribution 743 # build representation of current date under format $dformt 744 745 set dt [clock format [clock seconds] -format "%Y %m %d %H %M %S"] 746 scan $dt %d%0d%0d%0d%0d%0d y m d h mn s 747 return [eval FormatDate $dformt [LocalTimeAndUTC $y $m $d $h $mn $s Local]] 748} 749 750proc FormatDay {dformt y m d} { 751 # build representation of date (without time of day) under format $dformt 752 # (see proc FormatDate for date with time of day) 753 # $dformt in 754 # {YYYYMMDD, YYYY-MM-DD, MMDDYYYY, DDMMMYYYY, YYYY/MM/DD, ISO8601} 755 # changes in proc FormatDate are likely to affect this proc! 756 757 set cdate [FormatDate $dformt $y $m $d 0 0 0] 758 if { $dformt == "YYYY/MM/DD" } { 759 return [string range $cdate 0 9] 760 } 761 return [lindex $cdate 0] 762} 763 764proc FormatDate {dformt y m d h mn s} { 765 # build representation of date (including time of day) under format $dformt 766 # (see proc FormatDay for date without time of day) 767 # $dformt either in 768 # {YYYYMMDD, YYYY-MM-DD, MMDDYYYY, DDMMMYYYY, YYYY/MM/DD, ISO8601} 769 # or DDMMYY in which case the date will have no time of day 770 # changes here must be reflected in proc FormatDay and in proc ScanDate 771 # when adding new formats here, array DATEW must be updated in main.tcl 772 global MONTHNAMES 773 774 # make sure there are no leading 0's 775 foreach v "y m d h mn s" { 776 scan [set $v] %0d $v 777 } 778 set h [format "%02d:%02d:%02d" $h $mn $s] 779 switch $dformt { 780 DDMMYY { 781 # MGM contribution 782 return [format "%02d%02d%02d" $d $m [expr $y % 100]] 783 } 784 ISO8601 { 785 # VR contribution 786 return [format "%4d-%02d-%02dT%sZ" $y $m $d $h] 787 } 788 YYYYMMDD { return [format "%4d.%02d.%02d %s" $y $m $d $h] } 789 YYYY-MM-DD { return [format "%4d-%02d-%02d %s" $y $m $d $h] } 790 MMDDYYYY { return [format "%02d/%02d/%4d %s" $m $d $y $h] } 791 DDMMMYYYY { 792 incr m -1 ; set m [lindex $MONTHNAMES $m] 793 return [format "%02d-%s-%4d %s" $d $m $y $h] 794 } 795 YYYY/MM/DD { return [format "%4d/%02d/%02d-%s" $y $m $d $h] } 796 } 797} 798 799proc FormatTime {secs args} { 800 # build represention of seconds in hours:minutes:seconds format 801 # $secs may be a float or an integer 802 # if $args!="" it gives the number of decimal places for the seconds field 803 # that defaults to 0 with $secs rounded to integer 804 805 if { $args != "" && $args > 0 } { 806 set fss "%02.${args}f" ; set fss0 "%8.${args}f" 807 } else { 808 set secs [expr round($secs)] 809 set fss "%02d" ; set fss0 "%8d" 810 } 811 set s0 [expr int($secs)] 812 set s [expr $s0%60] ; set x [expr ($s0-$s)/60] 813 set s [expr $s+$secs-$s0] 814 set mn [expr $x%60] ; set h [expr ($x-$mn)/60] 815 if { $h > 0 } { return [format "%2d:%02d:$fss" $h $mn $s] } 816 if { $mn > 0 } { return [format "%5d:$fss" $mn $s] } 817 return [format $fss0 $s] 818} 819 820proc TimeToSecs {time} { 821 # convert represention of seconds in hours:minutes:seconds format 822 # to seconds (hours and minutes fields may be missing) 823 # return integer or float according to seconds field 824 825 set n [llength [set l [split $time ":"]]] 826 set s [lindex $l end] 827 if { [string first 0 $s] == 0 } { set s [string range $s 1 end] } 828 if { $n == 1 } { return $s } 829 set ll "" 830 foreach f $l { 831 scan $f %0d f 832 lappend ll $f 833 } 834 scan [lindex $l 0] %0d f 835 if { $n == 2 } { 836 return [expr $s+60*$f] 837 } 838 scan [lindex $l 1] %0d g 839 return [expr $s+60*$g+3600*$f] 840} 841 842proc NowTZ {} { 843 # current date under default format with time zone appended 844 global DateFormat 845 846 return "[Today $DateFormat] ([clock format 0 -format %Z])" 847} 848 849proc Now {} { 850 # current date under default format 851 global DateFormat 852 853 return [Today $DateFormat] 854} 855 856array set DAYSOF { 857 1 31 2 28 3 31 4 30 5 31 6 30 7 31 8 31 9 30 858 10 31 11 30 12 31 859} 860 861proc ScanDate {date} { 862 # convert date string to a list with year, month (in 1..12), day, 863 # hours, minutes and seconds 864 # if there is time zone information, return local time 865 # return empty list on error 866 # see proc FormatDate for possible date formats 867 # changes in proc FormatDate are likely to affect this proc! 868 869 set date [string trim $date] 870 if { [regexp {^([0-9][0-9])([0-9][0-9])([0-9][0-9])$} $date \ 871 x d m y] } { 872 # DDMMYY, no time of day 873 # YY taken as 2000+YY if YY<70, and 1900+YY otherwise 874 foreach v "d m y" { 875 scan [set $v] %02d $v 876 } 877 if { $y < 70 } { 878 incr y 2000 879 } else { incr y 1900 } 880 return [list $y $m $d 0 0 0] 881 } 882 set goodn -1 883 # includes VR contribution: decimal comma, empty tail 884 if { [regexp {^ *([-0-9]+)T([0-9:]+([.,][0-9]+)?)([^ ]*) *$} \ 885 $date m ymd hms m tail] } { 886 # coping with ISO8601 format as in: 887 # 2000-01-01T12:00:00Z 888 # 2000-01-01T12:00:00.00000Z 889 # 2004-06-07T13:51:27.2515650-07:00 890 # 2009-09-23T07:02:46.773 891 # 2009-09-23T07:02:46,773 892 # must replace decimal comma by period 893 regsub {,} $hms "." hms 894 set n 1 895 if { [scan $hms %02d:%02d:%0f h mn s] == 3 && \ 896 [scan $ymd %4d-%02d-%02d y m d] == 3 } { 897 if { [set s [expr round($s)]] > 59 } { 898 incr s -60 899 if { [incr mn] > 59 } { 900 incr mn -60 901 if { [incr h] > 23 } { 902 incr h -24 903 foreach "y m d" [NextDay $y $m $d] {} 904 } 905 } 906 } 907 if { [scan $tail %03d:%02d hz mz] == 2 } { 908 set goodn 1 909 # compute UTC date 910 if { $hz > 0 } { set mz [expr -$mz] } 911 set hz [expr -$hz] 912 if { [incr mn $mz] > 59 } { 913 incr h ; incr mn -60 914 } elseif { $mn < 0 } { 915 incr h -1 ; incr mn 60 916 } 917 if { [incr h $hz] > 23 } { 918 incr h -24 919 foreach "y m d" [NextDay $y $m $d] {} 920 } elseif { $h < 0 } { 921 incr h 24 922 foreach "y m d" [PreviousDay $y $m $d] {} 923 } 924 # convert to local date 925 foreach "y m d h mn s" \ 926 [LocalTimeAndUTC $y $m $d $h $mn $s UTC] {} 927 } elseif { $tail == "Z" } { 928 set goodn 1 929 # convert to local date 930 foreach "y m d h mn s" \ 931 [LocalTimeAndUTC $y $m $d $h $mn $s UTC] {} 932 } elseif { $tail == "" } { 933 # local time : should convert to UTC 934 set goodn 1 935 } 936 } 937 } else { 938 set l [split $date -] ; set n [llength $l] 939 if { $n == 1 } { 940 # no "-"s 941 set l [split $date /] ; set n [llength $l] 942 if { $n == 1 } { 943 # neither "-"s nor "/"s 944 set l [split $date .] ; set n [llength $l] 945 if { $n == 3 } { 946 # YYYY.MM.DD HH:MM:SS 947 set y [lindex $l 0] 948 set m [string trimleft [lindex $l 1] 0] 949 set n [scan [lindex $l 2] "%0d %0d:%0d:%0d" d h mn s] 950 set goodn 4 951 } 952 } elseif { $n == 3 } { 953 # MM/DD/YYYY HH:MM:SS 954 set m [string trimleft [lindex $l 0] 0] 955 set d [string trimleft [lindex $l 1] 0] 956 set n [scan [lindex $l 2] "%0d %0d:%0d:%0d" y h mn s] 957 set goodn 4 958 } 959 } elseif { $n == 3 } { 960 set mt [lindex $l 1] 961 if { [set m [Month $mt]] } { 962 # DD-MMM-YYYY HH:MM:SS 963 set d [string trimleft [lindex $l 0] 0] 964 set n [scan [lindex $l 2] "%d %0d:%0d:%0d" y h mn s] 965 set goodn 4 966 } elseif { [scan $mt %0d m] } { 967 # YYYY-MM-DD HH:MM:SS 968 set y [lindex $l 0] 969 set n [scan [lindex $l 2] "%0d %0d:%0d:%0d" d h mn s] 970 set goodn 4 971 } 972 } elseif { $n == 2 } { 973 # YYYY/MM/DD-HH:MM:SS 974 set n [scan [lindex $l 0] "%0d/%0d/%0d" y m d] 975 incr n [scan [lindex $l 1] "%0d:%0d:%0d" h mn s] 976 set goodn 6 977 } 978 } 979 if { $n == $goodn && [CheckDateEls $y $m $d $h $mn $s] } { 980 return [list $y $m $d $h $mn $s] 981 } 982 return "" 983} 984 985proc DateToSecs {y m d h mn s} { 986 # convert date to seconds ellapsed since beginning of $YEAR0, not 987 # necessarily a leap year 988 global YEAR0 989 990 return [DateToSecsFrom $y $m $d $h $mn $s $YEAR0] 991} 992 993proc DateToSecsFrom {y m d h mn s year0} { 994 # convert date to seconds ellapsed since beginning of $year0, not 995 # necessarily a leap year 996 global DAYSOF 997 998 # make sure there are no leading 0's 999 foreach v "y m d h mn s" { 1000 scan [set $v] %0d $v 1001 } 1002 # days in current month and in whole years except February 29ths 1003 set days [expr 365*($y-$year0)+$d-1] 1004 set yy [expr $year0+$year0%4] 1005 while { $yy < $y } { 1006 if { $yy%100!=0 || $yy%400==0 } { incr days } 1007 incr yy 4 1008 } 1009 # add days in whole months of current year 1010 if { $m>2 && $y%4==0 && ($y%100!=0 || $y%400==0) } { incr days } 1011 incr m -1 1012 while { $m > 0 } { incr days $DAYSOF($m) ; incr m -1 } 1013 return [expr $days*86400+$h*3600+$mn*60+$s] 1014} 1015 1016proc DateFromSecs {secs} { 1017 # build date from seconds ellapsed since beginning of $YEAR0, not 1018 # necessarily a leap year, using default format 1019 global DateFormat 1020 1021 return [eval FormatDate $DateFormat [DateIntsFromSecs $secs]] 1022} 1023 1024proc DateIntsFromSecs {secs} { 1025 # compute date from seconds ellapsed since beginning of $YEAR0, 1026 # not necessarily a leap year 1027 # return list with year, month, day, hour, minutes, seconds as integers 1028 global DAYSOF YEAR0 1029 1030 set s [expr $secs%60] ; set x [expr ($secs-$s)/60] 1031 set mn [expr $x%60] ; set x [expr ($x-$mn)/60] 1032 set h [expr $x%24] ; set x [expr ($x-$h)/24] 1033 set y [expr int($x/365)] ; set yd [expr $y*365] 1034 set yy [expr $YEAR0+$YEAR0%4] ; incr y $YEAR0 1035 while { $yy < $y } { 1036 if { $yy%100!=0 || $yy%400==0 } { incr yd } 1037 incr yy 4 1038 } 1039 if { $yd > $x } { 1040 incr y -1 1041 if { $y%4==0 && ($y%100!=0 || $y%400==0) } { incr yd -1 } 1042 incr yd -365 1043 } 1044 incr x [expr 1-$yd] ; set m 1 1045 while { $x > $DAYSOF($m) } { 1046 if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } { 1047 if { $x == 29 } { break } 1048 incr x -1 1049 } 1050 incr x -$DAYSOF($m) ; incr m 1051 } 1052 return [list $y $m $x $h $mn $s] 1053} 1054 1055proc DateFromSecsFmt {secs fmt} { 1056 # build date from seconds ellapsed since beginning of $YEAR0, not 1057 # necessarily a leap year, using given format 1058 global DAYSOF YEAR0 1059 1060 set s [expr $secs%60] ; set x [expr ($secs-$s)/60] 1061 set mn [expr $x%60] ; set x [expr ($x-$mn)/60] 1062 set h [expr $x%24] ; set x [expr ($x-$h)/24] 1063 set y [expr int($x/365)] ; set yd [expr $y*365] 1064 set yy [expr $YEAR0+$YEAR0%4] ; incr y $YEAR0 1065 while { $yy < $y } { 1066 if { $yy%100!=0 || $yy%400==0 } { incr yd } 1067 incr yy 4 1068 } 1069 if { $yd > $x } { 1070 incr y -1 1071 if { $y%4==0 && ($y%100!=0 || $y%400==0) } { incr yd -1 } 1072 incr yd -365 1073 } 1074 incr x [expr 1-$yd] ; set m 1 1075 while { $x > $DAYSOF($m) } { 1076 if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } { 1077 if { $x == 29 } { break } 1078 incr x -1 1079 } 1080 incr x -$DAYSOF($m) ; incr m 1081 } 1082 return [FormatDate $fmt $y $m $x $h $mn $s] 1083} 1084 1085proc NextDay {y m d} { 1086 # return list with year, month and day for the day following the 1087 # given date 1088 global DAYSOF 1089 1090 if { [incr d] > $DAYSOF($m) } { 1091 if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } { 1092 if { $d != 29 } { 1093 set d 1 ; set m 3 1094 } 1095 } else { 1096 set d 1 1097 if { [incr m] == 13 } { 1098 incr y ; set m 1 1099 } 1100 } 1101 } 1102 return [list $y $m $d] 1103} 1104 1105proc PreviousDay {y m d} { 1106 # return list with year, month and day for the day before the 1107 # given date 1108 global DAYSOF 1109 1110 if { [incr d -1] == 0 } { 1111 if { [incr m -1] == 0 } { 1112 set m 12 ; incr y -1 1113 } 1114 if { $m==2 && $y%4==0 && ($y%100!=0 || $y%400==0) } { 1115 set d 29 1116 } else { 1117 set d $DAYSOF($m) 1118 } 1119 } 1120 return [list $y $m $d] 1121} 1122 1123proc Month {name} { 1124 # from month name to 0 (error) or 1..12 1125 global ALLMONTH 1126 1127 for {set m 1} { $m<13 } { incr m } { 1128 if { [lsearch -exact $ALLMONTH($m) $name] != -1 } { return $m } 1129 } 1130 return 0 1131} 1132 1133proc LocalTimeAndUTC {y m d h mn s origin} { 1134 # convert from/to UTC to/from local time 1135 # $origin =="UTC" converts from UTC 1136 # =="local" converts to UTC 1137 # return list with same elements as in the args 1138 global TimeOffset DAYSOF 1139 1140 if { $origin != "UTC" } { 1141 incr h $TimeOffset 1142 } else { set h [expr $h-$TimeOffset] } 1143 if { $h < 0 } { 1144 incr h 24 1145 if { [incr d -1] == 0 } { 1146 if { [incr m -1] == 0 } { 1147 set m 12 ; set d 31 ; incr y -1 1148 } else { 1149 if { $m == 2 && $y%4 == 0 && \ 1150 ($y%100 != 0 || $y%400 == 0) } { 1151 set d 29 1152 } else { 1153 set d $DAYSOF($m) 1154 } 1155 } 1156 } 1157 } elseif { $h > 23 } { 1158 incr h -24 1159 if { [incr d] > $DAYSOF($m) } { 1160 if { $m == 2 } { 1161 if { $y%4 != 0 || ($y%100 == 0 && $y%400 != 0) } { 1162 set m 3 ; set d 1 1163 } 1164 } elseif { $m == 12 } { 1165 set m 1 ; set d 1 ; incr y 1166 } else { 1167 incr m ; set d 1 1168 } 1169 } 1170 } 1171 return [list $y $m $d $h $mn $s] 1172} 1173 1174proc CompatibleDates {y0 m0 d0 secs1} { 1175 # check whether two dates are compatible 1176 # $y0, $m0, $d0 give one date (year, month name and day) and may 1177 # be * (to match anything) 1178 # $secs1 is the other date in seconds since $YEAR0 1179 # return 1 if compatible 1180 1181 if { $m0 != "*" && [set m0 [Month $m0]] == 0 } { 1182 BUG calling CompatibleDates with bad month name 1183 } 1184 foreach "y1 m1 d1 h1 mn1 s1" [DateIntsFromSecs $secs1] {} 1185 if { [string match $d0-$m0-$y0 $d1-$m1-$y1] } { return 1 } 1186 return 0 1187} 1188 1189### varia 1190 1191proc VectorBearing {vx vy} { 1192 # compute bearing (angle from y-axis clockwise, in 0..359) of vector 1193 # with given components 1194 # return "_" if components are too small 1195 1196 if { abs($vy)+abs($vx) > 1e-20 } { 1197 set b [expr round(atan2($vx,$vy)*57.29577951308232087684)] 1198 if { $b < 0 } { incr b 360 } 1199 } else { set b "_" } 1200 return $b 1201} 1202 1203 1204