''Geodesy'' is the "science concerned with the study of the shape and size of the earth in the geometrical sense and with the study of certain physical phenomena, such as gravity, in seeking explanation of fine irregularities in the earth's shape." (D.H. Maling: Coordinate Systems and Map Projections, 2nd edition, Pergamon press 1992). So what has Tcl to do with that? All geographic applications like [GIS] need geodesy for their calculations. One of the important calculations is the measurement of the distance betwen two point on the earth giving the two geographical positions. This distance can be calculated using T. Vincenty's algorithm as published in ''Survey Review'' [http://www.surveyreview.org]: Vincenty T., (1975) "Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations", Survey Review 23, Number 176, p 88-93 The only thing we have to know beyond the two positions are the parameters of the ellipsoid in which the positions are expressed. This is important since the result will vary according to the used ellipsoid. So we define: # all known Ellipsoids: # each entry is a list with: ("major semi-axis","minor semi-axis","name") # axes are in metres set App(ellipsoid,wgs84) [list 6378137 6356750.321 "WGS84"] set App(ellipsoid,grs80) [list 6378137 6356752.314 "GRS80"] set App(ellipsoid,clarke) [list 6378206.4 6356583.8 "Clarke (1866) / NAD27"] set App(ellipsoid,int) [list 6378388 6356911.946 "International (1909/1924)"] set App(ellipsoid,krassovsky) [list 6378245 6356863.019 "Krassovsky (1940)"] set App(ellipsoid,bessel) [list 6377397.155 6356078.963 "Bessel (1841)"] set App(ellipsoid,wgs72) [list 6378135 6356750.52 "WGS72"] set App(ellipsoid,wgs66) [list 6378145 6356759.769 "WGS66"] set App(ellipsoid,airy) [list 6377563.396 6356256.909 "Airy 1830"] Often you will want to use the WGS84 or GRS80 ellipsoid since GPS data are expressed with reference to these. Now we need to define some constant values: set App(PI) [expr {atan(1)*4}] ;# this is pi set App(nm) 1852.0 ;# this is the number of meters that go into a nautical mile And then a conversion factor between degress (ranging from 0 to 360) and radians (ranging from 0 to 2*pi): set App(d2r) [expr {2*$App(PI)/360.0}] Now we are ready to go along with the algorithmic solution. It uses a loop to calculate the correct value iteratively until the result is good enough: # method -> of of the elements of the above ellipsoid array, e.g. 'ellipsoid,bessel' # x1 ... y2 -> position of the points 1 and 2 in decimal degrees # returns the distance between the points expressed in metres proc DistEllipsoid {method x1 y1 x2 y2} { global App # set our precision set eps 0.00000000005 # convert to radians: set lon1 [expr {$x1 * $App(d2r)}] set lon2 [expr {$x2 * $App(d2r)}] set lat1 [expr {$y1 * $App(d2r)}] set lat2 [expr {$y2 * $App(d2r)}] # get ellipsoid data: foreach {a0 b0 name} $App($method) break # flatness: set flat [expr {($a0 - $b0) / $a0}] set r [expr {1 - $flat}] # set tanu1 [expr {$r * tan($lat1)}] set tanu2 [expr {$r * tan($lat2)}] # set x [expr {atan($tanu1)}] set cosu1 [expr {cos($x)}] set sinu1 [expr {sin($x)}] # set x [expr {atan($tanu2)}] set cosu2 [expr {cos($x)}] set sinu2 [expr {sin($x)}] # set omega [expr {$lon2 - $lon1}] set lambda $omega # loop untile the result is good enough: while {1} { set testlambda $lambda set s2s [expr {pow($cosu2 * sin($lambda),2) + \ pow($cosu1 * $sinu2 - $sinu1 * $cosu2 * cos($lambda),2)}] set ss [expr {sqrt($s2s)}] set cs [expr {$sinu1 * $sinu2 + $cosu1 * $cosu2 * cos($lambda)}] set tansigma [expr {$ss / $cs}] set sinalpha [expr {$cosu1 * $cosu2 * sin($lambda) / $ss}] set x [expr {asin($sinalpha)}] set cosalpha [expr {cos($x)}] set c2sm [expr {$cs - (2 * $sinu1 * $sinu2 / pow($cosalpha,2))}] set c [expr {$flat / 16.0 * pow($cosalpha,2) * \ (4 + $flat * (4 - 3 * pow($cosalpha,2)))}] set lambda [expr {$omega + (1 - $c) * \ $flat * $sinalpha * (asin($ss) + $c * $ss * ($c2sm + \ $c * $cs * (-1 + 2 * pow($c2sm,2))))}] # result good enough? if {abs($testlambda - $lambda) <= 0.00000000000005} {break} } set u2 [expr {pow($cosalpha,2) * (pow($a0,2) - pow($b0,2)) / pow($b0,2)}] set a [expr {1 + ($u2 / 16384.0) * (4096 + $u2 * \ (-768 + $u2 * (320 - 175 * $u2)))}] set b [expr {($u2 / 1024.0) * (256 + $u2 * (-128 + $u2 * (74 - 47 * $u2)))}] set dsigma [expr {$b * $ss * ($c2sm + ($b / 4.0) * \ ($cs * (-1 + 2 * pow($c2sm,2)) - ($b / 6.0) * $c2sm * \ (-3 + 4 * pow($ss,2)) * (-3 + 4 * pow($c2sm,2))))}] return [expr {$b0 * $a * (asin($ss) - $dsigma)}] } Thats it! ---- Related pages: * [Canvas geometrical calculations] * [GIS] ---- [[ [CL] needs to make time to demonstrate how much simpler the calculation is for a sphere, and how little inaccuracy that introduces.]] [[Nice work above, though, Torsten.]] [TR] - you are right, calculation for a sphere is ''much'' simpler: set d [expr {acos(sin($y1)*sin($y2)+cos($y1)*cos($y2)*cos($x1-$x2))}] set distance_in_metres [expr {180*60.0/3.1415926535*$d*1852}] See [Canvas geometrical calculations] for details. Under the assumption that 1 nautical mile is the same as 1 Minute in latitude, I can give an example: Two points: x1=10 , y1=54 , x1=10.5 , y2=54 * spherical: 63735 metres * WGS84 (GPS): 63993 metres The difference is 258 metres. It depends on the application if this difference is significant enough.