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: projs_main.tcl
22#  Last change:  6 October 2013
23#
24
25# Includes contributions by
26#  - Sandor Laky (laky.sandor_at_freemail.hu) marked "SL contribution"
27#  - Miguel Filgueiras (to code by others) marked "MF contribution/change"
28
29### most of this code used to be at the end of file projections.tcl
30
31### general initialization procedure for all main projections
32# see also the comments before the procedures for UTM in projections.tcl
33
34proc ProjInit {proj data datum ps} {
35    # changes here must be reflected on proc BadProjSetup below
36    # initialize projection parameters
37    #  $ps is a list of latd,longd,datum to be projected
38    # should call proc ProjParams if its parameters may be changed by
39    #  the user and if $ASKPROJPARAMS is true
40    # should call proc to compute auxiliary parameters if any
41    # must return projection of first position in $ps
42    global $data TMPARAM LNTFPARAMS ASKPROJPARAMS
43
44    foreach "latd longd posdatum" [lindex $ps 0] { break }
45    if { $posdatum != $datum } {
46	foreach "latd longd" \
47		[ToDatum $latd $longd $posdatum $datum] { break }
48    }
49    set ${data}(datum) $datum
50    set askauxcomp 1
51    switch $proj {
52	UTM {
53	    # only the first position is used
54	    # no parameters can be changed by the user
55	    set askauxcomp 0
56	    set p [DegreesToUTM $latd $longd $datum]
57	    scan [lindex $p 0] %0d ze
58	    set zn [lindex $p 1]
59	    set ${data}(UTMzone) [list $ze $zn]
60	    set ${data}(m_0) [expr -183+6.0*$ze]
61	    set res [list [lindex $p 2] [lindex $p 3]]
62	}
63	TM {
64	    # take averages of latitudes and of longitudes as central
65	    #  coordinates, and use k0=0.9996
66	    foreach "${data}(lat0) ${data}(long0)" \
67		    [AverageLatLong $latd $longd $datum $ps] {}
68	    set ${data}(k0) 0.9996
69	    if { $ASKPROJPARAMS } {
70		ProjParams change TM $data
71	    }
72	    set res [ConvToTM $latd $longd [set ${data}(lat0)] \
73		          [set ${data}(long0)] [set ${data}(k0)] $datum]
74	    set askauxcomp 0
75	}
76	BMN -  CTR -  GKK -  KKJP -  TWG {
77	    # special cases of TM
78	    # take averages of longitudes to find long0, and use lat0=0
79	    set long0 [AverageLong $longd $datum $ps]
80	    foreach "pr lmin lmax delta k0" $TMPARAM($proj) {}
81	    if { $long0 < $lmin } {
82		set long0 $lmin
83	    } elseif { $long0 > $lmax } {
84		set long0 $lmax
85	    } else {
86		set long0 [expr int(1.0*($long0-$lmin+0.5*$delta)/ \
87			            $delta)*$delta+$lmin]
88	    }
89	    set ${data}(${pr}long0) $long0
90	    set ${data}(m_0) $k0
91	    if { $ASKPROJPARAMS } {
92		ProjParams change $proj $data
93	    }
94	    set res [ConvToTM $latd $longd 0 [set ${data}(${pr}long0)] \
95		              $k0 $datum]
96	    set askauxcomp 0
97	}
98	LCC1 -  Stereogr -  SOM {
99	    # Lambert Conformal Conic single parallel,
100	    #  Stereographic and Swiss Oblique Mercator projections
101	    # take averages of latitudes and of longitudes as central
102	    #  coordinates, and use k0=1
103	    # in Stereographic central latitude determines the case: polar,
104	    #  oblique, equatorial
105	    foreach "${data}(lat0) ${data}(long0)" \
106		    [AverageLatLong $latd $longd $datum $ps] {}
107	    set ${data}(k0) 1.0
108	}
109	LCC2 {
110	    # Lambert Conformal Conic projection, two parallels
111	    # take extremes of latitudes as latitudes of standard parallels,
112	    #  averages of latitudes and longitudes as false origin coordinates
113	    foreach "${data}(latF) ${data}(longF) \
114		     ${data}(lat1) ${data}(lat2)" \
115		    [AverageExtrLatLong $latd $longd $datum $ps] {}
116	}
117	LambNTF {
118	    # as LCC2 but with zones with specific parameters
119	    # only the first position is used
120	    # no parameters can be changed by the user
121	    set askauxcomp 0
122	    set p [DegreesToLambNTF $latd $longd $datum]
123	    if { [set z [lindex $p 0]] == "--" } {
124		set z II
125	    }
126	    global LNTFz$z
127	    array set $data [array get LNTFz$z]
128	    set res [ProjLambNTFPoint $data $latd $longd $datum]
129	}
130	AEA {
131	    # Albers Equal Area Conic projection
132	    # take extremes of latitudes as latitudes of standard parallels,
133	    #  average of longitudes as central longitude, 0 for central
134	    #  latitude
135	    foreach "x ${data}(long0) lamax lamin" \
136		    [AverageExtrLatLong $latd $longd $datum $ps] {}
137	    if { abs($lamax) < abs($lamin) } {
138		set x $lamax ; set lamax $lamin ; set lamin $x
139	    }
140	    set ${data}(lat1) $lamax ; set ${data}(lat2) $lamin
141	    set ${data}(lat0) 0
142	}
143	LEAC {
144	    # Lambert Equal Area Conic projection
145	    # take extreme of latitudes as latitude of standard parallel,
146	    #  average of longitudes as central longitude, North polar
147	    #  aspect if extreme latitude is nonnegative, 0 for central
148	    #  latitude
149	    foreach "x ${data}(long0) lamax lamin" \
150		    [AverageExtrLatLong $latd $longd $datum $ps] {}
151	    if { abs($lamax) < abs($lamin) } {
152		set lamax $lamin
153	    }
154	    set ${data}(lat1) $lamax
155	    if { $lamax < 0 } {
156		set ${data}(polasp) south
157	    } else { set ${data}(polasp) north }
158	    set ${data}(lat0) 0
159	}
160	Merc1 {
161	    # Mercator projection, single parallel, ellipsoidal case
162	    # take average of longitudes as central longitude, and use k0=1
163	    set ${data}(long0) [AverageLong $longd $datum $ps]
164	    set ${data}(k0) 1.0
165	}
166	Merc2 {
167	    # Mercator projection, two parallels
168	    # take average of longitudes as central longitude and extreme
169	    #  latitude (in absolute value) as first standard parallel
170	    foreach "x ${data}(long0) lamin lamax" \
171		    [AverageExtrLatLong $latd $longd $datum $ps] {}
172	    if { abs($lamin) > $lamax } { set lamax [expr abs($lamin)] }
173	    set ${data}(lat1) $lamax
174	}
175	SphMerc -  CS {
176	    # Mercator projection, single parallel, spherical case
177            #  and Cassini-Soldner projection
178	    # take average of longitudes and longitudes as central coordinates
179	    foreach "${data}(lat0) ${data}(long0)" \
180		    [AverageLatLong $latd $longd $datum $ps] {}
181	}
182	APOLY -  EqCyl {
183	    # American Polyconic and Equidistant Cylindrical projections
184	    # take averages of latitudes as standard latitude
185	    set ${data}(lat0) [AverageLat $latd $datum $ps]
186	}
187	Schreiber {
188	    # Schreiber projection
189	    # all parameters are fixed, except the datum that must use the
190	    #  "Bessel 1841" ellipsoid
191	    if { [EllipsdOf $datum] != "Bessel 1841" } {
192		GMMessage "Datum: Rijks Driehoeksmeting"
193		set datum "Rijks Driehoeksmeting"
194	    }
195	    ProjSchreiberComputeAux $data $datum
196	    set res [ProjSchreiberPoint $data $latd $longd $datum]
197	    set askauxcomp 0
198	}
199	EOV {
200	    # Hungarian grid projection: no parameters needed
201	    set res [ProjEOVPoint {} $latd $longd $datum]
202	}
203	default {
204	    BUG "no initialization proc for projection"
205	    set askauxcomp 0
206	}
207    }
208    if { $askauxcomp } {
209	if { $ASKPROJPARAMS } {
210	    ProjParams change $proj $data
211	}
212	Proj${proj}ComputeAux $data $datum
213	set res [Proj${proj}Point $data $latd $longd $datum]
214    }
215    return $res
216}
217
218proc BadProjSetup {data proj latd longd datum} {
219    # this must be kept in line with proc ProjInit (above)
220    # initialize projection for use in command mode
221    #  $data is array for the parameters with external parameters and
222    #    datum defined, except in the case of "UTM" and "EOV" that only
223    #    have datum
224    # return 0 on success, 1 on error
225    global ASKPROJPARAMS $data
226
227    set ASKPROJPARAMS 0
228    set aux 0
229    switch $proj {
230	UTM {
231	    if { [catch {set ${data}(UTMzone)}] } {
232		set p [DegreesToUTM $latd $longd $datum]
233		scan [lindex $p 0] %0d ze
234		set zn [lindex $p 1]
235		set ${data}(UTMzone) [list $ze $zn]
236		set ${data}(m_0) [expr -183+6.0*$ze]
237	    }
238	}
239	LambNTF {
240	    # as LCC2 but with zones with specific parameters
241	    # no parameters can be changed by the user
242	    set p [DegreesToLambNTF $latd $longd $datum]
243	    if { [set z [lindex $p 0]] == "--" } {
244		set z II
245	    }
246	    global LNTFz$z
247	    catch {unset $data}
248	    array set $data [array get LNTFz$z]
249	}
250	Schreiber {
251	    # Schreiber projection
252	    # all parameters are fixed, except the datum that must use the
253	    #  "Bessel 1841" ellipsoid
254	    if { [EllipsdOf [set ${data}(datum)]] != "Bessel 1841" } {
255		set ${data}(datum) "Rijks Driehoeksmeting"
256	    }
257	}
258	TM -  BMN -  CTR -  EOV -  GKK -  KKJP -  TWG {
259	}
260	LCC1 -  Stereogr -  SOM -  LCC2 -  AEA -
261	LEAC -  Merc1 -  Merc2 -  SphMerc -  CS -  APOLY -  EqCyl {
262	    set aux 1
263	}
264	default {
265	    return 1
266	}
267    }
268    if { $aux } {
269	Proj${proj}ComputeAux $data [set ${data}(datum)]
270    }
271    return 0
272}
273
274proc AverageExtrLatLong {latd longd datum ps} {
275    # compute the averages of latitudes and longitudes and the minimum and
276    #  maximum latitude for $latd,$longd and all positions in $ps except
277    #  the first
278    # return list with the two averages followed by the maximum and the minimum
279
280    set n 1
281    set sla [set lamx [set lamn $latd]] ; set slo $longd
282    foreach p [lreplace $ps 0 0] {
283	foreach "la lo posdatum" $p { break }
284	if { $posdatum != $datum } {
285	    foreach "la lo" \
286		    [ToDatum $la $lo $posdatum $datum] { break }
287	}
288	set sla [expr $sla+$la] ; set slo [expr $slo+$lo]
289	if { $la > $lamx } { set lamx $la }
290	if { $la < $lamn } { set lamn $la }
291	incr n
292    }
293    return [list [expr 1.0*$sla/$n] [expr 1.0*$slo/$n] $lamx $lamn]
294}
295
296proc AverageLatLong {latd longd datum ps} {
297    # compute the averages of latitudes and longitudes for $latd,$longd and
298    #  all positions in $ps except the first
299    # return list with the two averages
300
301    set n 1
302    foreach p [lreplace $ps 0 0] {
303	foreach "la lo posdatum" $p { break }
304	if { $posdatum != $datum } {
305	    foreach "la lo" \
306		    [ToDatum $la $lo $posdatum $datum] { break }
307	}
308	set latd [expr $latd+$la] ; set longd [expr $longd+$lo]
309	incr n
310    }
311    return [list [expr 1.0*$latd/$n] [expr 1.0*$longd/$n]]
312}
313
314proc AverageLat {latd datum ps} {
315    # compute the average of latitudes for $latd and
316    #  all positions in $ps except the first
317
318    set n 1
319    foreach p [lreplace $ps 0 0] {
320	set la [lindex $p 0] ; set pdatum [lindex $p 2]
321	if { $pdatum != $datum } {
322	    set la [lindex [ToDatum $la [lindex $p 1] $pdatum $datum] 0]
323	}
324	set latd [expr $la+$latd]
325	incr n
326    }
327    return [expr 1.0*$latd/$n]
328}
329
330proc AverageLong {longd datum ps} {
331    # compute the average of longitudes for $longd and
332    #  all positions in $ps except the first
333
334    set n 1
335    foreach p [lreplace $ps 0 0] {
336	set lo [lindex $p 1] ; set pdatum [lindex $p 2]
337	if { $pdatum != $datum } {
338	    set lo [lindex [ToDatum [lindex $p 0] $lo $pdatum $datum] 1]
339	}
340	set longd [expr $longd+$lo]
341	incr n
342    }
343    return [expr 1.0*$longd/$n]
344}
345
346## converting longitudes to -180, 180 (inclusive) or -360, 360 (inclusive)
347
348proc NormalLong {long} {
349    # convert longitude in decimal degrees to the -180,180 (inclusive) range
350
351    if { abs([set x [expr abs($long)*0.00555555555555555556]])-1 > 1e-14 } {
352	set x [expr 0.5*($x+1)]
353	set x [expr (($x-floor($x))-0.5)*360]
354    } else { set x [expr abs($long)] }
355    if { $long < 0 && abs($x-180) > 1e-14 } { return [expr -$x] }
356    return $x
357}
358
359proc NormalLongCentred {long longc} {
360    # convert longitude in decimal degrees to the closest value to $longc in
361    #  the -360,360 (inclusive) range
362
363    set long [NormalLong $long]
364    if { $long <= $longc-180 } {
365	set long [expr $long+360]
366    } elseif { $long > $longc+180 } {
367	set long [expr $long-360]
368    }
369    return $long
370}
371
372## projection procedures for UTM projection
373
374# there may be auxiliary parameters to be computed from the main ones
375
376proc ProjUTMComputeAux {data args} {
377    # compute auxiliary parameters from main parameters
378    #  $data is name of global array for the parameters
379    #  $args is not used but is needed because of uniform way of calling
380    #   these procedures
381    global $data
382
383    set z [set ${data}(UTMzone)]
384    # set the zone to a list of number and letter if it is a string
385    if { [regexp {^([0-9]+)([A-Z])$} $z x ze zn] } {
386	set ${data}(UTMzone) [list $ze $zn]
387    } else { set ze [lindex $z 0] }
388    set ${data}(m_0) [expr -183+6.0*$ze]
389    return
390}
391
392# projection proc is a function of latd,longd,datum to planar Cartesian
393#  coordinates x,y in the terrain
394
395proc ProjUTMPoint {data latd longd datum} {
396    # compute planar Cartesian coordinates for given position
397    global $data
398
399    set prdatum [set ${data}(datum)]
400    if { $datum != $prdatum } {
401	foreach "latd longd" \
402		[ToDatum $latd $longd $datum $prdatum] { break }
403    }
404    foreach "pze pzn px py" [DegreesToUTM $latd $longd $prdatum] {}
405    foreach "mze mzn" [set ${data}(UTMzone)] {}
406    if { $mze!=$pze || $mzn != $pzn } {
407	foreach "px py" [CompUTMOnZone $latd $longd [set ${data}(m_0)]] {
408	    break
409	}
410    }
411    return [list $px $py]
412}
413
414# inverse projection proc is a function of planar Cartesian coordinates in the
415#  terrain to list with latitude and longitude in the datum of the projection
416
417proc ProjUTMInvert {data x y} {
418    # return list with latitude and longitude for point at $x,$y in the terrain
419    global $data
420
421    foreach "ze zn" [set ${data}(UTMzone)] {}
422    set datum [set ${data}(datum)]
423    set p [UTMToDegrees $ze $zn $x $y $datum]
424    return $p
425}
426
427## Transverse Mercator projection
428
429proc ProjTMPoint {data latd longd datum} {
430    # compute planar Cartesian coordinates for given position
431    global $data
432
433    set prdatum [set ${data}(datum)]
434    if { $datum != $prdatum } {
435	foreach "latd longd" \
436		[ToDatum $latd $longd $datum $prdatum] { break }
437    }
438    return [ConvToTM $latd $longd [set ${data}(lat0)] \
439	    [set ${data}(long0)] [set ${data}(k0)] $prdatum]
440}
441
442proc ProjTMInvert {data x y} {
443    # return list with latitude and longitude for point at $x,$y in the terrain
444    global $data
445
446    set datum [set ${data}(datum)]
447    set p [ConvFromTM $x $y [set ${data}(lat0)] [set ${data}(long0)] \
448	    [set ${data}(k0)] $datum]
449    return $p
450}
451
452## German Grid (Gauss-Krueger-Koordinatensystem) projection
453# Information provided by Andreas Lange (Andreas.C.Lange_at_GMX.de)
454# This is a Transverse Mercator projection with zone in [0-6], lat0=0,
455#  lon0=zone*3 (0, 3, 6, 9, 12, 15E), and scale factor at central meridian
456#  of k0=1.0
457
458# Basic Finnish Grid (KKJP) projection
459# Similar with zone in [1-4], lon0=zone*3+18 (21, 24, 27, 30E)
460
461# Taiwan Grid projection
462# Information provided by Dan Jacobson (jidanni_at_yahoo.com.tw)
463# Similar with zone in [1-6], lon0=zone*2+113 (115, 117, ..., 125E),
464#  and k0=0.9999
465
466# Carta Tecnica Regionale (CTR), Italian projection
467# "Le carte topografiche CTR ed il loro uso GPS"
468# (http://www.gpscomefare.com/guide/tutorialgps/mapdatum.htm)
469# May 2003 (information kindly sent by Alessandro Palmas)
470# Similar with zone in [1-2], lon0=zone*6+3 (9, 15), k0=0.9996
471
472# Austrian BMN grid projection
473# information kindly sent by Alessandro Palmas, July 2003
474# Similar with zones M28, M31, M34, lon0=zone_index*3+10.333333, k0=1
475
476proc ProjGKKComputeAux {data args} {
477    # compute auxiliary parameters from main parameters
478    #  $data is name of global array for the parameters
479    #  $args is not used but is needed because of uniform way of calling
480    #   these procedures
481    global $data
482
483    set ${data}(m_0) 1.0
484    return
485}
486
487proc ProjKKJPComputeAux {data args} {
488    # compute auxiliary parameters from main parameters
489    #  $data is name of global array for the parameters
490    #  $args is not used but is needed because of uniform way of calling
491    #   these procedures
492    global $data
493
494    set ${data}(m_0) 1.0
495    return
496}
497
498proc ProjTWGComputeAux {data args} {
499    # compute auxiliary parameters from main parameters
500    #  $data is name of global array for the parameters
501    #  $args is not used but is needed because of uniform way of calling
502    #   these procedures
503    global $data
504
505    set ${data}(m_0) 0.9999
506    return
507}
508
509proc ProjCTRComputeAux {data args} {
510    # compute auxiliary parameters from main parameters
511    #  $data is name of global array for the parameters
512    #  $args is not used but is needed because of uniform way of calling
513    #   these procedures
514    global $data
515
516    set ${data}(m_0) 0.9996
517    return
518}
519
520proc ProjBMNComputeAux {data args} {
521    # compute auxiliary parameters from main parameters
522    #  $data is name of global array for the parameters
523    #  $args is not used but is needed because of uniform way of calling
524    #   these procedures
525    global $data
526
527    set ${data}(m_0) 1.0
528    return
529}
530
531proc ProjGKKPoint {data latd longd datum} {
532    # compute planar Cartesian coordinates for given position
533
534    return [Proj_TMZ_Point gkk $data $latd $longd $datum]
535}
536
537proc ProjKKJPPoint {data latd longd datum} {
538    # compute planar Cartesian coordinates for given position
539
540    return [Proj_TMZ_Point kkjp $data $latd $longd $datum]
541}
542
543proc ProjTWGPoint {data latd longd datum} {
544    # compute planar Cartesian coordinates for given position
545
546    return [Proj_TMZ_Point twg $data $latd $longd $datum]
547}
548
549proc ProjCTRPoint {data latd longd datum} {
550    # compute planar Cartesian coordinates for given position
551
552    return [Proj_TMZ_Point ctr $data $latd $longd $datum]
553}
554
555proc ProjBMNPoint {data latd longd datum} {
556    # compute planar Cartesian coordinates for given position
557
558    return [Proj_TMZ_Point bmn $data $latd $longd $datum]
559}
560
561proc Proj_TMZ_Point {pr data latd longd datum} {
562    # compute planar Cartesian coordinates for given position
563    global $data
564
565    set prdatum [set ${data}(datum)]
566    if { $datum != $prdatum } {
567	foreach "latd longd" \
568		[ToDatum $latd $longd $datum $prdatum] { break }
569    }
570    return [ConvToTM $latd $longd 0 [set ${data}(${pr}long0)] \
571	             [set ${data}(m_0)] $prdatum]
572}
573
574proc ProjGKKInvert {data x y} {
575    # return list with latitude and longitude for point at $x,$y in the terrain
576
577    return [Proj_TMZ_Invert gkk $data $x $y]
578}
579
580proc ProjKKJPInvert {data x y} {
581    # return list with latitude and longitude for point at $x,$y in the terrain
582
583    return [Proj_TMZ_Invert kkjp $data $x $y]
584}
585
586proc ProjTWGInvert {data x y} {
587    # return list with latitude and longitude for point at $x,$y in the terrain
588
589    return [Proj_TMZ_Invert twg $data $x $y]
590}
591
592proc ProjCTRInvert {data x y} {
593    # return list with latitude and longitude for point at $x,$y in the terrain
594
595    return [Proj_TMZ_Invert ctr $data $x $y]
596}
597
598proc ProjBMNInvert {data x y} {
599    # return list with latitude and longitude for point at $x,$y in the terrain
600
601    return [Proj_TMZ_Invert bmn $data $x $y]
602}
603
604proc Proj_TMZ_Invert {pr data x y} {
605    # return list with latitude and longitude for point at $x,$y in the terrain
606    global $data
607
608    set datum [set ${data}(datum)]
609    set p [ConvFromTM $x $y 0 [set ${data}(${pr}long0)] \
610	              [set ${data}(m_0)] $datum]
611    return $p
612}
613
614## Lambert Conic Conformal projections
615
616# single standard parallel
617
618proc ProjLCC1ComputeAux {data datum} {
619    # compute auxiliary parameters from main parameters
620    #  $data is name of global array for the parameters
621    global $data
622
623    set d [EllipsdData $datum]
624    set a [lindex $d 0] ; set f [lindex $d 1]
625    set es [expr $f*(2-$f)]
626    set e [set ${data}(m_e) [expr sqrt($es)]]
627    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
628    if { [set n [set ${data}(lcc_n) [expr sin($phi0)]]] < 0 } {
629	set ${data}(lcc_sn) -1
630    } else { set ${data}(lcc_sn) 1 }
631    set m0 [expr cos($phi0)/sqrt(1-$es*$n*$n)]
632    set t0n [expr pow([ExpRedLat $phi0 $n $e],$n)]
633    set F [expr $m0/($n*$t0n)]
634    set aFk0 [set ${data}(m_a) [expr $a*$F*[set ${data}(k0)]]]
635    set ${data}(lcc_rho0F) [expr $aFk0*$t0n]
636    return
637}
638
639proc ProjLCC1Point {data latd longd datum} {
640    # compute planar Cartesian coordinates for given position
641    global $data
642
643    set prdatum [set ${data}(datum)]
644    if { $datum != $prdatum } {
645	foreach "latd longd" \
646		[ToDatum $latd $longd $datum $prdatum] { break }
647    }
648    set long0 [set ${data}(long0)]
649    set longd [NormalLongCentred $longd $long0]
650    set phi [expr $latd*0.01745329251994329576]
651    set theta [expr [set ${data}(lcc_n)]* \
652	            ($longd-$long0)*0.01745329251994329576]
653    set t [ExpRedLat $phi [expr sin($phi)] [set ${data}(m_e)]]
654    set rho [expr [set ${data}(m_a)]* \
655	          pow($t,[set ${data}(lcc_n)])]
656    set x [expr $rho*sin($theta)]
657    set y [expr [set ${data}(lcc_rho0F)]-$rho*cos($theta)]
658    return [list $x $y]
659}
660
661proc ProjLCC1Invert {data x y} {
662    # return list with latitude and longitude for point at $x,$y in the terrain
663    global $data
664
665    set z [expr [set ${data}(lcc_rho0F)]-$y]
666    set thetap [expr atan(1.0*$x/$z)]
667    set rhop [expr [set ${data}(lcc_sn)]*[Hypot $x $z]]
668    set tp [expr pow($rhop/([set ${data}(m_a)]), \
669	             1.0/[set ${data}(lcc_n)])]
670    set longd [expr $thetap/[set ${data}(lcc_n)]*57.29577951308232087684+ \
671	            [set ${data}(long0)]]
672    set longd [NormalLong $longd]
673    set latd [expr [LatFromRedLat $tp [set ${data}(m_e)]]* \
674	           57.29577951308232087684]
675    return [list $latd $longd]
676}
677
678# two standard parallels
679
680proc ProjLCC2ComputeAux {data datum} {
681    # compute auxiliary parameters from main parameters
682    #  $data is name of global array for the parameters
683    global $data MESS
684
685    set d [EllipsdData $datum]
686    set a [lindex $d 0] ; set f [lindex $d 1]
687    set es [expr $f*(2-$f)]
688    set e [set ${data}(m_e) [expr sqrt($es)]]
689    set phi1 [expr [set ${data}(lat1)]*0.01745329251994329576]
690    set phi2 [expr [set ${data}(lat2)]*0.01745329251994329576]
691    set phiF [expr [set ${data}(latF)]*0.01745329251994329576]
692    set s1 [expr sin($phi1)] ; set s2 [expr sin($phi2)]
693    set t1 [ExpRedLat $phi1 $s1 $e]
694    set t2 [ExpRedLat $phi2 $s2 $e]
695    if { [catch {set m1 [expr cos($phi1)/sqrt(1-$es*$s1*$s1)]}] || \
696	    [catch {set m2 [expr cos($phi2)/sqrt(1-$es*$s2*$s2)]}] || \
697	    [catch {set n [expr (log($m1)-log($m2))/(log($t1)-log($t2))]}] } {
698	GMMessage $MESS(badProjargs)
699	ProjParams change LCC2 $data
700	ProjLCC2ComputeAux $data $datum
701	return
702    }
703    set ${data}(lcc_n) $n
704    if { $n < 0 } {
705	set ${data}(lcc_sn) -1
706    } else { set ${data}(lcc_sn) 1 }
707    set F [expr $m1/($n*pow($t1,$n))]
708    set aF [set ${data}(m_a) [expr $a*$F]]
709    set tF [ExpRedLat $phiF [expr sin($phiF)] $e]
710    set ${data}(lcc_rho0F) [expr $aF*pow($tF,$n)]
711    return
712}
713
714proc ProjLCC2Point {data latd longd datum} {
715    # compute planar Cartesian coordinates for given position
716    global $data
717
718    set prdatum [set ${data}(datum)]
719    if { $datum != $prdatum } {
720	foreach "latd longd" \
721		[ToDatum $latd $longd $datum $prdatum] { break }
722    }
723    set longF [set ${data}(longF)]
724    set longd [NormalLongCentred $longd $longF]
725    set phi [expr $latd*0.01745329251994329576]
726    set theta [expr [set ${data}(lcc_n)]* \
727	            ($longd-$longF)*0.01745329251994329576]
728    set t [ExpRedLat $phi [expr sin($phi)] [set ${data}(m_e)]]
729    set rho [expr [set ${data}(m_a)]* \
730	          pow($t,[set ${data}(lcc_n)])]
731    set x [expr $rho*sin($theta)]
732    set y [expr [set ${data}(lcc_rho0F)]-$rho*cos($theta)]
733    return [list $x $y]
734}
735
736proc ProjLCC2Invert {data x y} {
737    # return list with latitude and longitude for point at $x,$y in the terrain
738    global $data
739
740    set z [expr [set ${data}(lcc_rho0F)]-$y]
741    set thetap [expr atan(1.0*$x/$z)]
742    set rhop [expr [set ${data}(lcc_sn)]*[Hypot $x $z]]
743    set tp [expr pow($rhop/([set ${data}(m_a)]), \
744	             1.0/[set ${data}(lcc_n)])]
745    set longd [expr $thetap/[set ${data}(lcc_n)]*57.29577951308232087684+ \
746	            [set ${data}(longF)]]
747    set longd [NormalLong $longd]
748    set latd [expr [LatFromRedLat $tp \
749	             [set ${data}(m_e)]]*57.29577951308232087684]
750    return [list $latd $longd]
751}
752
753# French NTF projection
754
755proc ProjLambNTFComputeAux {data datum} {
756    global $data
757
758    set z [set $data(NTFzone)]
759    global LNTFz$z
760    array set $data [array get LNTFz$z]
761    return
762}
763
764proc ProjLambNTFPoint {data latd longd datum} {
765
766    return [ProjLCC2Point $data $latd $longd $datum]
767}
768
769proc ProjLambNTFInvert {data x y} {
770
771    return [ProjLCC2Invert $data $x $y]
772}
773
774# ancillary procs for LCC projections
775
776proc ExpRedLat {phi sinphi e} {
777    # compute exp of reduced latitude for $phi
778    #  $e is first eccentricity of ellipsoid
779
780    set se [expr $sinphi*$e]
781    return [expr tan((1.5707963267948966-$phi)/2.0) / \
782	    pow((1.0-$se)/(1.0+$se), $e/2.0)]
783}
784
785proc LatFromRedLat {ts e} {
786    # compute latitude from reduced latitude
787
788    set eccnth [expr $e/2.0]
789    set phi [expr 1.5707963267948966-2*atan($ts)]
790    set i 16 ; set d 1
791    while { abs($d) > 1e-10 && [incr i -1] } {
792	set con [expr $e*sin($phi)]
793	set d [expr 1.5707963267948966- \
794		    2*atan($ts*pow((1-$con)/(1+$con),$eccnth))-$phi]
795	set phi [expr $phi+$d]
796    }
797    return $phi
798}
799
800## Albers Equal Area Conic and Lambert Equal Area Conic projection
801# adapted from PROJ4.3.3
802
803proc ProjAEAComputeAux {data datum} {
804    # compute auxiliary parameters from main parameters
805    #  $data is name of global array for the parameters
806    global $data
807
808    set phi1 [expr [set ${data}(lat1)]*0.01745329251994329576]
809    set phi2 [expr [set ${data}(lat2)]*0.01745329251994329576]
810    ComputeAux_AEA_LEAC AEA $data $datum $phi1 $phi2
811    return
812}
813
814proc ProjAEAPoint {data latd longd datum} {
815    # compute planar Cartesian coordinates for given position
816    global $data MESS
817
818    set prdatum [set ${data}(datum)]
819    if { $datum != $prdatum } {
820	foreach "latd longd" \
821		[ToDatum $latd $longd $datum $prdatum] { break }
822    }
823    set long0 [set ${data}(long0)]
824    set longd [NormalLongCentred $longd $long0]
825    set a [set ${data}(m_a)] ; set n [set ${data}(m_1)]
826    set phi [expr $latd*0.01745329251994329576]
827    set lam [expr $n*($longd-$long0)*0.01745329251994329576]
828    set rho [expr [set ${data}(m_3)]-$n* \
829	            [SmallQ [expr sin($phi)] [set ${data}(m_e)] \
830		            [set ${data}(m_0)]]]
831    if { $rho < 0 } {
832	GMMessage $MESS(outrngproj)
833	return [list 0 0]
834    }
835    set rho [expr [set ${data}(m_4)]*sqrt($rho)]
836    set x [expr $a*$rho*sin($lam)]
837    set y [expr $a*([set ${data}(m_5)]-$rho*cos($lam))]
838    return [list $x $y]
839}
840
841proc ProjAEAInvert {data x y} {
842    # return list with latitude and longitude for point at $x,$y in the terrain
843    global $data MESS TXT
844
845    set a [set ${data}(m_a)] ; set n [set ${data}(m_1)]
846    set datum [set ${data}(datum)]
847    set x [expr $x/$a]
848    set y [expr [set ${data}(m_5)]-$y/$a]
849    if { [set rho [Hypot $x $y]] == 0 } {
850	set lam 0
851	if { $n > 0 } {
852	    set phi 1.5707963267948966
853	} else { set phi -1.5707963267948966 }
854    } else {
855	if { $n < 0 } {
856	    set rho [expr -$rho] ; set x [expr -$x] ; set y [expr -$y]
857	}
858	set phi [expr 1.0*$rho/[set ${data}(m_4)]]
859	set phi [expr 1.0*([set ${data}(m_3)]-$phi*$phi)/$n]
860	if { abs([set ${data}(m_2)]-abs($phi)) > 1e-7 } {
861	    set phi [LatAngle $phi [set ${data}(m_e)] \
862		              [set ${data}(m_0)]]
863	    if { $phi > 1e19 } {
864		## return 0 180 on error...
865		GMMessage [format $MESS(badinvproj) $TXT(PRJAEA)/$TXT(PRJLEAC)]
866		return [list 0 180]
867	    }
868	} elseif { $phi < 0 } {
869	    set phi -1.5707963267948966
870	} else { set phi 1.5707963267948966 }
871	set lam [expr 1.0*atan2($x,$y)/$n]
872    }
873    set latd [expr $phi*57.29577951308232087684]
874    set longd [expr ($lam*57.29577951308232087684)+[set ${data}(long0)]]
875    return [list $latd [NormalLong $longd]]
876}
877
878proc ProjLEACComputeAux {data datum} {
879    # compute auxiliary parameters from main parameters
880    #  $data is name of global array for the parameters
881    global $data
882
883    set phi2 [expr [set ${data}(lat1)]*0.01745329251994329576]
884    if { [set ${data}(polasp)] == "south" } {
885	set phi1 -1.5707963267948966
886    } else { set phi1 1.5707963267948966 }
887    ComputeAux_AEA_LEAC LEAC $data $datum $phi1 $phi2
888    return
889}
890
891proc ProjLEACPoint {data latd longd datum} {
892
893    return [ProjAEAPoint $data $latd $longd $datum]
894}
895
896proc ProjLEACInvert {data x y} {
897
898    return [ProjAEAInvert $data $x $y]
899}
900
901# ancillary procedures for AEA and LEAC projections
902
903proc ComputeAux_AEA_LEAC {proj data datum phi1 phi2} {
904    global $data MESS
905
906    if { abs($phi1+$phi2) < 1e-10 } {
907	GMMessage $MESS(badProjargs)
908	ProjParams change $proj $data
909	ProjAEAComputeAux $data $datum
910	return
911    }
912    set d [EllipsdData $datum]
913    set ${data}(m_a) [lindex $d 0] ; set f [lindex $d 1]
914    set es [expr $f*(2-$f)]
915    set e [set ${data}(m_e) [expr sqrt($es)]]
916    set oes [set ${data}(m_0) [expr 1-$es]]
917    set n [set ${data}(m_1) [expr sin($phi1)]]
918    set cosp [expr cos($phi1)]
919    set m1 [expr $cosp/sqrt(1.0-$es*$n*$n)]
920    set ml1 [SmallQ $n $e $oes]
921    if { abs($phi1-$phi2) > 1e-10 } {
922	# secant
923	set s2 [expr sin($phi2)]
924	set m2 [expr cos($phi2)/sqrt(1.0-$es*$s2*$s2)]
925	set ml2 [SmallQ $s2 $e $oes]
926	set n [set ${data}(m_1) [expr ($m1*$m1-$m2*$m2)/($ml2-$ml1)]]
927    }
928    set ${data}(m_2) [expr 1.0-0.5*$oes*log((1.0-$e)/(1.0+$e))/$e]
929    set c [set ${data}(m_3) [expr $m1*$m1+$n*$ml1]]
930    set dd [set ${data}(m_4) [expr 1.0/$n]]
931    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
932    set ${data}(m_5) [expr $dd*sqrt($c-$n*[SmallQ sin($phi0) $e $oes])]
933    return
934}
935
936proc LatAngle {qs te toes} {
937    # compute latitude angle phi-1
938
939    set phi [Aasin [expr 0.5*$qs]]
940    if { $te < 1e-7 } { return $phi }
941    for { set i 15 } { $i > 0 } { incr i -1 } {
942	set sinpi [expr sin($phi)] ; set cospi [expr cos($phi)]
943	set con [expr $te*$sinpi]
944	set com [expr 1.0-$con*$con]
945	set dphi [expr 0.5*$com*$com/$cospi* \
946		       ($qs/$toes-$sinpi/$com+0.5/$te* \
947		        log((1.0-$con)/(1.0+$con)))]
948	set phi [expr $phi+$dphi]
949	if { abs($dphi) < 1e-10 } { return $phi }
950    }
951    return 1e20
952}
953
954proc SmallQ {sinphi e oes} {
955    # compute small q
956
957    if { $e < 1e-7 } { return [expr $sinphi+$sinphi] }
958    set con [expr $e*$sinphi]
959    return [expr $oes*($sinphi/(1.0-$con*$con)- \
960	               (0.5/$e)*log((1.0-$con)/(1.0+$con)))]
961}
962
963## Mercator projection: a special case of Lambert Conic Conformal
964#   with the equator as standard parallel
965
966# single standard parallel
967
968proc ProjMerc1ComputeAux {data datum} {
969    # compute auxiliary parameters from main parameters
970    #  $data is name of global array for the parameters
971    global $data
972
973    set d [EllipsdData $datum]
974    set a [lindex $d 0] ; set f [lindex $d 1]
975    set ${data}(m_a) [expr 1.0*$a*[set ${data}(k0)]]
976    set es [expr $f*(2-$f)]
977    set ${data}(m_e) [expr sqrt($es)]
978    set e4 [expr $es*$es] ; set e6 [expr $e4*$es]
979    set ${data}(m_1) [expr $es* \
980	    (0.5+0.20833333333333333333*$es+0.08333333333333333333*$e4+ \
981             0.03611111111111111111*$e6)]
982    set ${data}(m_2) [expr $es* \
983	    (0.14583333333333333333*$es+0.12083333333333333333*$e4+ \
984	     0.07039930555555555555*$e6)]
985    set ${data}(m_3) [expr $es* \
986	    (0.05833333333333333333*$e4+0.07232142857142857142*$e6)]
987    set ${data}(m_4) [expr 0.02653149801587301587*$es*$e6]
988    return
989}
990
991proc ProjMerc1Point {data latd longd datum} {
992    # compute planar Cartesian coordinates for given position
993    global $data
994
995    set prdatum [set ${data}(datum)]
996    if { $datum != $prdatum } {
997	foreach "latd longd" \
998		[ToDatum $latd $longd $datum $prdatum] { break }
999    }
1000    set long0 [set ${data}(long0)]
1001    set longd [NormalLongCentred $longd $long0]
1002    set phi [expr $latd*0.01745329251994329576]
1003    set se [expr sin($phi)*[set ${data}(m_e)]]
1004    set mm [expr tan((1.5707963267948966+$phi)/2.0)* \
1005     	         pow((1.0-$se)/(1.0+$se), [set ${data}(m_e)]/2.0)]
1006    set x [expr [set ${data}(m_a)]*($longd-$long0)*0.01745329251994329576]
1007    set y [expr [set ${data}(m_a)]*log($mm)]
1008    return [list $x $y]
1009}
1010
1011proc ProjMerc1Invert {data x y} {
1012    # return list with latitude and longitude for point at $x,$y in the terrain
1013    global $data
1014
1015    set ksi [expr 1.5707963267948966- \
1016	          2*atan(exp(-$y/[set ${data}(m_a)]))]
1017    set latd [expr ($ksi+[set ${data}(m_1)]*sin($ksi+$ksi)+ \
1018	            [set ${data}(m_2)]*sin(4*$ksi)+ \
1019	            [set ${data}(m_3)]*sin(6*$ksi)+ \
1020	            [set ${data}(m_4)]*sin(8*$ksi))* \
1021		   57.29577951308232087684]
1022    set longd [expr $x/[set ${data}(m_a)]*57.29577951308232087684+ \
1023	            [set ${data}(long0)]]
1024    return [list $latd [NormalLong $longd]]
1025}
1026
1027# two standard parallels (symmetrical with respect to the equator)
1028
1029proc ProjMerc2ComputeAux {data datum} {
1030    # compute auxiliary parameters from main parameters
1031    #  $data is name of global array for the parameters
1032    global $data
1033
1034    set d [EllipsdData $datum]
1035    set a [lindex $d 0] ; set f [lindex $d 1]
1036    set es [expr $f*(2-$f)]
1037    set ${data}(m_e) [expr sqrt($es)]
1038    set phi1 [expr abs([set ${data}(lat1)])*0.01745329251994329576]
1039    set s [expr sin($phi1)]
1040    set k0 [expr cos($phi1)/sqrt(1-$es*$s*$s)]
1041    set ${data}(m_a) [expr 1.0*$a*$k0]
1042    set e4 [expr $es*$es] ; set e6 [expr $e4*$es]
1043    set ${data}(m_1) [expr $es* \
1044	    (0.5+0.20833333333333333333*$es+0.08333333333333333333*$e4+ \
1045             0.03611111111111111111*$e6)]
1046    set ${data}(m_2) [expr $es* \
1047	    (0.14583333333333333333*$es+0.12083333333333333333*$e4+ \
1048	     0.07039930555555555555*$e6)]
1049    set ${data}(m_3) [expr $es* \
1050	    (0.05833333333333333333*$e4+0.07232142857142857142*$e6)]
1051    set ${data}(m_4) [expr 0.02653149801587301587*$es*$e6]
1052    return
1053}
1054
1055proc ProjMerc2Point {data latd longd datum} {
1056    # compute planar Cartesian coordinates for given position
1057    global $data
1058
1059    return [ProjMerc1Point $data $latd $longd $datum]
1060}
1061
1062proc ProjMerc2Invert {data x y} {
1063    # return list with latitude and longitude for point at $x,$y in the terrain
1064    global $data
1065
1066    return [ProjMerc1Invert $data $x $y]
1067}
1068
1069## Mercator projection: a special case of Lambert Conic Conformal
1070#   with the equator as standard parallel
1071
1072# single standard parallel, spherical case
1073# adapted from libproj 4.6.1
1074
1075proc ProjSphMercComputeAux {data datum} {
1076    # compute auxiliary parameters from main parameters
1077    #  $data is name of global array for the parameters
1078    global $data
1079
1080    set a [lindex [EllipsdData $datum] 0]
1081    set lat0 [set ${data}(lat0)]
1082    if { $lat0 > 1e-15 } {
1083	set k0 [expr cos($lat0*0.01745329251994329576)]
1084    } else { set k0 1 }
1085    set ${data}(m_a) [expr 1.0*$a*$k0]
1086    return
1087}
1088
1089proc ProjSphMercPoint {data latd longd datum} {
1090    # compute planar Cartesian coordinates for given position
1091    global $data
1092
1093    set prdatum [set ${data}(datum)]
1094    if { $datum != $prdatum } {
1095	foreach "latd longd" \
1096		[ToDatum $latd $longd $datum $prdatum] { break }
1097    }
1098    set long0 [set ${data}(long0)]
1099    set longd [NormalLongCentred $longd $long0]
1100    set x [expr [set ${data}(m_a)]*($longd-$long0)*0.01745329251994329576]
1101    set y [expr [set ${data}(m_a)]* \
1102	       log(tan(0.78539816339744830962+$latd*0.00872664625997164788))]
1103    return [list $x $y]
1104}
1105
1106proc ProjSphMercInvert {data x y} {
1107    # return list with latitude and longitude for point at $x,$y in the terrain
1108    global $data
1109
1110    set latd [expr (1.5707963267948966-2*atan(exp(-$y/[set ${data}(m_a)]))) * \
1111		57.29577951308232087684]
1112    set longd [expr $x/[set ${data}(m_a)]*57.29577951308232087684+ \
1113	            [set ${data}(long0)]]
1114    return [list $latd [NormalLong $longd]]
1115}
1116
1117## Cassini-Soldner projection
1118
1119proc ProjCSComputeAux {data datum} {
1120    # compute auxiliary parameters from main parameters
1121    #  $data is name of global array for the parameters
1122    global $data
1123
1124    set d [EllipsdData $datum]
1125    set a [set ${data}(m_a) [lindex $d 0]] ; set f [lindex $d 1]
1126    set es [set ${data}(m_e) [expr $f*(2-$f)]]
1127    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
1128    set e4 [expr $es*$es] ; set e6 [expr $es*$e4]
1129    set m1 [set ${data}(m_1) [expr 1-0.25*$es-0.046875*$e4-0.01953125*$e6]]
1130    set m2 [set ${data}(m_2) [expr 0.375*$es+0.09375*$e4+0.0439453125*$e6]]
1131    set m3 [set ${data}(m_3) [expr 0.05859375*$e4+0.0439453125*$e6]]
1132    set m4 [set ${data}(m_4) [expr 0.01139322916666666667*$e6]]
1133    set ${data}(m_0) [expr $a*($m1*$phi0-$m2*sin($phi0+$phi0)+ \
1134	                          $m3*sin(4*$phi0)-$m4*sin(6*$phi0))]
1135    set se [expr sqrt(1-$es)]
1136    set e1 [expr (1.0-$se)/(1+$se)] ; set e12 [expr $e1*$e1]
1137    set ${data}(m_5) [expr $e1*(1.5-0.84375*$e12)]
1138    set ${data}(m_6) [expr $e12*(1.3125-1.71875*$e12)]
1139    set ${data}(m_7) [expr 1.57291666666666666667*$e12*$e1]
1140    set ${data}(m_8) [expr 2.142578125*$e12*$e12]
1141    return
1142}
1143
1144proc ProjCSPoint {data latd longd datum} {
1145    # compute planar Cartesian coordinates for given position
1146    global $data
1147
1148    set prdatum [set ${data}(datum)]
1149    if { $datum != $prdatum } {
1150	foreach "latd longd" \
1151		[ToDatum $latd $longd $datum $prdatum] { break }
1152    }
1153    set long0 [set ${data}(long0)]
1154    set longd [NormalLongCentred $longd $long0]
1155    set a [set ${data}(m_a)] ; set es [set ${data}(m_e)]
1156    set phi [expr $latd*0.01745329251994329576]
1157    set cphi [expr cos($phi)] ; set sphi [expr sin($phi)]
1158    set A [expr ($longd-$long0)*0.01745329251994329576*$cphi]
1159    set A2 [expr $A*$A] ; set A4 [expr $A2*$A2]
1160    set t [expr tan($phi)] ; set T [expr $t*$t]
1161    set C [expr $es*$cphi*$cphi/(1.0-$es)]
1162    set v [expr $a/sqrt(1-$es*$sphi*$sphi)]
1163    set M [expr $a*([set ${data}(m_1)]*$phi- \
1164	            [set ${data}(m_2)]*sin($phi+$phi)+ \
1165	            [set ${data}(m_3)]*sin(4*$phi)- \
1166                    [set ${data}(m_4)]*sin(6*$phi))]
1167    set x [expr $v*$A*(1-$T*$A2/6.0-(8-$T+8*$C)*$T*$A4/120.0)]
1168    set y [expr $M-[set ${data}(m_0)]+$v*$t*(0.5*$A2+(5-$T+6*$C)*$A4/24.0)]
1169    return [list $x $y]
1170}
1171
1172proc ProjCSInvert {data x y} {
1173    # return list with latitude and longitude for point at $x,$y in the terrain
1174    global $data
1175
1176    set a [set ${data}(m_a)] ; set es [set ${data}(m_e)]
1177    set miu1 [expr ([set ${data}(m_0)]+$y)/($a*[set ${data}(m_1)])]
1178    set phi1 [expr $miu1+[set ${data}(m_5)]*sin($miu1+$miu1)+ \
1179	           [set ${data}(m_6)]*sin(4*$miu1)+ \
1180	           [set ${data}(m_7)]*sin(6*$miu1)+ \
1181	           [set ${data}(m_8)]*sin(8*$miu1)]
1182    set sphi1 [expr sin($phi1)]
1183    set t [expr 1-$es*$sphi1*$sphi1]
1184    set v1 [expr $a/sqrt($t)] ; set rho1 [expr $v1*(1-$es)/$t]
1185    set t [expr tan($phi1)] ; set T1 [expr $t*$t]
1186    set D [expr $x/$v1] ; set D2 [expr $D*$D]
1187    set tt [expr (1+3*$T1)*$D2]
1188    set lat [expr ($phi1-$v1*$t/$rho1*$D2*(0.5-$tt/24.0))* \
1189	          57.29577951308232087684]
1190    set long [expr $D*(1-$T1*$D2/3.0+$tt*$D2/15.0)/cos($phi1)* \
1191	          57.29577951308232087684+[set ${data}(long0)]]
1192    return [list $lat [NormalLong $long]]
1193}
1194
1195## American polyconic projection
1196# adapted from PROJ4.0
1197
1198proc ProjAPOLYComputeAux {data datum} {
1199    # compute auxiliary parameters from main parameters
1200    #  $data is name of global array for the parameters
1201    global $data
1202
1203    set d [EllipsdData $datum]
1204    set ${data}(m_a) [lindex $d 0]
1205    set f [lindex $d 1]
1206    set es [set ${data}(m_5) [expr $f*(2-$f)]]
1207    MeridDistParams $data $es
1208    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
1209    set ${data}(m_6) \
1210	    [MeridionalDist $data $phi0 [expr sin($phi0)] [expr cos($phi0)]]
1211    return
1212}
1213
1214proc ProjAPOLYPoint {data latd longd datum} {
1215    # compute planar Cartesian coordinates for given position
1216    global $data
1217
1218    set prdatum [set ${data}(datum)]
1219    if { $datum != $prdatum } {
1220	foreach "latd longd" \
1221		[ToDatum $latd $longd $datum $prdatum] { break }
1222    }
1223    set a [set ${data}(m_a)]
1224    set phi [expr $latd*0.01745329251994329576]
1225    set lam [expr $longd*0.01745329251994329576]
1226    if { abs($phi) <= 1e-10 } {
1227	set x [expr $lam*$a] ; set y [expr -$a*[set ${data}(m_6)]]
1228    } else {
1229	set sp [expr sin($phi)]
1230	if { abs([set cp [expr cos($phi)]]) > 1e-10 } {
1231	    set ms [expr $cp/sqrt(1-[set ${data}(m_5)]*$sp*$sp)/$sp]
1232	} else { set ms 0 }
1233	set t [expr $lam*$sp]
1234	set x [expr $a*$ms*sin($t)]
1235	set y [expr $a*(([MeridionalDist $data $phi $sp $cp]- \
1236		[set ${data}(m_6)])+$ms*(1-cos($t)))]
1237    }
1238    return [list $x $y]
1239}
1240
1241proc ProjAPOLYInvert {data x y} {
1242    # return list with latitude and longitude for point at $x,$y in the terrain
1243    global $data MESS TXT
1244
1245    set a [set ${data}(m_a)]
1246    set x [expr 1.0*$x/$a] ; set y [expr 1.0*$y/$a]
1247    set y [expr $y+[set ${data}(m_6)]]
1248    if { abs($y) <= 1e-10 } {
1249	set phi 0 ; set lam $x
1250    } else {
1251	set r [expr $x*$x+$y*$y] ; set phi $y
1252	set es [set ${data}(m_5)]
1253	for { set i 0 } { $i < 20 } { incr i } {
1254	    set sp [expr sin($phi)] ; set cp [expr cos($phi)]
1255	    set s2ph [expr $sp*$cp]
1256	    if { abs($cp) < 1e-12 } {
1257		GMMessage [format $MESS(badinvproj) $TXT(PRJAPOLY)]
1258		set i 0
1259		break
1260	    }
1261	    set mlp [expr sqrt(1-$es*$sp*$sp)]
1262	    set c [expr $sp*$mlp/$cp]
1263	    set ml [MeridionalDist $data $phi $sp $cp]
1264	    set mlb [expr $ml*$ml+$r]
1265	    set mlp [expr (1-$es)/($mlp*$mlp*$mlp)]
1266	    set dphi [expr ($ml+$ml+$c*$mlb-2.0*$y*($c*$ml+1))/ \
1267		    ($es*$s2ph*($mlb-2*$y*$ml)/$c+ \
1268		     2*($y-$ml)*($c*$mlp-1.0/$s2ph)-$mlp-$mlp)]
1269	    set phi [expr $phi+$dphi]
1270	    if { abs($dphi) <= 1e-12 } { break }
1271	}
1272	if { $i == 0 } {
1273	    GMMessage [format $MESS(badinvproj) $TXT(PRJAPOLY)]
1274	}
1275	set c [expr sin($phi)]
1276	set lam [expr [Aasin [expr $x*tan($phi)*sqrt(1-$es*$c*$c)]]/sin($phi)]
1277    }
1278    set longd [expr $lam*57.29577951308232087684]
1279    set latd [expr $phi*57.29577951308232087684]
1280    return [list $latd $longd]
1281}
1282
1283## Equidistant Cylindrical projection
1284# adapted from libproj 4.6.1
1285
1286proc ProjEqCylComputeAux {data datum} {
1287    # compute auxiliary parameters from main parameters
1288    #  $data is name of global array for the parameters
1289    global $data
1290
1291    set a [set ${data}(m_a) [lindex [EllipsdData $datum] 0]]
1292    set ${data}(m_0) [expr cos([set ${data}(lat0)]*0.01745329251994329576)*$a]
1293    return
1294}
1295
1296proc ProjEqCylPoint {data latd longd datum} {
1297    # compute planar Cartesian coordinates for given position
1298    global $data
1299
1300    set prdatum [set ${data}(datum)]
1301    if { $datum != $prdatum } {
1302	foreach "latd longd" \
1303		[ToDatum $latd $longd $datum $prdatum] { break }
1304    }
1305    set x [expr [set ${data}(m_0)]*$longd*0.01745329251994329576]
1306    set y [expr [set ${data}(m_a)]*($latd-[set ${data}(lat0)])* \
1307	       0.01745329251994329576]
1308    return [list $x $y]
1309}
1310
1311proc ProjEqCylInvert {data x y} {
1312    # return list with latitude and longitude for point at $x,$y in the terrain
1313    global $data
1314
1315    set latd [expr $y/[set ${data}(m_a)]*57.29577951308232087684+ \
1316		  [set ${data}(lat0)]]
1317    set longd [expr $x/[set ${data}(m_0)]*57.29577951308232087684]
1318    return [list $latd [NormalLong $longd]]
1319}
1320
1321## auxiliary procs for polyconic projection
1322# meridional distance for ellipsoid and inverse
1323
1324proc MeridDistParams {data es} {
1325    # compute parameters for meridional distance and inverse functions
1326    #  saving them at indices m_0, ..., m_4 of array $data
1327    global $data
1328
1329    set t [expr $es*(0.046875+$es*(0.01953125+$es*0.01068115234375))]
1330    set ${data}(m_0) [expr 1-$es*(0.25+$t)]
1331    set ${data}(m_1) [expr $es*(0.75-$t)]
1332    set t [expr $es*$es]
1333    set ${data}(m_2) [expr $t*(0.46875-$es* \
1334	    (0.01302083333333333333+$es*0.00712076822916666666))]
1335    set t [expr $es*$t]
1336    set ${data}(m_3) \
1337	    [expr $t*(0.36458333333333333333-$es*0.00569661458333333333)]
1338    set ${data}(m_4) [expr $t*$es*0.3076171875]
1339    return
1340}
1341
1342proc MeridionalDist {data phi sphi cphi} {
1343    # compute meridional distance for given latitude
1344    # assume parameters for transformation are in array $data at indices
1345    #  m_0, ..., m_4
1346    global $data
1347
1348    set cphi [expr $cphi*$sphi]
1349    set sphi [expr $sphi*$sphi]
1350    return [expr [set ${data}(m_0)]*$phi-$cphi* \
1351	    ([set ${data}(m_1)]+$sphi* \
1352	     ([set ${data}(m_2)]+$sphi* \
1353	      ([set ${data}(m_3)]+$sphi*[set ${data}(m_4)])))]
1354}
1355
1356proc LatFromMeridDist {data d} {
1357    # compute latitude from given meridional distance
1358    # assume parameters for transformation and e^2 are in array $data at
1359    #  indices m_0, ..., m_4 and m_5
1360    global $data MESS
1361
1362    set es [set ${data}(m_5)]
1363    set k [expr 1.0/(1-$es)]
1364    set phi $d
1365    for { set i 0 } { $i < 10 } { incr i } {
1366	set s [expr sin($phi)]
1367	set t [expr 1.0-$es*$s*$s]
1368	set dphi [expr [MeridionalDist $data $phi $s [expr cos($phi)]]* \
1369		($t*sqrt($t))*$k]
1370	set phi [expr $phi-$dphi]
1371	if { abs($dphi) < 1e-11 } { return $phi }
1372    }
1373    GMMessage $MESS(badinvmdist)
1374    return $phi
1375}
1376
1377## Stereographic projection
1378#
1379
1380proc ProjStereogrComputeAux {data datum} {
1381    # compute auxiliary parameters from main parameters
1382    #  $data is name of global array for the parameters
1383    global $data
1384
1385    set d [EllipsdData $datum]
1386    set ${data}(m_a) [lindex $d 0] ; set f [lindex $d 1]
1387    set e [set ${data}(m_e) [expr sqrt($f*(2-$f))]]
1388    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
1389    set k0 [set ${data}(k0)]
1390    set t [expr abs($phi0)]
1391    if { abs($t-1.5707963267948966) < 1e-10 } {
1392	if { $phi0 < 0 } {
1393	    set mode s_pole
1394	} else { set mode n_pole }
1395	set ${data}(m_0) [expr 2*$k0/sqrt(pow(1+$e,1+$e)*pow(1-$e,1-$e))]
1396    } elseif { $t > 1e-10 } {
1397	set mode obliq
1398	set t [expr sin($phi0)]
1399	set X [expr 2*atan([SSFunc_ $phi0 $t $e])-1.5707963267948966]
1400	set t [expr $t*$e]
1401	set ${data}(m_0) [expr 2*$k0*cos($phi0)/sqrt(1.0-$t*$t)]
1402	set ${data}(m_1) [expr sin($X)]
1403	set ${data}(m_2) [expr cos($X)]
1404    } else {
1405	set mode equat
1406	set ${data}(m_0) [expr 2*$k0]
1407	set ${data}(m_1) 0
1408	set ${data}(m_2) 1.0
1409    }
1410    set ${data}(m_3) $mode
1411    return
1412}
1413
1414proc ProjStereogrPoint {data latd longd datum} {
1415    # compute planar Cartesian coordinates for given position
1416    global $data
1417
1418    set prdatum [set ${data}(datum)]
1419    if { $datum != $prdatum } {
1420	foreach "latd longd" \
1421		[ToDatum $latd $longd $datum $prdatum] { break }
1422    }
1423    set long0 [set ${data}(long0)]
1424    set longd [NormalLongCentred $longd $long0]
1425    set e [set ${data}(m_e)]
1426    set akm1 [set ${data}(m_0)]
1427    set phi [expr $latd*0.01745329251994329576]
1428    set lam [expr ($longd-$long0)*0.01745329251994329576]
1429    set coslam [expr cos($lam)] ; set sinlam [expr sin($lam)]
1430    set sinphi [expr sin($phi)]
1431    switch [set ${data}(m_3)] {
1432	obliq -
1433	equat {
1434	    set X [expr 2*atan([SSFunc_ $phi $sinphi $e])-1.5707963267948966]
1435	    set sinX [expr sin($X)] ; set cosX [expr cos($X)]
1436	    set sinX1 [set ${data}(m_1)]
1437	    set cosX1 [set ${data}(m_2)]
1438	    set A [expr $akm1/($cosX1*(1+$sinX1*$sinX+$cosX1*$cosX*$coslam))]
1439	    set y [expr $A*($cosX1*$sinX-$sinX1*$cosX*$coslam)]
1440	    set x [expr $A*$cosX]
1441	}
1442	n_pole {
1443	    set x [expr $akm1*[ExpRedLat $phi $sinphi $e]]
1444	    set y [expr -$x*$coslam]
1445	}
1446	s_pole {
1447	    set x [expr $akm1*[ExpRedLat [expr -$phi] [expr -$sinphi] $e]]
1448	    set y [expr $x*$coslam]
1449	}
1450    }
1451    set a [set ${data}(m_a)]
1452    set x [expr $x*$sinlam*$a] ; set y [expr $y*$a]
1453    return [list $x $y]
1454}
1455
1456proc ProjStereogrInvert {data x y} {
1457    # return list with latitude and longitude for point at $x,$y in the terrain
1458    global $data MESS TXT
1459
1460    set e [set ${data}(m_e)] ; set a [set ${data}(m_a)]
1461    set akm1 [set ${data}(m_0)]
1462    set long0 [set ${data}(long0)] ; set datum [set ${data}(datum)]
1463    set x [expr $x/$a] ; set y [expr $y/$a]
1464    if { abs($x)<1e-6 && abs($y)<1e-6 } {
1465	return [list [set ${data}(lat0)] $long0]
1466    }
1467    set rho [Hypot $x $y]
1468    set spole 0
1469    switch [set m [set ${data}(m_3)]] {
1470	obliq -
1471	equat {
1472	    set sinX1 [set ${data}(m_1)]
1473	    set cosX1 [set ${data}(m_2)]
1474	    set tp [expr 2*atan2($rho*$cosX1,$akm1)]
1475	    set cosphi [expr cos($tp)] ; set sinphi [expr sin($tp)]
1476	    if { $rho == 0.0 } {
1477		set phi_l [Aasin [expr $cosphi*$sinX1]]
1478	    } else {
1479		set phi_l [Aasin [expr $cosphi*$sinX1+($y*$sinphi*$cosX1/$rho)]]
1480	    }
1481	    set tp [expr tan(0.5*(1.5707963267948966+$phi_l))]
1482	    set x [expr $x*$sinphi]
1483	    set y [expr $rho*$cosX1*$cosphi-$y*$sinX1*$sinphi]
1484	    set halfpi 1.5707963267948966
1485	    set halfe [expr 0.5*$e]
1486	}
1487	n_pole {
1488	    set y [expr -$y]
1489	    set tp [expr -1.0*$rho/$akm1]
1490	    set phi_l [expr 1.5707963267948966-2*atan($tp)]
1491	    set halfpi -1.5707963267948966
1492	    set halfe [expr -0.5*$e]
1493	}
1494	s_pole {
1495	    set spole 1
1496	    set tp [expr -1.0*$rho/$akm1]
1497	    set phi_l [expr 1.5707963267948966-2*atan($tp)]
1498	    set halfpi -1.5707963267948966
1499	    set halfe [expr -0.5*$e]
1500	}
1501    }
1502    for { set i 8 } { [incr i -1] } { set phi_l $phi } {
1503	set sinphi [expr $e*sin($phi_l)]
1504	set phi [expr 2*atan($tp*pow((1+$sinphi)/(1-$sinphi),$halfe))-$halfpi]
1505	if { abs($phi_l-$phi) < 1e-10 } {
1506	    if { $spole } { set phi [expr -$phi] }
1507	    if { $x == 0 && $y == 0 } {
1508		set longd 0
1509	    } else { set longd [expr atan2($x,$y)*57.29577951308232087684] }
1510	    set latd [expr $phi*57.29577951308232087684]
1511	    set longd [expr $longd+$long0]
1512	    return [list $latd $longd]
1513	}
1514    }
1515    ## return 0 180 on error...
1516    GMMessage [format $MESS(badinvproj) $TXT(PRJStereogr)]
1517    return [list 0 180]
1518}
1519
1520# Schreiber double projection: used for the RDG (Netherlands) grid
1521#  this is a special case of the oblique Stereographic projection
1522#
1523
1524proc ProjSchreiberComputeAux {data datum} {
1525    # compute auxiliary parameters from main parameters
1526    #  $data is name of global array for the parameters
1527    global $data
1528
1529    set phi0 [set ${data}(m_0) 0.91029672689323934677]
1530    set lam0 [set ${data}(m_1) 0.09403203751960005358]
1531    set e [set ${data}(m_2) 0.08169683122252750299]
1532    set k0 0.9999079
1533    # 2*k*R
1534    set ${data}(m_3) [expr 12765289.142*$k0]
1535    set n 1.00047585668 ; set m 0.003773953832
1536    set tau0 [expr log([SSFunc_ $phi0 [expr sin($phi0)] $e])]
1537    set B0 [expr 2*atan(pow(2.7182818284590452354,$n*$tau0+$m))- \
1538	    1.5707963267948966]
1539    set ${data}(m_4) [expr sin($B0)]
1540    set ${data}(m_5) [expr cos($B0)]
1541    return
1542}
1543
1544proc ProjSchreiberPoint {data latd longd datum} {
1545    # compute planar Cartesian coordinates for given position
1546    global $data
1547
1548    set prdatum [set ${data}(datum)]
1549    if { $datum != $prdatum } {
1550	foreach "latd longd" \
1551		[ToDatum $latd $longd $datum $prdatum] { break }
1552    }
1553    # just for the sake of those not knowing that this projection
1554    #  was designed for use in the Netherlands...
1555    set long0 5.38763888888888888564
1556    set longd [NormalLongCentred $longd $long0]
1557    set n 1.00047585668 ; set m 0.003773953832
1558    set e [set ${data}(m_2)]
1559    set phi [expr $latd*0.01745329251994329576]
1560    set lam [expr $n*($longd*0.01745329251994329576-[set ${data}(m_1)])]
1561    set coslam [expr cos($lam)] ; set sinlam [expr sin($lam)]
1562    set sinphi [expr sin($phi)]
1563    set tau [expr log([SSFunc_ $phi $sinphi $e])]
1564    set B [expr 2*atan(pow(2.7182818284590452354,$n*$tau+$m))- \
1565	    1.5707963267948966]
1566    set sinB [expr sin($B)] ; set cosB [expr cos($B)]
1567    set sinB0 [set ${data}(m_4)]
1568    set cosB0 [set ${data}(m_5)]
1569    set tkR [set ${data}(m_3)]
1570    set A [expr $tkR/(1+$sinB*$sinB0+$cosB*$cosB0*$coslam)]
1571    set x [expr $A*$sinlam*$cosB]
1572    set y [expr $A*($sinB*$cosB0-$cosB*$sinB0*$coslam)]
1573    return [list $x $y]
1574}
1575
1576proc ProjSchreiberInvert {data x y} {
1577    # return list with latitude and longitude for point at $x,$y in the terrain
1578    global $data MESS TXT
1579
1580    set datum [set ${data}(datum)]
1581    set long0 [expr [set ${data}(m_1)]*57.29577951308232087684]
1582    if { abs($x)<1e-6 && abs($y)<1e-6 } {
1583	set lat0 [expr [set ${data}(m_0)]*57.29577951308232087684]
1584	return [list $lat0 $long0]
1585    }
1586    set n 1.00047585668 ; set m 0.003773953832
1587    set sinB0 [set ${data}(m_4)]
1588    set cosB0 [set ${data}(m_5)]
1589    set tkR [set ${data}(m_3)]
1590    set rho [Hypot $x $y]
1591    set Psi [expr 2*atan2($rho,$tkR)]
1592    set sinPsi [expr sin($Psi)]
1593    set B [Aasin [expr cos($Psi)*$sinB0+$y*$sinPsi*$cosB0/$rho]]
1594    set lam [Aasin [expr $x*$sinPsi/$rho/cos($B)]]
1595    set longd [expr $lam/$n*57.29577951308232087684+$long0]
1596    set tau [expr (log(tan(0.5*(1.5707963267948966+$B)))-$m)/$n]
1597    set etau [expr pow(2.7182818284590452354,$tau)]
1598    set phi [expr 2*atan($etau)-1.5707963267948966]
1599    set e [set ${data}(m_2)] ; set he [expr 0.5*$e]
1600    for { set i 0 } { $i < 50 } { incr i } {
1601	set esp [expr $e*sin($phi)]
1602	set phi_ [expr 2*atan($etau*pow((1+$esp)/(1-$esp),$he))- \
1603		1.5707963267948966]
1604	if { abs($phi_-$phi) < 1e-10 } {
1605	    set latd [expr $phi*57.29577951308232087684]
1606	    return [list $latd $longd]
1607	}
1608	set phi $phi_
1609    }
1610    ## return 0 180 on error...
1611    GMMessage [format $MESS(badinvproj) $TXT(PRJSchreiber)]
1612    return [list 0 180]
1613}
1614
1615## ancillary procs for Stereographic and Schreiber projections
1616
1617proc SSFunc_ {phit sinphi eccen} {
1618
1619    set sinphi [expr $sinphi*$eccen]
1620    return [expr tan(0.5*(1.5707963267948966+$phit))* \
1621	    pow((1-$sinphi)/(1+$sinphi),0.5*$eccen)]
1622}
1623
1624# Swiss Oblique Mercator projection
1625
1626proc ProjSOMComputeAux {data datum} {
1627    # compute auxiliary parameters from main parameters
1628    #  $data is name of global array for the parameters
1629    global $data
1630
1631    set d [EllipsdData $datum]
1632    set a [lindex $d 0] ; set f [lindex $d 1]
1633    set es [expr $f*(2-$f)]
1634    set e [set ${data}(m_e) [expr sqrt($es)]]
1635    set he [set ${data}(m_0) [expr 0.5*$e]]
1636    set oes [expr 1-$es]
1637    set roes [set ${data}(m_6) [expr 1.0/$oes]]
1638    set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576]
1639    set cp [expr cos($phi0)] ; set cp [expr $cp*$cp]
1640    set c [set ${data}(m_1) [expr sqrt(1+$es*$cp*$cp*$roes)]]
1641    set sp [expr sin($phi0)]
1642    set sp0 [set ${data}(m_3) [expr $sp/$c]]
1643    set p0 [Aasin $sp0]
1644    set ${data}(m_2) [expr cos($p0)]
1645    set sp [expr $sp*$e]
1646    set ${data}(m_4) [expr log(tan(0.78539816339744833+0.5*$p0))- \
1647	    $c*(log(tan(0.78539816339744833+0.5*$phi0))- \
1648	        $he*log((1.0+$sp)/(1.0-$sp)))]
1649    set ${data}(m_5) [expr $a*[set ${data}(k0)]*sqrt($oes)/(1.0-$sp*$sp)]
1650    return
1651}
1652
1653proc ProjSOMPoint {data latd longd datum} {
1654    # compute planar Cartesian coordinates for given position
1655    global $data
1656
1657    set prdatum [set ${data}(datum)]
1658    if { $datum != $prdatum } {
1659	foreach "latd longd" \
1660		[ToDatum $latd $longd $datum $prdatum] { break }
1661    }
1662    # just for the sake of those not knowing that this projection
1663    #  was designed for use in Switzerland...
1664    set long0 [set ${data}(long0)]
1665    set longd [NormalLongCentred $longd $long0]
1666    set phi [expr $latd*0.01745329251994329576]
1667    set sp [expr [set ${data}(m_e)]*sin($phi)]
1668    set c [set ${data}(m_1)]
1669    set phip [expr 2*atan(exp($c*(log(tan(0.78539816339744833+0.5*$phi))- \
1670	    [set ${data}(m_0)]*log((1.0+$sp)/(1.0-$sp)))+ \
1671	    [set ${data}(m_4)]))-1.5707963267948966]
1672    set lam [expr $c*($longd-$long0)*0.01745329251994329576]
1673    set cp [expr cos($phip)]
1674    set phipp [Aasin [expr [set ${data}(m_2)]*sin($phip)- \
1675	    [set ${data}(m_3)]*$cp*cos($lam)]]
1676    set akR [set ${data}(m_5)]
1677    set x [expr $akR*[Aasin [expr $cp*sin($lam)/cos($phipp)]]]
1678    set y [expr $akR*log(tan(0.78539816339744833+0.5*$phipp))]
1679    return [list $x $y]
1680}
1681
1682proc ProjSOMInvert {data x y} {
1683    # return list with latitude and longitude for point at $x,$y in the terrain
1684    global $data MESS TXT
1685
1686    set datum [set ${data}(datum)]
1687    set akR [set ${data}(m_5)]
1688    set phipp [expr 2*(atan(exp(1.0*$y/$akR))-0.78539816339744833)]
1689    set lampp [expr 1.0*$x/$akR]
1690    set cp [expr cos($phipp)]
1691    set phip [Aasin [expr [set ${data}(m_2)]*sin($phipp)+ \
1692	    [set ${data}(m_3)]*$cp*cos($lampp)]]
1693    set lamp [Aasin [expr $cp*sin($lampp)/cos($phip)]]
1694    set e [set ${data}(m_e)]
1695    set he [set ${data}(m_0)] ; set c [set ${data}(m_1)]
1696    set roes [set ${data}(m_6)]
1697    set con [expr ([set ${data}(m_4)]- \
1698	    log(tan(0.78539816339744833+0.5*$phip)))/$c]
1699    set i 10
1700    while { [incr i -1] } {
1701	set esp [expr $e*sin($phip)]
1702	set delp [expr ($con+log(tan(0.78539816339744833+0.5*$phip))- \
1703		$he*log((1.0+$esp)/(1.0-$esp)))*(1.0-$esp*$esp)*cos($phip)* \
1704		$roes]
1705	set phip [expr $phip-$delp]
1706	if { abs($delp) < 1e-10 } { break }
1707    }
1708    if { $i } {
1709	set latd [expr $phip*57.29577951308232087684]
1710	set longd [expr 57.29577951308232087684*$lamp/$c+ \
1711		[set ${data}(long0)]]
1712	return [list $latd $longd]
1713    }
1714    ## return 0 180 on error...
1715    GMMessage [format $MESS(badinvproj) $TXT(PRJSOM)]
1716    return [list 0 180]
1717}
1718
1719# ancillary procedures for Swiss Oblique Mercator and others
1720
1721proc Aasin {v} {
1722    # arcsin function
1723    global MESS
1724
1725    set av [expr abs($v)]
1726    if { $av >= 1 } {
1727	if { $av > 1.00000000000001 } {
1728	    GMMessage [format $MESS(badargtofunc) asin]
1729	}
1730	if { $v < 0 } {
1731	    return -1.5707963267948966
1732	}
1733	return 1.5707963267948966
1734    }
1735    return [expr asin($v)]
1736}
1737
1738# ancillary procedures for LCC1, LCC2, Stereographic and Schreiber projections
1739
1740proc Hypot {x y} {
1741    # computation of sqrt($x*$x+$y*$y) avoiding overflows
1742
1743    if { $x < 0 } {
1744	set x [expr -$x]
1745    } elseif { $x == 0 } {
1746	if { $y < 0 } { return [expr -$y] }
1747	return $y
1748    }
1749    if { $y < 0 } {
1750	set y [expr -$y]
1751    } elseif { $y == 0 } { return $x }
1752    if { $x < $y } {
1753	set x [expr 1.0*$x/$y]
1754	return [expr $y*sqrt(1+$x*$x)]
1755    }
1756    set y [expr 1.0*$y/$x]
1757    return [expr $x*sqrt(1+$y*$y)]
1758}
1759
1760#
1761# SL contribution
1762#
1763# Implementation of the EOV (Hungarian National Projection)
1764#
1765# 2008-07-15
1766#
1767# Based on publications by Jozsef Varga Phd.
1768#   http://www.agt.bme.hu/staff_h/varga/Osszes/Dok3uj.htm
1769#   http://www.agt.bme.hu/staff_h/varga/gps/kezdoknek.html
1770#
1771
1772proc ProjEOVPoint {data latd longd datum} {
1773    # compute planar Cartesian coordinates for given position
1774    # although this projection computes x for northing and
1775    #  y for easting, for compatibility with the mapping procedures
1776    #  this procedure returns the usual (easting, northing) pair,
1777    #  and the grid coordinates are computed by inversing it
1778    #  $data is not used
1779
1780    # MF change: use the 3 parameter conversion to the HD72 if needed
1781    if { $datum != "Hungarian Datum 1972" } {
1782	foreach "latd longd" \
1783		[ToDatum $latd $longd $datum "Hungarian Datum 1972"] { break }
1784    }
1785    # MF changes after this point: for efficiency, original code commented out
1786    # HD72 Phi, Lambda to radians
1787    set Phi [expr $latd*0.01745329251994329576]
1788    set Lambda [expr $longd*0.01745329251994329576]
1789    #---
1790
1791    # HD72 Phi, Lambda to spherical phi, lambda
1792    ## GRS67 semi-major / semi-minor axis
1793    # set a 6378160.0
1794    # set b 6356774.516
1795    ## EOV parameters
1796    set k 1.003110007693
1797    set n 1.000719704936
1798    # set Lambda_0 [expr (19.0 + 02.0/60.0 + 54.8584/3600.0) / 360.0 * 2.0*$pi]
1799    # set Lambda_0 0.33246029532469185632
1800    ## GRS67 eccentricity
1801    # set epsilon [expr sqrt(($a*$a-$b*$b)/($a*$a))]
1802    set epsilon 0.0818205680555
1803    # set A [expr pow(tan($pi/4.0+$Phi/2.0), $n)]
1804    set A [expr pow(tan(0.785398163397448+0.5*$Phi), $n)]
1805    # set B [expr pow((1.0-$epsilon*sin($Phi))/(1.0+$epsilon*sin($Phi)),
1806    #                 $n*$epsilon/2.0)]
1807    set sinPhi [expr sin($Phi)]
1808    set B [expr pow((1.0-$epsilon*$sinPhi)/(1.0+$epsilon*$sinPhi), \
1809		    0.5*$n*$epsilon)]
1810    # set phi [expr 2.0 * atan($k*$A*$B) - $pi/2.0]
1811    set phi [expr 2.0 * atan($k*$A*$B) - 1.570796326794897]
1812    # set lambda [expr $n*($Lambda-$Lambda_0)]
1813    set lambda [expr $n*($Lambda-0.33246029532469185632)]
1814
1815    # Spherical phi, lambda to auxiliary phi, lambda
1816    ## EOV starting point
1817    # set phi_0 [expr (47.0 + 06.0/60.0 + 00.0000/3600.0) / 360.0 * 2.0*$pi]
1818    set phi_0 0.8220500776893292307101144
1819    # set phi_a [expr asin(sin($phi)*cos($phi_0)
1820    #                       - cos($phi)*sin($phi_0)*cos($lambda))]
1821    set cosphi [expr cos($phi)]
1822    set phi_a [expr asin(sin($phi)*cos($phi_0)- \
1823			     $cosphi*sin($phi_0)*cos($lambda))]
1824    # set lambda_a [expr asin((cos($phi)*sin($lambda))/(cos($phi_a)))]
1825    set lambda_a [expr asin(($cosphi*sin($lambda))/cos($phi_a))]
1826
1827    # Auxiliary phi, lambda to cartesian x, y (x == northing  &&  y == easting)
1828    ## Radius of the "new" Gauss-sphere
1829    # set R 6379743.001
1830    ## Reduction (scale factor)
1831    # set m_0 0.99993
1832    set Rm_0 6379296.41898993
1833    # set x [expr $R * $m_0 * log(tan($pi/4+$phi_a/2))]
1834    # set y [expr $R * $m_0 * $lambda_a]
1835    set x [expr $Rm_0 * log(tan(0.78539816339744830961+0.5*$phi_a))]
1836    set y [expr $Rm_0 * $lambda_a]
1837
1838    # MF change: the false easting and false northing are added by the
1839    #  grid procedures; return the "raw" easting and northing
1840    return [list $y $x]
1841}
1842
1843proc ProjEOVInvert {data x y} {
1844    # return list with latitude and longitude for point at $x,$y in the terrain
1845    #  where $x is the easting and $y the northing (against the EOV convention)
1846    #  with no false easting/northing added
1847    #  $data is not used
1848    # return 0, 180 if computation of latitude diverges
1849    global MESS TXT
1850
1851    # Return to the EOV convention ($y == easting  &&  $x == northing)
1852    set z $x ; set x $y ; set y $z
1853
1854    # MF changes after this point: for efficiency, original code commented out
1855    # Cartesian y, x to auxiliary phi, lambda
1856    # set R 6379743.001
1857    ## Reduction (scale factor)
1858    # set m_0 0.99993
1859    set Rm_0 6379296.41898993
1860    ## Pi :-P
1861    # set pi 3.1415926535897932384626434
1862    # set phi_a [expr 2*atan(exp($x/($R*$m_0)))-$pi/2]
1863    set phi_a [expr 2*atan(exp($x/$Rm_0))-1.57079632679489661923]
1864    # set lambda_a [expr $y/($R*$m_0)]
1865    set lambda_a [expr $y/($Rm_0)]
1866
1867    # Auxiliary phi, lambda to spherical phi, lambda
1868    ## EOV starting point
1869    # set phi_0 [expr (47.0 + 06.0/60.0 + 00.0000/3600.0) / 360.0 * 2.0*$pi]
1870    # set phi_0 0.82205007768932923029
1871    set cosphi_a [expr cos($phi_a)]
1872    # set phi [expr asin( sin($phi_a)*cos($phi_0) +
1873    #                     cos($phi_a)*sin($phi_0)*cos($lambda_a) )]
1874    set phi [expr asin( sin($phi_a)*0.680720868959 + \
1875			    $cosphi_a*0.732542898787*cos($lambda_a) )]
1876    # set lambda [expr asin( (cos($phi_a)*sin($lambda_a)) / cos($phi) )]
1877    set lambda [expr asin( ($cosphi_a*sin($lambda_a)) / cos($phi) )]
1878
1879    # Spherical phi, lambda to HD72 Phi, Lambda
1880    ## GRS67 semi-major / semi-minor axis
1881    # set a 6378160.0
1882    # set b 6356774.516
1883    ## EOV parameters
1884    set k 1.003110007693
1885    set n 1.000719704936
1886    # set Lambda_0 [expr (19.0 + 02.0/60.0 + 54.8584/3600.0) / 360.0 * 2.0*$pi]
1887    # set Lambda_0 0.33246029532469185632
1888    ## GRS67 eccentricity
1889    # set epsilon [expr sqrt(($a*$a-$b*$b)/($a*$a))]
1890    set epsilon 0.0818205680555
1891    ## Iteration to calculate Phi
1892    set Phi $phi
1893    set i 0
1894    # set A [expr tan($pi/4+$phi/2)]
1895    set A [expr tan(0.785398163397448+0.5*$phi)]
1896    set invn [expr 1/$n]
1897    while 1 {
1898        set B [expr $k *pow( (1-$epsilon*sin($Phi)) / (1+$epsilon*sin($Phi)), \
1899				 $n*$epsilon/2)]
1900        # set Phi_new [expr 2 * atan(pow($A/$B, 1/$n)) - $pi/2]
1901	set Phi_new [expr 2 * atan(pow(1.0*$A/$B, $invn)) - 1.570796326794897]
1902	if { abs($Phi-$Phi_new) <= 1.5E-9 } {
1903            break
1904	}
1905	if { [incr i] > 20} {
1906            GMMessage [format $MESS(badinvproj) $TXT(PRJEOV)]
1907            return [list 0 180]
1908	}
1909	set Phi $Phi_new
1910    }
1911    ## Calculate Lambda
1912    # set Lambda [expr $Lambda_0 + $lambda/$n]
1913    set Lambda [expr 0.33246029532469185632 + $lambda/$n]
1914
1915    # MF change: results will be for the HD72
1916    # Convert the result from radian to DD.DDDDDDDD...
1917    set la [expr $Phi*57.29577951308232087684]
1918    set lo [expr $Lambda*57.29577951308232087684]
1919    #---
1920    return [list $la $lo]
1921}
1922
1923
1924
1925
1926