UTM <-> Lat/Lon for WP-34S
Message #1 Posted by Tom Grydeland on 5 Sept 2013, 1:31 p.m.
Hi all,
Since I saw WGS-84 constants in the catalogue, I thought this handy little machine should be able to do the conversion between Lat/Lon and UTM. After a short detour in Snyder's series, I realised that The Wikipedia entry had a development that was both terser and more accurate (I verified with a small Python implementation that I can transform to UTM and back with an introduced inaccuracy in the micrometer range even when the transformation was forced to use UTM zone 5 away from its "true" zone, once I had gone to n^4 in the expansions.)
I have a working implementation on my WP-34s, and I believe the listing below is the same, modulo local label numbers. In a fully loaded calc, I must reduce to 80 global registers to run the code, but I'm sure it can be shortened. The '34 version goes only to order n^3, which means that the inaccuracy is around 0.2 mm instead of micrometers. Good enough for my purposes.
Issues:
- Does not handle southern hemisphere (users must add/subtract 1e7 meters themselves)
- If the first conversion is UTM->Lat/Lon, zone must be prestored [C], but the code does not check for it.
- No prompting or instructions
Usage:
- XEQ 'UTM' // Runs intialization and stops.
At this point, the labels A, B and C are active as follows
- [A] transforms latitude (Y) and longitude (X), both in decimal degrees, to Northing (Y) and Easting (X) in metres. Rest of stack is preserved. If zone was set manually (with [C]), conversion takes place in the given zone, otherwise zone is given by coordinates.
- [B] transforms Northing (Y) and Easting (X) in metres to latitude (Y) and longitude (X), both in decimal degrees. Rest of stack is preserved. The zone must have been set (either with [C] or with a transformation to UTM [A]) before this function can be used, but the program does not check.
- [C] is used to recall UTM zone, or to store/force a zone (when ENTRY? is true). Store 0 to return to automatically determined zone.
Any number of transformations can be performed while the program is active. Global registers 10-23 should not be touched between transformations, and the remainder of the registers 00-36 will be modified during transformations.
Flags 00-03 are also used by the program.
Suggestions are welcome, as always.
What an amazing little machine to have in one's pocket! Thanks to its creators!
--T
Edit: had to fix the last bug before reposting the code
/*
* Allocation of registers
* r00 - r02 scratch vector T1
* r03 - r05 scratch vector T2
* r06 - r08 scratch vector T3
* r10 - r12 alpha_i
* r13 - r15 beta_i
* r16 - r18 delta_i
* r19 - UTM zone number
* r20 - lambda_0, UTM zone center longitude, in radians
* r21 - k0 (0.9996)
* r22 - A = a/(1+n) * (1 + n**2/4 + n**4/64)
* r23 - n (0.00167922039)
*/
/*
* Labels:
* 'UTM' : main entry point, initializes fixed vectors (alpha, beta, delta)
* and zone if the X value is an integer 1 <= X <= 60
*
* A : lat lon -> utm_n utm_e
* B : utm_n utm_e -> lat lon
* C : sto/rcl UTM zone
*
* dot2 : two-vector product-sum (dot product) VDy VDx -> X = Y [dot] Z
* dot3 : three-vector product-sum VDz VDy VDx -> (X [dot] Y [dot] Z)
* v_trig : create vector of trig(i*x) terms; xi Vi -> xi
* where 'trig' is one of [cos sin cosh sinh] for flags 1/2 of [0/0 1/0 0/1 1/1]
* X points at first register to modify before call,
* and returns the angle in X (in case next vector uses same)
*/
/*
* Flags:
* 00: set when UTM zone is chosen manually, cleared (default) otherwise
* 01: set when v_trig computes sin(h), cleared for cos(h)
* 02: set when v_trig uses hyperbolic versions
* 03: set when dot3 is used for two vectors only, cleared otherwise
*
* During init
* 01: set when power series in n has order-1 term
* 02: set when power series in n has order-2 term
*/
LBL'UTM'
// Initialization:
RAD
LocR 10
STOS .00
// TODO: Initialize UTM zone
// r21: k0 - scale factor at central latitude
.
9
9
9
6
STO 21
// r23: n - third flattening
# Sf[^-1]
2
*
DEC X
1/x
STO 23
// r22: A - length of elliptic arc from equator to pole
// n is on stack already
STO Z
INC Z
x^2
STO Y
# 64
1/x
*
# 4
1/x
+
*
INC X
RCL/ Y
# Sa
*
STO 22
CF 00 // no UTM zone set
SF 01 // n-term present
SF 02 // n**2-term present
CF 03
// alpha_1
# 1/2
# 2
# 3
/
+/-
STOS 00 // will be used again in beta_1 and delta_1
# 5
# 16
/
XEQ npow
STO 10
// beta_1
RCLS 00
# 37
# 96
/
XEQ npow
STO 13
//delta_1
RCLS 00
# 2
STO Z
+/-
XEQ npow
STO 16
CF 01 // n-term absent
// alpha_2
# 13
# 48
/
# 3
# 5
/
+/-
XEQ npow
STO 11
// beta_2
# 48
1/x
# 15
1/x
XEQ npow
STO 14
// delta_2
# 7
# 3
/
# 8
# 5
/
+/-
XEQ npow
STO 17
CF 02 // n**2-term absent
// alpha_3
# 61
# 240
/
XEQ npow
STO 12
// beta_3
# 17
# 48
/
# 10
/
XEQ npow
STO 15
// delta_3
# 56
# 15
/
XEQ npow
STO 18
RCLS .00
RTN // initialization done, wait for orders
/* Convert lat/lon to UTM E/N and zone
*
* lat/lon to UTM register allocations
* r30 - phi (latitude in radians)
* r31 - dla (delta lambda in radians)
* r32 - t (intermediate value)
* r33 - xi' (intermediate angle)
* r34 - eta' (intermediate angle)
* r.00-.07 saved stack
* r.00 - Easting (overwrites longitude)
* r.01 - Northing (overwrites latitude)
*
* On entry:
* latitude in degrees in Y
* longitude in degrees in X
*
*/
LBL A
LocR 10
STOS .00
FC? 00
XEQ zone
// At this point, zone is in R19 and lambda0 in r20
DEG[->] // X: longitude in radians
RCL- 20
STO 31
x<> Y
DEG[->]
STO 30
// Compute t
SIN
ATANH
RCL L
RCL 23
[sqrt]
RCL L
INC X
/
# 2
*
<> XYXZ
*
ATANH
*
-
SINH
STO 32
// Compute xi'
RCL 31
COS
/
ATAN
STO 33
// Compute cos(2*j*xi') -> T1 [0-2]
# 2
*
CF 01
CF 02
# 0
XEQ v_trig
// Compute sin(2*j*xi') -> T2 [3-5]
SF 01 // sin
# 3
XEQ v_trig
// Compute eta'
RCL 31
SIN
RCL 32
x^2
INC X
[sqrt]
/
ATANH
STO 34
// Compute sinh(2*j*eta') -> T3 [6-8]
# 2
*
SF 02
# 6
XEQ v_trig
// Compute easting
# 10 // alpha
# 0 // T1 = cos(2*j*xi')
# 6 // T3 = sinh(2*j*eta')
XEQ dot3 // sum(...)
RCL+ 34
RCL* 21
RCL* 22 // A*k0*(eta + sum(...))
5
EEX
5
+
STO .00
// Compute cosh(2*j*eta') -> T1 [0-2]
# 2
RCL* 34
CF 01 // cosh
# 0 // T1
XEQ v_trig
// Compute northing
# 10 // alpha
# 3 // T2 = sin(2*j*xi')
# 0 // T1 = cosh(2*j*eta')
XEQ dot3 // sum(...)
RCL+ 33
RCL* 21
RCL* 22 // A*k0*(xi + sum(...))
// TODO: Add false northing in southern hemisphere
STO .01 // place northing in Y in returned stack
// Restore stack and return!
RCLS .00
RTN
/* Convert UTM E/N and zone to lat/lon
*
* UTM to lat/lon register allocations
* r30 - phi (latitude in radians)
* r31 - dla (delta lambda in radians)
* r32 - chi (intermediate value)
* r33 - xi' (intermediate angle)
* r34 - eta' (intermediate angle)
* r35 - xi (rescaled latitude)
* r36 - eta (rescaled longitude)
* r.00-.07 saved stack
*
* On entry:
* northing in metres in Y
* easting in metres in X
*
*/
LBL B
LocR 10
STOS .00
// FC? 00
// prompt for zone
5
EEX
5
-
RCL/ 21
RCL/ 22
STO 36 // eta = (E-E0)/(k0*A)
x<> Y
// TODO: subtract false northing in southern hemisphere
RCL/ 21
RCL/ 22
STO 35 // xi = (N-N0)/(k0*A)
// Compute cos(2*j*xi) -> T1 [0-2]
# 2
*
CF 01
CF 02
# 0
XEQ v_trig
// Compute sin(2*j*xi) -> T2 [3-5]
SF 01
# 3
XEQ v_trig
// Compute sinh(2*j*eta) -> T3 [6-8]
# 2
RCL* 36
SF 02 // sinh
# 6
XEQ v_trig
// compute eta'
# 13 // beta
# 0 // T1 = cos(2*j*xi)
# 6 // T3 = sinh(2*j*eta)
XEQ dot3
RCL- 36
+/-
STO 34 // eta' = eta - sum(...)
// Compute cosh(2*j*eta) -> T1 [0-2]
# 2
RCL* 36
CF 01 // cosh
# 0 // T1
XEQ v_trig
// Compute xi'
# 13 // beta
# 3 // T2 = sin(2*j*xi)
# 0 // T1 = cosh(2*j*eta)
XEQ dot3
RCL- 35
+/-
STO 33 // xi' = xi - sum(...)
// Compute chi
SIN
RCL 34
COSH
/
ASIN
STO 32
// Compute phi, which is latitude in radians
# 2
*
SF 01
CF 02 // sin
# 0
XEQ v_trig
# 16 // delta
# 0 // T1 = sin(2*j*chi)
XEQ dot2
RCL+ 32 // phi = chi + sum(...)
rad[->][degree] // latitude
STO .01
// delta lambda, then longitude
RCL 34
SINH
RCL 33
COS
/
ATAN // delta lambda
RCL+ 20 // Translate to correct UTM zone
rad[->][degree] // longitude
STO .00
RCLS .00
RTN
LBL C
ENTRY? // Setting zone
JMP set_zone
RCL 19
RTN
set_zone:: STO 19
x[!=]0?
SKIP 3
STO 20
CF 00
RTN
# 6
*
# 183
-
DEG[->]
STO 20
DROP
RCL 19
SF 00
RTN
// Determine correct UTM zone
// Y: lat (in degrees)
// X: lon (in degrees)
// Stack is left undisturbed, zone in r19, lambda_0 in r20
zone:: LocR 10 // STOS bug, could be 8 (or 4)
STOS .00 // lon in .00, lat in .01
# 6
/
FLOOR
# 31
+
STO 19
RCL L // look for SW Norway anomaly
x[!=]? 19
JMP zone_next
# 56
x>? .01 // lat < 56?
JMP zone_next
# 64
x<? .01 // lat > 64?
JMP zone_next
INC 19
JMP zone_done
zone_next:: # 72 // look for Svalbard anomaly
x>? .01 // lat < 72?
JMP zone_done
# 6
x>? .00 // lon < 6?
JMP zone_done
# 36
x<? .00 // lon > 36?
JMP zone_done
RCL .00
# 3
+
# 12
/
FLOOR
# 2
*
# 31
+
STO 19
zone_done:: RCL 19
# 6
*
# 183
-
DEG[->]
STO 20
RCLS .00
RTN
// Z: label to execute for xi -> trig(xi)
// Y: xi
// X: first register (of 3)
v_trig:: STO I // start of vector in I, theta in J
x<> Y
STO J
XEQ trig
STO[->]I
INC I
# 2
RCL* J
XEQ trig
STO[->]I
INC I
# 3
RCL* J
XEQ trig
STO[->]I
RCL J
RTN
trig:: FS? 01
JMP trig_sin
FS? 02
COSH
FC? 02
COS
RTN
trig_sin:: FS? 02
SINH
FC? 02
SIN
RTN
// takes three power series coefficients from the stack, n^3 in x, n^2 in y n in z
// used to compute alpha_i, beta_i, delta_i
npow:: RCL* 23
FS? 02
+
RCL* 23
FS? 01
+
RCL* 23
RTN
// Takes vector start numbers in X, and Y, drops them both
// off the stack and returns their dot product in X.
// Uses only local registers.
dot2:: SF 03
// Takes vector start numbers in X, Y, and Z, drops them all
// off the stack and returns their triple-dot product in X.
// Uses only local registers.
dot3:: LocR 4
STO .00
DROP
STO .01
DROP
FS? 03
SKIP 2
STO .02
DROP
# 3
STO .03
DROP
# 0
dot3_top:: RCL[->].00
RCL*[->].01
FC? 03
RCL*[->].02
+
INC .00
INC .01
FC? 03
INC .02
DSE .03
JMP dot3_top
CF 03
END
Edited: 5 Sept 2013, 5:27 p.m.
| 
Go back to the main exhibit hall