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)