ROUTINES and PROCEDURES for PERFORMING
TRANSVERSE MERCATOR PROJECTIONS
james@gentlesSPAM.info
February 1992
July 1997
July 2000
February 2003
1. Introduction
2. Basics
3. Pseudo Code for Main Routines
4. HP48 Routines: Descriptions
5. HP48 Routines: Programs (ASCII DIR object)
6. APPENDIX and Revision History
1. Introduction
This file contains a collection of HP48 programs used for translation
between Latitude/Longitude and Transverse Mercator Projection. The programs
were written specifically to translate between the Universal Transverse
Mercator standard, and specifically for the Great Britain and Ireland
"National Grid" systems, however they could be used for translation between
latitude/longitude and ANY mapping system that uses this projection by
altering the constants used to define the projection conveniently contained
in a list. In July 2000, a number of other routines were added that can be
used to calculate, Convergence, Scale Factor, and (t-T) Correction.
These programs were written for Amateur Radio use in the UK. Some amateurs
use positioning reference systems based on Latitude and Longitude, where
others use a system based on the UK National Grid which is a Transverse
Mercator Projection. This file only contains the projection translation part
of this work and as such may be useful to a wider audience.
2. Basics
"Mercator Projection" is one of the projections commonly used to show the
whole world in atlases. It's main disadvantage for 'whole world' maps is the
distortion near the poles, making Antarctica and Greenland appear much larger
than they actually are. For small parts of the globe (e.g. the UK) or slices
of the globe based around a latitude meridian this projection provides a very
accurate model, which is used in many countries.
This system allows a map to contain a regular grid for reference and short
distance calculations. However, because the map is a projection of the surface
of a shape approximating to an oblate spheroid onto a plane (or cylindrical
shaped) piece of paper then these grid lines have a very complex relationship
to latitude and longitude. In "Mercator" projection the map of the world is
projected onto a cylinder of paper from the globe, where the only point of
contact between the two is the equator. Transverse Mercator uses the same
principle, except a line of latitude is used. Absolute accuracy is only
achieved where the globe touches the map, however by choosing scalings and
sphere models correctly these errors can be minimised.
2.1 Universal Transverse Mercator (UTM)
UTM is in effect a set of 60 projections, called Zones, covering strips of
the globe from pole to pole, each 6 degrees wide. Zone 1 is centered around
longitude 177 degrees west, through to Zone 60 at 177 degrees east. (For
practical purposes, the projection is only used between -80 & +84 deg lat)
Each Zone is further subdivided into a number of Zone designators by 8degrees
of latitude. This goes from C (80 degrees south) through X (84 degrees north)
(omitting I and O, with X being 12 degrees). The Zone Designator is only
really of minor importance, since the projection itself only needs to know
whether you are in the northern or southern hemisphere.In this implementation
the convention of positive Zone for the northern hemisphere and negative
Zone for the southern hemisphere is adopted.
(The Zone designator is returned as a TAG to the Zone number when translating
from Lat/Long to UTM, but this is not required when translating backwards)
The constants required for the translation are all included in a list, this
is based on list WORLD, but values are changed depending on which of the
60 areas are being used, and whether you are in the northern or southern
hemisphere. This list is then passed to ->TMP or TMP-> to perform the
actual translation. This structure, separating the raw projection from
the specifics of various systems, allows the easy addition of other
projections to the programs.
2.2 British Isles National Grids (NGR, IGR, UK)
Although the National Grid Systems are peculiar to the British Isles, they
employ the "Transverse Mercator Projection" system as used elsewhere. The
constants required for the translation are all included in a list, called
AIRY or IRISH, and do not need such complex pre-processing as does the UTM
system. They do however require more post-processing. In the UK systems
100Km squares are given grid letters, rather than numbers. In addition the
accuracy is normally truncated to 100m for recreational use.
The false origin of the UK grid is off the SW corner of England, and there
are 4 main squares, each with 500Km sides. The square nearest the origin is
called S and covers most of England. The squares north of that are N and H.
Finally the eastern extremity of England just pokes into square T to the east
of S. These letters are chosen as part of a large 5x5 grid of 500Km squares,
similar to the subdivision of these 500Km squares themselves:
1500Km --------------
| |
G | H . | J
| / | (Northern Isles)
| . |
1000Km ---------/----
| ./---- |
M | | / N /. | O
(Scotland) | / /| ./_ |
| //| | |
500Km ------|----\--------------
| ._| \| |
R | |_ S |_ T |
(England | __- | | |
& Wales) | _-__________/ |
0 --------------------------
0 500Km 1000Km
Each of these 500Km squares is divided into 25 100Km squares as follows:
A B C D E
F G H J K
L M N O P
Q R S T U
V W X Y Z
Note that this is A to Z from top left to bottom right with letter "I"
omitted. So for example I live at "NT212752" to the nearest 100metres.
The full reference from datum for this is Eastings 321200, Northings
675200. The program XX6-> is used to do the translation from "NT212752"
to the eastings and northings in metres. ->XX6 does the inverse.
The Irish grid only covers one of the above described 500Km squares, so
it is identical to the British system except the leading N, S, H or T is
not required.
2.3 Implementation
This file contains programs written for a HP48SX calculator, but they have
been tested on a HP48GX. The programs were originally developed on a HP28S,
and in translating the code they have been improved in accuracy, with
reduced execution time, smaller memory occupancy, and taken opportunity of
the HP48s advanced features, e.g. TAGing. This work included the removal of
one error which limited accuracy in some cases. Accuracy now appears to be
better than +-3 counts of the HP48's least significant digit.
For those more interested in the mathematics of the projection, I have
included pseudo-code descriptions of the routines to aid understanding.
2.4 Credits
The routines are all based on Reference Transverse Mercator Projection
Constants, Formulae and Methods, Ordnance Survey Information March 1983
Thanks also to Peter Dana for his web pages describing UTM projection
(http://www/utexas.edu/depts/grg/gcraft/notes/coordsys/coordsys.html), and
Brian Stefanuik who needed to use UTM on the HP48 in anger in the Canadian
Northwest Territory.
3. Pseudo Code for Main Calculations
If you're only interested in using these routines, please skip this section,
and go directly to section 4 and 5.
This pseudo-code is based on the following calculator routines and on the
Ordnance and Survey Document, Transverse Mercator Projection, March 1983.
The routines use mostly two letter variable names, this does not aid
readability however these are the calculator variables. Leaving the
variables like this means that no errors have been introduced in the
translation to pseudo-code from HP48 code. Pseudo code is only provided for
the main translations, and has not been provided for distance and bearings.
Latitude/Longitude to Grid Reference
Note: All calculations performed in radians
k = Latitude of point
l = Longitude of point
Constants for appropriate spheroid:
UK Airy Irish Airy UTM etc...
a = Major semi-axis m 6377563.396 6377563.396 6378137
b = Minor semi-axis m 6356256.91 6356256.91 6356752.314
k0 = True Origin N in RADIANS 49 degrees 53.5 degrees 0 degrees
l0 = E IN RADIANS -2 degrees -8 degrees variable
n0 = False Origin: m S of true origin -100000 250000 0(N)or10^7(S)
e0 = m W of true origin 400000 200000 500000
f0 = Scale factor on Central meridian 0.9996012717 1.000035 0.9996
p = l - l0
k3 = k - k0 difference in latitude
k4 = k + k0 sum of latitudes
call subroutine ABNE, uses a b f0 and returns a1 b1 n1 e2
call subroutine AOM, uses n1 k3 k4 b1 and returns m
call subroutine VRH2, uses a1 e2 k and returns v r h2
J3 = 'm+n0' Formula 1
J4 = 'v/2*SIN(k)*COS(k)' Formula 2
J5 = 'v/24*SIN(k)*COS(k)^3*(5-TAN(k)^2+9*h2)' Formula 3
J6 = 'v/720*SIN(k)*COS(k)^5*(61-58*TAN(k)^2+TAN(k)^4)' Formula 3a
northings = 'J3+p^2*J4+p^4*J5+p^6*J6' COMPUTE NORTHINGS
(answer in metres from origin)
J7 = COS(k) * v Formula 4
J8 = 'v/6* COS(k)^3*(v/r-TAN(k)^2)' Formula 5
J9 = 'v/120*COS(k)^5*(5-18*TAN(k)^2+TAN(k)^4+14*h2-58*TAN(k)^2*h2)' Formula 6
eastings = 'e0+p*J7+p^3*J8+p^5*J9' COMPUTE EASTINGS
(answer in metres from origin)
Grid reference can now be easily obtained from eastings and northings.
Grid Reference to Latitude/Longitude
Note: All calculations performed in radians
Initially decode grid reference into metres from origin, to obtain
e = eastings of point (in metres from origin)
n = northings of point (in metres from origin)
Constants for appropriate spheroid:
UK Airy Irish Airy UTM etc...
a = Major semi-axis m 6377563.396 6377563.396 6378137
b = Minor semi-axis m 6356256.91 6356256.91 6356752.314
k0 = True Origin N in RADIANS 49 degrees 53.5 degrees 0 degrees
l0 = E IN RADIANS -2 degrees -8 degrees variable
n0 = False Origin: m S of true origin -100000 250000 0(N)or10^7(S)
e0 = m W of true origin 400000 200000 500000
f0 = Scale factor on Central meridian 0.9996012717 1.000035 0.9996
y1 = e - e0
call subroutine ABNE, uses a b f0 and returns a1 b1 n1 e2
call subroutine PHI, uses n n0 n1 k0 a1 b1 and returns k k3 k4
call subroutine VRH2, uses a1 e2 k and returns v r h2
J3 = TAN(k)/(2*r*v) Formula 7
J4 = 'TAN(k)/(24*r*v^3)*(5+3*TAN(k)^2+h2-9*TAN(k)^2*h2)' Formula 8
J5 = 'TAN(k)/(720*r*v^5)*(61+90*TAN(k)^2+45*TAN(k)^4)' Formula 9
Latitude = 'k-y1^2*J3+y1^4*J4-y1^6*J5' COMPUTE LATITUDE
J6 = 1/(COS(k)*v) Formula10
J7 = '1/(COS(k)*6*v^3)*(v/r+2*TAN(k)^2)' Formula11
J8 = '1/(COS(k)*120*v^5)*(5+28*TAN(k)^2+24*TAN(k)^4)' Formula12
J9 = '1/(COS(k)*5040*v^7)*(61+662*TAN(k)^2+1320*TAN(k)^4+720*TAN(k)^6)' 12a
Longitude = 'l0+y1*J6-y1^3*J7+y1^5*J8-y1^7*J9' COMPUTE LONGITUDE
Finally translate Latitude and longitude back into degrees, mins and secs
Grid Reference to Convergence ( Grid Bearing + Convergence = True Bearing)
Note: All calculations performed in radians
Initially decode grid reference into metres from origin, to obtain
e = eastings of point (in metres from origin)
n = northings of point (in metres from origin)
Constants for appropriate spheroid as per previous examples
y1 = e - e0
call subroutine ABNE, uses a b f0 and returns a1 b1 n1 e2
call subroutine PHI, uses n n0 n1 k0 a1 b1 and returns k k3 k4
call subroutine VRH2, uses a1 e2 k and returns v r h2
J3 = TAN(k)/(v) Formula 16
J4 = 'TAN(k)/(3*v^3)*(1+TAN(k)^2-h2-2*h2^2)' Formula 17
J5 = 'TAN(k)/(15*v^5)*(2+5*TAN(k)^2+3*TAN(k)^4)' Formula 18
Convergence = 'y1*J3-y1^3*J4+y1^5*J5' COMPUTE CONVERGENCE
Grid Reference to Scale Factor (The local scale factor for a given point)
Note: All calculations performed in radians
Initially decode grid reference into metres from origin, to obtain
e = eastings of point (in metres from origin)
n = northings of point (in metres from origin)
Constants for appropriate spheroid as per previous examples
y1 = e - e0
call subroutine ABNE, uses a b f0 and returns a1 b1 n1 e2
call subroutine PHI, uses n n0 n1 k0 a1 b1 and returns k k3 k4
call subroutine VRH2, uses a1 e2 k and returns v r h2
J3 = 1/(2*r*v) Formula 21
J4 = '(1+4*h2)/(24*r^2*v^2)' Formula 22
Factor = 'f0*(1+y1^2*j3+y1^4*j4)' COMPUTE SCALE FACTOR
Grid Reference to (T-t) Correction (for 2 eastings & northings Ea,Eb,Na,Nb)
Note: All calculations performed in radians
Initially decode grid reference into metres from origin, to obtain
ya = Ea-e0
yb = Eb-e0
n = (na+nb)/2
Constants for appropriate spheroid as per previous examples
call subroutine ABNE, uses a b f0 and returns a1 b1 n1 e2
call subroutine PHI, uses n n0 n1 k0 a1 b1 and returns k k3 k4
call subroutine VRH2, uses a1 e2 k and returns v r h2
J3 = 1/(6*r*v) Formula 23
Ga = (2*ya+yb)*(na-nb)*J3
Gb = (2*yb+ya)*(nb-na)*J3
SUBROUTINES used in Grid Ref / Latitude & Longitude TRANSFORMS
VRH2 (uses a1 e2 k and returns v r h2)
v = 'a1/SQRT(1-e2*SIN(k)^2)'
r = 'v*(1- e2)/(1-e2*SIN(k)^2)'
h2 = '(v/r)-1'
AOM (uses n1 k3 k4 b1 and returns m), Compute Developed Arc of Meridian
J3 = '(1+n1+5/4*n1^2+ 5/4*n1^3)*k3'
J4 = '(3*n1+3*n1^2+21/8*n1^3)*SIN(k3)*COS(k4)'
J5 = '(15/8*n1^2+15/ 8*n1^3)*SIN(2*k3)*COS(2*k4)'
J6 = '35/ 24*n1^3*SIN(3*k3)*COS(3*k4)'
m = 'b1*(J3 - J4 + J5 - J6)'
PHI (uses n n0 n1 k0 a1 b1 and returns k k3 k4) Iteratively Calc Lat k
k = '(n-n0)/a1+k0' First approximation to Latitude
REPEAT
k3 = 'k-k0'
k4 = 'k+k0'
call subroutine AOM, uses n1 k3 k4 b1 and returns m
Test= 'n-n0-m'
k = 'k+(n-n0-m)/a1'
UNTIL ABS(Test)<0.001
ABNE (uses a b f0 and returns a1 b1 n1 e2)
a1 = 'a*f0'
b1 = 'b*f0' Scale Major & Minor axes by f0
n1 = '(a1-b1)/(a1+b1)' calculate n from scaled axes
e2 = '(a1^2-b1^2)/a1^2' calculate e^2 (e2) from scaled axes
4. HP48 Utilities: Descriptions
Transverse Mercator Projection Routines, to change to Latitude and Longitude,
all latitude and longitude numbers are entered, and returned in DD.MMSSsss
format.
The routines fall into four categories:
4. Constants: Lists of constants for particular transformations.
e.g. AIRY.
3. Subroutines: Program element that performs a task the user wont normally
need to access, e.g. ->XX6 does string formatting and AOM
calculates the Arc of Meridian.
2. Calculations: Performs the actual translation, but the input and output
may be poor, e.g. you need to provide lists of constants,
to make it work.
1. Utilities: Performs some pre-processing to input data to make it easier
to call the level 2 routines. e.g. the currect constants
are put on the stack for the transformation you are making.
You should be able to write your own utilities to use the calculations. e.g.
UKD calculates the distance between two points assuming the UK projection
model. This could be changes to assume another model, or a particular UTM
slice.
U T I L I T I E S
->UTM Lat and Long to Universal Grid. Put Latitude on stack (north is +ve),
then Longitude on stack (east is +ve), returns the Zone Number in
stack level 3 (southern hemisphere is -ve, and the Zone designator
is in the TAG) full eastings and northings in stack 2 and 1 (both
in metres).
UTM-> Universal Grid to Lat and Long. The inverse of ->UTM. The Zone
Designator is not required, however southern hemisphere zones must
have negative Zone numbers.
->NGR Lat and Long to National Grid Reference (Great Britain). Put
Latitude on stack (north is +ve), then Longitude on stack
(east is +ve), returns a string of format "NT119779". The bottom
left corner of the 100m square defined.
NGR-> National Grid to Lat and Long. The inverse of ->NGR. Returns the
Lat and Long of the center of the 100m square defined.
->UK UK-> As above, but exact location to a fraction of a metre is returned.
This demonstrates how the routines can be used in different ways.
->IGR Lat and Long to Irish Grid Reference. Put Latitude on stack
(north is +ve), then Longitude on stack (east is +ve), returns
a string of format "T119779". The bottom left corner of the 100m
square defined.
IGR-> Irish Grid to Lat and Long. The inverse of ->IGR. Returns the
Lat and Long of the center of the 100m square defined.
UKA Takes near station eastings and northings, and far station eastings
and northings (Level 4,3,2,1 respectively) and computes the true
bearing including Convergence and (t-T), returning the bearing in
D.MSsss format in stack level 1.
UKD Takes 2 station eastings and northings (as per UKA) and computes the
true distance between them using Simpson's Rule on the stations and
the mid point between them. Distance returned in metres on the stack.
C A L C U L A T I O N S
->TMP Performs the Transverse Mercator Projection.On entry the stack should
contain Lat,Long (both in radians),then all the constants in the list
(e.g. IRISH) depending on what projection is being used. The routine
returns eastings and northings on the stack in metres.
TMP-> Performs the inverse Transverse Mercator Projection. On entry the
stack should contain eastings, northings, then all the constants in
the list (e.g. IRISH) depending on what projection is being used. The
routine returns latitude and longitude on the stack in radians.
EN->C Compute the Convergence for a given easting and northing in metres.
On entry the stack should contain eastings, northings, then all the
constants in the list (e.g. IRISH) depending on what projection is
being used. The routine returns Convergence on the stack in radians.
EN->F Compute the Scale Factor for a given easting and northing in metres.
On entry the stack should contain eastings, northings, then all the
constants in the list (e.g. IRISH) depending on what projection is
being used. The routine returns a number close to one on the stack.
EN->T Compute the (T-t) Correction for a given easting and northing in mtrs.
On entry the stack should contain eastings A, northings A, eastings B,
northings B,then all the constants in the list (e.g. IRISH) depending
on what projection is being used. The routine returns a (T-t)a and
(T-t)b on the stack in radians.
S U B R O U T I N E S
->XX6 Used by ->NGR and ->IGR to translate eastings and northings in metres
to a national grid reference, e.g. 311800 678485 become "NT118784".
Truncates accuracy to 100m,so result refers to the bottom left corner
of the grid reference 100metre square produced. If the coordinates do
not fall into one of the 4 500Km squares a "#" is returned in the
first character of the string. Used by routines ->NGR & ->IGR.
XX6-> Performs the inverse process of ->XX6. Returns the full eastings and
northings of the square center referred to in metres, i.e. +50m is
added to the eastings and northings. Used by routines NGR-> & IGR->.
VRH2 Used by ->TMP and TMP-> routines. Calculates v r & h2.
AOM Used by ->TMP and PHI routines. Calculates Arc of Meridian, m.
PHI Used by TMP-> routine.Iteratively calculates latitude data k k3 & k4.
ABNE Used by ->TMP and TMP-> routines. Calculates a1 b1 n1 & e2.
PAD0 Used by ->XX6 routine to pad leading zeros to a length of 3.
C O N S T A N T S
AIRY Lists of constants used by
IRISH ->TMP, TMP-> & other routines
WORLD to define the spheroid & scales.
External programs used (included in the listing for completeness):
PRESERVE Returns flags to initial state on program exit. Used in ->TMP, TMP->,
EN->C, EN->F, EN->T, & ->XX6 routines. See HP48SX Manual II page 555.
5. HP48 Routines: Programs (ASCII DIR Object) xxxxbytes #xxxxhex
----------------------------cut here------------------------------------
%%HP: T(3)A(D)F(.);
DIR
\->UTM
\<< \-> la lo @ take lat and long from stack
\<< lo 180 + 360 MOD 6 / IP 1 + @ Compute zone number from longitude
IF la 0 < THEN NEG END @ If southern hemisphere negative
la 8 / 77 + IP @ Compute Zone designator letter
IF DUP 72 > THEN 1 + END @ correct for missing letter I
IF DUP 78 > THEN 1 + END @ correct for missing letter O
66 MAX 89 MIN @ limit letters to range B - Y
IF DUP 89 == la 84 \<= AND @ artificially make X 12 degrees
THEN 1 - END
CHR "Zone " SWAP + \->TAG @ build tag for Zone
DUP ABS 6 * 183 - D\->R @ Compute central meridian
WORLD 4 3 ROLL PUTI @ put central meridian in const list
la IF 0 > THEN 0 PUTI END @ If northern hemi, put 0 in false Ns
DROP lo HMS\-> D\->R @ Put long on stack and translate
la HMS\-> D\->R @ Put lat on stack and translate
3 ROLL OBJ\-> DROP \->TMP @ Roll stack and run translation
"n" \->TAG SWAP
"e" \->TAG SWAP
\>>
\>>
@
UTM\->
\<< \-> zn e n @ take params from stack
\<< e n WORLD
4 zn ABS 6 * 183 - D\->R PUTI @ Compute central meridian,put in const
IF zn 0 > THEN 0 PUTI END @ If northern hemi, put 0 in false Ns
DROP OBJ\-> DROP TMP\-> SWAP @ drop stack and run translation
R\->D \->HMS "Lat\^o" \->TAG @ tag outputs
SWAP R\->D \->HMS "Long\^o" \->TAG
\>>
\>>
@
\->NGR @ Lat & Long to XX123456
\<< HMS\-> D\->R SWAP HMS\-> D\->R @ Change to degs decimal then rads
AIRY OBJ\-> DROP @ Put projection consts on stack
\->TMP @ Do projection, returns e&n on stack
\->XX6 @ translate to grid letters + 6 nums
"NGR" \->TAG @ tag output
\>>
@
NGR\-> @ XX123456 to Lat & Long
\<< DTAG XX6\-> @ translate grid to full coordinates
AIRY OBJ\-> DROP @ put projection consts on stack
TMP\-> @ do projection, returns Lat & Long
SWAP R\->D \->HMS "Lat\^o" \->TAG @ translate to degrees mins secs
SWAP R\->D \->HMS "Long\^o" \->TAG @ and tag output
\>>
@
\->UK @ Lat & Long to e & n
\<< HMS\-> D\->R SWAP HMS\-> D\->R @ Change to degs decimal then rads
AIRY OBJ\-> DROP @ Put projection consts on stack
\->TMP @ Do projection, returns e&n on stack
SWAP "e" \->TAG SWAP "n" \->TAG @ tag output
\>>
@
UK\-> @ e & n to Lat & Long
\<< AIRY OBJ\-> DROP @ put projection consts on stack
TMP\-> @ do projection, returns Lat & Long
SWAP R\->D \->HMS "Lat\^o" \->TAG @ translate to degrees mins secs
SWAP R\->D \->HMS "Long\^o" \->TAG @ and tag output
\>>
@
\->TMP @ Reference Transverse Mercator Proj'n
\<< @ Constants, Formulae and Methods.
\<< RAD \-> l k a b k0 l0 n0 e0 f0 @ Ordnance Survey Information March 83
\<< l l0 - k k0 - k k0 + \-> p k3 k4
\<< a b f0 ABNE \-> a1 b1 n1 e2
\<< n1 k3 k4 b1 AOM a1 e2 k VRH2
l l0 - k SIN k COS k TAN @ Compute SIN(k) COS(k) & TAN(k) once
\-> m v r h2 p s c t
\<< m n0 + @ Formula I
s c v 2 / * * p SQ * @ Formula II *p^2
'v/24*s*c^3*(5-t^2+9*h2)'
\->NUM p 4 ^ * @ Formula III *p^4
'v/720*s*c^5*(61-58*SQ(t)+t^4)'
\->NUM p 6 ^ * @ Formula IIIA*p^6
+ + + @ compute northings
c v * p * @ Formula IV *p
'v/6*c^3*(v/r-SQ(t))' \->NUM p 3 ^ * @ Formula V *p^3
'v/120*c^5*(5-18*SQ(t)+t^4+14*h2
-58* SQ(t)*h2)' \->NUM p 5 ^ * @ Formula VI *p^5
+ + e0 + @ compute eastings
SWAP
\>>
\>>
\>>
\>> @ Return to original flag settings
\>> PRESERVE @ See HP48 Manual II page 555
\>>
@
TMP\-> @ Reference Transverse Mercator Proj'n
\<< @ Constants, Formulae and Methods.
\<< RAD \-> e n a b k0 l0 n0 e0 f0 @ Ordnance Survey Information March 83
\<< a b f0 ABNE e e0 - \-> a1 b1 n1 e2 y1
\<< n n0 n1 k0 a1 b1 PHI \-> k k3 k4
\<< a1 e2 k VRH2 k COS k TAN @ Calculate COS(k)&TAN(k) for use later
\-> v r h2 c t
\<< k t 2 r v * * / y1 SQ * @ Formula VII *y1^2
't/(24*r*v^3)*(5+3*SQ(t)+h2
-9*SQ(t)*h2)' \->NUM y1 4 ^ * @ Formula VIII *y1^4
't/(720*r*v^5)*(61+90*SQ(t)
+45*t^4)' \->NUM y1 6 ^ * @ Formula IX *y1^6
- - - @ Compute Latitude
c v * INV y1 * @ Formula X *y1
'1/(c*6*v^3)*(v/r+2*SQ(t))'
\->NUM y1 3 ^ * @ Formula XI *y1^3
'1/(c *120*v^5)*(5+28*SQ(t)+24
*t^4)' \->NUM y1 5 ^ * @ Formula XII *y1^5
'1/(c*5040*v^7)*(61+662*SQ(t)+1320
*t^4+720*t^6)' \->NUM y1 7 ^ * @ Formula XIIA *y1^7
- - - l0 + @ Compute Longitude
\>>
\>>
\>>
\>> @ Return to original flag settings
\>> PRESERVE @ See HP48 Manual II page 555
\>>
@
UKA @ TRUE BEARING BETWEEN 2 POINTS
\<< @ for UK Grid and speroid
\<< RAD AIRY \-> en nn ef nf sp @ Es/Ns Near&Far + Spheriod
\<< en nn ef nf sp OBJ\-> DROP EN\->T DROP @ (T-t)a (throw away (t-T)b)
en nn sp OBJ\-> DROP EN\->C @ Convergence C
nf nn - DUP SQ ef en - SQ + \v/ / ACOS @ Plane angle calculation
IF en ef > THEN -1 ELSE 1 END * @ angles > 180 degrees negated
+ SWAP - R\->D \->HMS "UKA" \->TAG @ Combine & translate
\>>
\>> PRESERVE
\>>
@
UKD @ TRUE DISTANCE BETWEEN 2 POINTS
\<< @ for UK Grid and speroid
\<< RAD AIRY \-> en nn ef nf sp @ Es/Ns Near&Far + Spheriod
\<< en nn sp OBJ\-> DROP EN\->F @ Near scale factor
ef nf sp OBJ\-> DROP EN\->F @ Far scale factor
en ef + 2 / nn nf + 2 /
sp OBJ\-> DROP EN\->F @ mid scale factor
4 * + + 6 / @ Simpsons Rule
en ef - SQ nn nf - SQ + \v/ @ raw distance
* "UKD" \->TAG @ SCALE distance
\>>
\>> PRESERVE
\>>
@
EN\->C @ Reference Transverse Mercator Proj'n
\<< @ Constants, Formulae and Methods.
\<< RAD \-> e n a b k0 l0 n0 e0 f0 @ Ordnance Survey Information March 83
\<< a b f0 ABNE e e0 - \-> a1 b1 n1 e2 y1
\<< n n0 n1 k0 a1 b1 PHI \-> k k3 k4
\<< a1 e2 k VRH2 k TAN @ Calculate TAN(k) for use later
\-> v r h2 t
\<< t v / y1 * @ Formula XVI *y1
't/(3*v^3)*(1+SQ(t)-h2
-2*SQ(h2))' \->NUM y1 3 ^ * @ Formula XVII *y1^3
't/(15*v^5)*(2+5*SQ(t)
+3*t^4)' \->NUM y1 5 ^ * @ Formula XVIII*y1^6
- - @ Compute Convergence
\>>
\>>
\>>
\>> @ Return to original flag settings
\>> PRESERVE @ See HP48 Manual II page 555
\>>
@
EN\->F @ Reference Transverse Mercator Proj'n
\<< @ Constants, Formulae and Methods.
\<< RAD \-> e n a b k0 l0 n0 e0 f0 @ Ordnance Survey Information March 83
\<< a b f0 ABNE e e0 - \-> a1 b1 n1 e2 y1
\<< n n0 n1 k0 a1 b1 PHI \-> k k3 k4
\<< a1 e2 k VRH2 \-> v r h2
\<< 2 r v * * INV y1 SQ * @ Formula XXI *y1^2
'(1+4*h2)/(24*r^2*v^2' \->NUM y1 4 ^ * @ Formula XXII *y1^4
+ 1 + f0 * @ Compute Scale Factor
\>>
\>>
\>>
\>> @ Return to original flag settings
\>> PRESERVE @ See HP48 Manual II page 555
\>>
@
EN\->T @ Reference Transverse Mercator Proj'n
\<< @ Constants, Formulae and Methods.
\<< RAD \-> ea na eb nb a b k0 l0 n0 e0 f0 @ Ordnance Survey Information March 83
\<< a b f0 ABNE e e0 - \-> a1 b1 n1 e2 y1
\<< na nb + 2 / n0 n1 k0 a1 b1 PHI \-> k k3 k4
\<< a1 e2 k VRH2
ea e0 - eb e0 - \-> v r h2 ya yb
\<< 6 r v * * INV DUP @ Formula XXIII
'(2*ya+yb)*(na-nb)' \->NUM * SWAP @ Ga
'(2*yb+ya)*(nb-na)' \->NUM * @ Gb
\>>
\>>
\>>
\>> @ Return to original flag settings
\>> PRESERVE @ See HP48 Manual II page 555
\>>
@
\->IGR @ Lat & Long to X123456
\<< HMS\-> D\->R SWAP HMS\-> D\->R @ Same as \->NGR but since Irish grid
IRISH OBJ\-> DROP @ is only one 500Km square then the
\->TMP \->XX6 @ initial grid letter "S" is dropped.
2 8 SUB "IGR" \->TAG
\>>
@
IGR\-> @ X123456 to Lat & Long
\<< DTAG "S" SWAP + @ Same as NGR\-> but since Irish grid
XX6\-> IRISH OBJ\-> DROP @ is only one 500Km square then an
TMP\-> @ extra "S" is added so the standard
SWAP R\->D \->HMS "Lat\^o" \->TAG @ routines can be used.
SWAP R\->D \->HMS "Long\^o" \->TAG
\>>
@
\->XX6
\<<
\<< STD 100 / IP SWAP 100 / IP @ throw away accuracy better than 100m
DUP2 5000 / IP SWAP 5000 / IP @ break eastings and northings into
\-> n e es ns @ 100 metre units and 500Km units
\<<
IF es 0 ==
THEN @ for eastings 0 to 500000
CASE ns 2 == @ for northings 1000K to 1500K square H
THEN "H"
END
ns 1 == @ for northings 500K to 1000K square N
THEN "N"
END
ns 0 == @ for northings 0 to 500K square S
THEN "S"
END
"#" @ else illegal square denoted by a #
END
ELSE
ns 0 == es 1 == AND @ for eastings 500K to 1000K AND
"T" @ northings 0 to 500K square T
"#" @ else illegal square denoted by a #
IFTE
END
e 1000 / IP 5 MOD 1 + @ Compute a number for each 100Km square
4 n 1000 / IP 5 MOD - @ within the 500Km square in the range
5 * + DUP @ 1 to 25 top left to bottom right
IF 8 \<= @ skip the number 8, gives 0-7,9-25
THEN 1 -
END 65 + CHR + @ translate to A to Z with I missing
e PAD0 n PAD0 @ use PAD0 to pad e & n with leading 0's
+ +
\>>
\>> PRESERVE @ Return to original flag settings
\>> @ See HP48 Manual II page 555
@
XX6\->
\<< DTAG DUP 3 5 SUB @ Split incoming string into 4 parts
SWAP DUP 6 8 SUB @ and put them on stack
SWAP DUP NUM @ e.g. "NJ123456"
SWAP 2 2 SUB NUM 66 - @ "123" "456" NUM(N) NUM(J)-66
DUP @ e n f s
IF 7 < @ NUM(J)-66 are translated to be
THEN 1 + @ contiguous top left to bottom right
END
CASE OVER 72 == @ H square msdigits
THEN 0 2
END
OVER 78 == @ N square msdigits
THEN 0 1
END
OVER 83 == @ S square msdigits
THEN 0 0
END
OVER 84 == @ T square msdigits
THEN 1 0
END
"XX Error" DOERR @ else display error
END
3 PICK 5 MOD 4 5 PICK 5 / IP 5 MOD -
\-> e n f s e1 n1 e2 n2
\<< e OBJ\-> e1 5000 *
e2 1000 * + + 100 * 50 +
n OBJ\-> n1 5000 * @ recombine elements into full N&E
n2 1000 * + + 100 * 50 + @ adding 50 metres to return square ctr
\>>
\>>
@
VRH2 @ Compute v r and h2
\<< \-> a1 e2 k
\<< 'a1/\v/(1-e2* SIN(k)^2)' \->NUM
DUP
'(1-e2)/(1-e2*SIN(k)^2)' \->NUM *
DUP2 / 1 -
\>>
\>>
@
AOM @ Compute Developed Meridian Arc m
\<< 4 PICK DUP SQ SWAP 3 ^ @ Compute n1^2 and n1^3 only once
\-> n1 k3 k4 b1 n2 n3
\<< '(1+n1+5/4*n2+5/4*n3)*k3' \->NUM
'(3*n1+3*n2+21/8*n3)*SIN(k3)*COS(k4)' \->NUM
'(15/8*n2+15/8*n3)*SIN(2*k3)*COS(2*k4)' \->NUM
'35/24*n3*SIN(3*k3)*COS(3*k4)' \->NUM
- - - b1 * @ on exit m is on stack
\>>
\>>
@
PHI @ iteratively calculate Latitude k
\<< \-> n n0 n1 k0 a1 b1 @ stack to local variables
\<< '(n-n0)/a1+k0' \->NUM @ first approximation to Latitude
0 0
DO
DROP2 DUP DUP k0 - SWAP k0 +
\-> k k3 k4
\<< n1 k3 k4 b1 AOM \-> m @ re-compute arc of meridian
\<< n n0 - m - ABS DUP
IF .001 <
THEN k
ELSE n n0 - m - a1 / k +
END
k3 k4
\>>
\>> 4 ROLL
UNTIL .001 < @ until error is less than 0.001
END
\>> @ returns k k3 k4 on stack
\>>
@
ABNE
\<< DUP ROT * 3 ROLLD * \-> b1 a1 @ scale a & b by f0 and call a1 and b1
\<< a1 b1 DUP2 - a1 b1 + / @ calculate n from scaled axes
a1 SQ b1 SQ - a1 SQ / @ calculate e^2 (e2) from scaled axes
\>>
\>>
@
PAD0
\<< 1000 MOD \->STR @ take string from stack and
WHILE @ while it is shorter than 3 characters
DUP SIZE 3 <
REPEAT "0" SWAP + @ add a leading "0"
END
\>>
@
@ Constants defining Transverse Mercator
AIRY { @ Projection for Airy Spheroid
6377563.396 @ Major semi-axis m
6356256.91 @ Minor semi-axis m
.855211333477 @ True Origin 49 deg N for
-3.490658E-2 @ 2 deg W GREAT BRITAIN
-100000 @ False Origin in m S of true origin
400000 @ m W of true origin
.9996012717 } @ Scale factor on Central meridian Fo
@
@ Constants defining Transverse Mercator
IRISH { @ Projection for Airy Spheroid
6377563.396 @ Major semi-axis m
6356256.91 @ Minor semi-axis m
.933751149817 @ True Origin 53.5 deg N for
-.139626340159 @ 8 deg W IRELAND
250000 @ False Origin in m S of true origin
200000 @ m W of true origin
1.000035 } @ Scale factor on Central meridian Fo
@
@ Constants defining Transverse Mercator
WORLD { @ Projection for Universal TMP
6378137 @ Major semi-axis m
6356752.314 @ Minor semi-axis m
0 @ True Origin 0 rads Latitude
0 @ True Origin Long (varies for each Zone)
10000000 @ False Origin m S (zero for north hemi)
500000 @ m W of true origin
.9996 } @ Scale factor on Central meridian Fo
@
PRESERVE @ See HP48SX Manual II page 555.
\<< RCLF \-> f @ May be more suitably located in HOME
\<< EVAL f STOF
\>>
\>>
END
APPENDIX
OTHER SPHEROID DATA
Spheroid: | Semi-major axis |Eccentricity sqrd(e),| Commonly used for:
|(Equatorial Radius)|Flattening (f), |
| (a): |or Polar Radius (b): |
_____________|___________________|_____________________|____________________
airy | a=6377563.396 | e=.006670540 or |Assumes scale factor
| | b=6356256.91 |
australian | a=6378160 | f=1/298.25 | Australia
bessel | a=6377397.155 | e=.006674372 | Japan
clark66 | a=6378206.4 | b=6356583.8 | N. America
everest | a=6377276.345 | e=.0066378466 | India, Burma
grs80 | a=6378137 | f=1/298.257 |
hayford | a=6378388 | f=1/297 |
international| a=6378388 | f=1/297 | Europe
krasovsky | a=6378245 | f=1/298.3 |
wgs66 | a=6378145 | f=1/298.25 | worldwide coverage
wgs72 | a=6378135 | f=1/298.26 | worldwide coverage
wgs84 | a=6378137 | f=1/298.257223563 | worldwide coverage
utm?? | a=6378137 | b=6356752.314 |Scale Factor: .99960
REVISION HISTORY
1989-1990: Original Document for the HP28 Calculator
February 1992: Originally published for the HP48.
Modified from the original HP28 version written around 1990.
July 1997: Various alteraition, key ones to the numerical integration
which improved accuracy.
July 2000: Minor changes to ensure compatability with HP48GX
February 2003: Error in Pseudo Code PHI routine corrected
(this error was not in the HP48 routines themselves)