Mercator Pri Ject I e
Mercator Pri Ject I e
Mercator Pri Ject I e
MERCATOR
PROJECTIONS
M ERCATOR PROJECTIONS ON
PETER OSBORNE
EDINBURGH
2008
The Transverse Mercator projection is the basis of many maps covering individual countries, such as Australia and Great Britain, as well
as the set of American UTM projections covering the whole world
(other than the polar regions). Such maps are invariably covered by a
set of grid lines. It is important to appreciate the following two facts
about the Transverse Mercator projection and the grids covering it:
1. Only one grid line runs true northsouth. Thus in Britain only
the grid line coincident with the meridian at 2 W is true: all
others are slightly distorted. The UTM series is a set of 60
projections covering a width of of 6 in latitude: the grid lines
run true northsouth only on the central meridians at 3 E, 9 E,
15 E, . . .
2. The scale on the maps derived from Transverse Mercator projections is not uniform: it is a function of position. For example the Landranger maps of the Ordnance Survey of Great
Britain have a nominal scale of 1:50000: this value is only exact on two slightly curved lines almost parallel to the central
meridian at 2 W.
The above facts are unknown to the majority of map users. They are
the subject of this article together with the presentation of formulae
relating latitude and longitude to grid coordinates.
Preface
For many years I had been intrigued by the the statement on the (British) Ordnance
Survey maps pointing out that the grid lines are not exactly aligned with meridians and parallels: four precise figures give the magnitude of the deviation at each corner of the map
sheets. My first retirement project has been to find out exactly how these figures are calculated and this led to an exploration of all aspects of the Transverse Mercator projection on
an ellipsoid of revolution (TME). This projection is also used for the Universal Transverse
Mercator series of maps covering the whole of the Earth, except for the polar regions.
The formulae for TME are given in many books and web pages but the full derivations
are only to be found in original publications which are not readily accessible: therefore I decided to write a short article explaining the derivation of the formulae. Pedagogical reasons
soon made it apparent that it would be necessary to start with the normal and transverse
Mercator projection on the sphere before going on to discuss the normal and transverse
Mercator projection on the ellipsoid. As a result the length of this document has doubled
and redoubled, but I have resisted the temptation to cut out the details which would be
straightforward for a professional but daunting for a layman. The mathematics involved
is not difficult (depending on your point of view) but it does require the rudiments of complex analysis for the crucial steps. On the other hand the algebra gets fairly heavy at times;
Redfearn (see bibliography) talks of a a particularly tough spot of work and Hotine talks
of reversing series by brute force and algebraso be warned. To make this article as
self-contained as possible I have added a number of appendices covering the required mathematics.
My sources for the TME formulae are to be found in Empire Survey Review dating from
the nineteen forties to sixties. The actual papers are fairly terse, as is normal for papers by
professionals for their peers, and their perusal will certainly not add to the details presented
here. Books on mathematical cartography are fairly thin on the ground, moreover they
usually try to cover all types of projections whereas we are concerned only with Mercator
projections. The few that I found to be of assistance are listed in the bibliography.
I would like to thank Harry Kogon for reading, commenting on and even checking the
mathematics outlined in these pages. Any remaining errors (and typographical slips) must
be attributed to myselfwhen you find them please send an email to the address below.
Peter Osborne
Edinburgh, 2008
peter.1@mercator.myzen.co.uk
Contents
1
Introduction
1.1
2.1
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
3.1
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4.1
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1
4.2
4.3
4.4
5.1
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
6.1
6.1
6.2
6.3
6.4
6.5
6.6
6.7
7.1
7.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1
7.2
7.3
8.1
8.1
8.2
8.3
8.4
8.5
8.6
8.7
8.8
8.9
9.1
9.1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.1
9.2
9.3
9.4
9.5
9.6
9.7
9.8
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3
. . . . . . . . . . . . . . . . . . . . . . . . . . . 9.5
. . . . . . . . . . . . . . . . . . . . . . . . . . 9.7
10.1
11.1
A.1
B.1
C.1
D Spherical trigonometry
D.1
E.1
F.1
G.1
R.1
X Index
X.1
Chapter
Introduction
1.2
Chapter 1. Introduction
The earliest accurate determinations of the figure of the earth were made by comparing
two high precision meridian arc surveys, each of which provided a measure of the distance
along the meridian per unit degree at a latitude in the middle of each arc. Two such measurements, preferably at very different latitudes, are sufficient to determine two parameters
which specify the ellipsoidthe major axis a together with the minor axis b or, more usually, the combination of the major axis with the flattening f (defined below). For example,
in the first half of the eighteenth century, French scientists measured a meridian arc of about
one degree of latitude in Lapland (crossing the Arctic circle) and a second arc of about
three degrees of latitude in Peru (crossing the equator) and confirmed for the first time the
oblateness of the ellipsoid. In 1830 Everest calculated an ellipsoid using what he took to
be the best two arcs, an earlier Indian Arc surveyed by his predecessor Lambton and once
again the arc of Peru. As more and longer arcs were measured the results were combined
to give more accurate ellipsoids. For example Airy discussed sixteen arcs before arriving at
the result he published in 1830:
a = 6377563.4m
b = 6356256.9m
f = 1/299.32
[Airy1830]
(1.1)
where the flattening f , defined as (a b)/a, gives a measure of the departure from the
sphere. Similarly Clarke used eight arcs to arrive at his 1866 ellipsoid:
a = 6378206.4m
b = 6356583.8m
f = 1/294.98
[Clarke1866]
(1.2)
Modern satellite methods have introduced global ellipsoid fits to the geoid, that for the
World Geodetic System of 1972 being
a = 6378135m
b = 6356750.5m
f = 1/298.26
[WGS72]
(1.3)
[GRS80]
(1.4)
and that for the Geodetic Reference System of 1980 (GRS80) being
a = 6378137m
b = 6356752.3m
f = 1/298.26
Topographic surveying
The aim of a topographic survey is to provide highly accurate maps of some region referenced to a specific datum. By this we mean a choice of a definite reference ellipsoid
together with a precise statement as to how the ellipsoid is related to the area under survey.
Chapter 1. Introduction
1.3
For example we could specify how the centre of the selected ellipsoid is related to the chosen origin of the survey and also how the orientation of the axes of the ellipsoid are related
to the vertical and meridian at the origin. It is very important to realize that the choice of
datum for any such survey work is completely arbitrary as long as it is a reasonable fit to
the geoid in the region of the survey. The chosen datum is usually stated on the final maps.
As an example, the maps produced by the Ordnance Survey of Great Britain (OSGB) are
defined with respect to a datum known as OSGB36 (established for the 1936 re-triangulation)
which is still based on the Airy 1830 ellipsoid which was chosen at the start of the original
triangulation in the first half of the nineteenth century. This ellipsoid is indeed a good fit to
the geoid under Britain but it is a poor fit everywhere else on the globe so it is not used for
mapping any other country. The OSGB36 datum defines how the Airy ellipsoid is related
to the the ground stations of the survey. Originally, in the nineteenth century, the origin
was chosen at Greenwich observatory but, for the 1936 re-triangulation no single origin
was chosen, rather the survey was adjusted so that the latitude and longitude of 11 control
stations remained as close as possible to their values established in the original nineteenth
century triangulation.
Until 1983, the United States, Canada and Mexico used the North American datum
established in 1927, namely NAD27. This is based on the Clarke (1866) ellipsoid tied to
an origin at Meades Ranch in Kansas where the latitude, longitude, elevation above the
ellipsoid and azimuth toward a second station (Waldo) were all fixed. Likewise, much
of south east Asia uses the Indian datum, ID1830, which is based on the Everest (1830)
ellipsoid tied to an origin at Kalianpur. The modern satellite ellipsoids used in datums such
as WGS72, GRS80, WGS84 are defined with respect to the Earths centre of mass and a
defined orientation of axes.
In all, there are two or three hundred datums in use over the world, each with a chosen
reference ellipsoid attached to some origin. The ellipsoids used in the datums do not agree
in size or position and a major problem for geodesy (and military planners in particular)
is how to tie these datums together so that we have an integrated picture of the worlds
topography. In the past datums were tied together where they overlapped but now we can
relate each datum to a single geocentric global datum determined by satellite.
Once the datum for a survey has been chosen we would traditionally have proceeded
with a high precision triangulation from which, by using the measured angles and baseline,
we can calculate the latitude and longitude of every triangulation station from assumed
values of latitude and longitude at the origin. Note that it is the latitude and longitude values
on the reference ellipsoid beneath every triangulation station that are calculated and used
as input data for the map projections. It is important to realise that once a datum has been
chosen for a survey in some region of the Earth (such as Britain or North America) then
it should not be altered, otherwise the latitude and longitude of every feature in the survey
region would have to be changed (by recalculating the triangulation data). But this has
already happened and it will happen again. For example the North American datum NAD27
was replaced by a new datum NAD83 necessitating the recalculation of all coordinates, with
resulting changes in position ranging from 10m to 200m. If (when) we use one of the new
global datums fitted by satellite technology as the basis for new maps then the latitude and
longitude values of every feature will change slightly again.
1.4
Chapter 1. Introduction
Cartography
A topographic survey produces a set of geographical locations (latitude and longitude) referenced to some specified ellipsoid. We are then faced with the problem of cartography,
the representation of the latitudelongitude data on the ellipsoid by a two-dimensional map.
There are an infinite number of projections which address this problem but in this article we
consider only the normal (N) and transverse (T) Mercator projections, first on the sphere (S)
and then on the ellipsoid (E). We shall abbreviate these projections as NMS, TMS, NME
and TME: they are considered in full detail in Chapters 2, 3, 6 and 7 respectively. At a
later date Chapter 10 may cover the oblique Mercator projection. Details of all these projections are given, without derivations, in the book by Snyder entitled Map ProjectionsA
Working Manual. (See bibliography).
We define a map projection by giving two functions x(, ) and y(, ) which specify
the plane Cartesian coordinates (x, y) corresponding to the latitude and longitude coordinates (, ). For the above projections, other than the oblique Mercator, the fundamental
origin is taken as a point O on the equator, the positive x-axis is taken as the eastward direction of the projected equator and the positive y-axis is taken as the northern direction of the
projected meridian through O. This convention agrees with that used in Snyders book but
beware other conventions! Many older texts, as well as most current continental sources,
adopt a convention with the x-axis as north and the y-axis as sometimes east and sometimes
west!
Chapter 1. Introduction
1.5
3. Isotropy of point scale. Ideally the scale factor would be isotropic (independent of
direction) at any point and as a corollary the shape of any small region would be
unalteredsuch a projection is said to be orthomorphic (right shape). By small we
mean that, at some level of measurement accuracy, the magnitude of the scale does
not vary over the small region. This condition is satisfied by the Mercator projections.
4. Conformal representation. Consider any two lines on the surface of the Earth which
intersect at a point P at an angle . Let P 0 and 0 be the corresponding point and angle
on the map projection. The map is said to be conformal if = 0 at all non-singular
points of the map. This has the consequence that the shape of a local feature (such
as a short stretch of coastline or a river) is well represented even though there will be
distortion over large areas. All of the Mercator projections satisfy this criterion.
5. Equal area. We may wish to demand that equal areas on the Earth have equal areas
on the projection. This is considered to be politically correct by many proponents
of the Peters projection but the downside is that such equal area projections distort
shapes in the large. The Mercator projections are not equal area projections and they
also distort shapes.
In summary the normal Mercator projection has the properties: (a) there are singular
points at the poles, (b) the point scale is isotropic (so the map is orthomorphic) but the
magnitude of the scale varies with latitude, being true on two parallels at most, (c) the
projection is conformal, (d) the projection does not preserve area. The transverse Mercator
projection has the properties: (a) there are singular points on the equator, (b) the scale
is isotropic (so the map is orthomorphic) with magnitude varying with both latitude and
longitude, being true on at most two curved lines which cannot be identified with parallels
or meridians, (c) the projection is conformal, (d) the projection does not preserve area.
1.6
Chapter 1. Introduction
give a scale factor of unity on the equator (and increasing as sec with latitude). Or again,
for a projection of Britain, the super-map would be 600km by 1200km and, if we demand
that the scale be unity on the central meridian at 2 W, we shall find that the scale factor
nowhere exceeds 1.001.
Once we have our mathematical super-map embodying a varying scale (but close to
unity) we can construct the actual maps for printing by a uniform scaling of the projection
coordinates by a constant representation factor (RF), constant that is for a given series of
maps. For example the OSGB uses 1:25000, 1:50000, 1:625000 and other values. Thus
the RF only arises at the printing stage and we can forget all about it and work with the
mathematical super-map for the theoretical analysis of the projections.
Finally, note the usage that a printed map is large scale when the RF, considered as
a mathematical fraction, is large and the map covers a small area. The OSGB 1:50000
maps are considered to be in this category and the 1:5000 series are of even larger scale.
Conversely small scale maps having a small RF, say 1:1000000 (or simply 1:1M), are used
to cover greater regions.
Chapter 1. Introduction
1.7
Historical note
Let us state at the outset that Gerardus Mercator (15121594) did not develop the mathematics that we shall present for his projection (NMS); moreover he had nothing at all
to do with three other projections that now carry his nameTMS, NME, TME. In 1569 he
published his map-chart entitled Nova et aucta orbis terrae descriptio ad usum navigantium
ementate accommadata which may be translated as A new and enlarged description of the
Earth with corrections for use in navigation. His full explanation is given on the map-chart:
In this mapping of the world we have [desired] to spread out the surface of the globe
into a plane that the places should everywhere be properly located, not only with respect to their true direction and distance from one another, but also in accordance with
their true longitude and latitude; and further, that the shape of the lands, as they appear
on the globe, shall be preserved as far as possible. For this there was needed a new
arrangement and placing of the meridians, so that they shall become parallels, for the
maps produced hereto by geographers are, on account of the curving and bending of
the meridians, unsuitable for navigation. Taking all this into consideration, we have
somewhat increased the degrees of latitude toward each pole, in proportion to the increase of the parallels beyond the ratio they really have to the equator. (Translation
from Fite and Freemansee bibliography).
This is an admirably clear statement and the last two sentences make clear his approach. In
order that the meridians should be perpendicular to the equator, and parallel to each other,
it is first necessary to increase the map length of a parallel as one moves away from the
equator. Now at latitude the circumference of a parallel is 2a cos and this must be
scaled up by a factor of sec so that the parallel and the equator have the same map length
(2a). Thus to preserve the shape of say a small rectangle at some latitude, projected from
ground to map, it is necessary to increase the meridian scale at that latitude by a factor of
sec . Exactly how Mercator produced his map is not known. He had had a good mathematical education but in 1559 he would not have had access to tables of the secant function
to aid him. So perhaps he simply drew rhumb lines on the globe from various points on the
equator, and at various azimuths, and took note of which locations on the globe lay on these
rhumb lines. He could then adjust the ordinate scale of his projection so that all the locations on any rhumb line on the globe lay on a straight line on the projection. Alternatively,
it has been suggested that he modified the parallels at ten degree intervals so that projected
rhumb lines were approximately straight for each ten degree interval. Whatever method he
used it is clear that he had grasped what was required, but his projection may have lacked
high accuracy.
Mercators prime aim was to construct a useful navigators chart but note that in the
above statement he was also concerned that shapes shall be preserved as far as possible;
he understood the desirability of an orthomorphic projection. He would have noticed the
distortion at high latitudes but he was probably satisfied that the appearance at temperate
latitudes (Europe) was really quite good. He may have regretted the distortions to landmass
shapes but this was far outweighed by the utility of his projection for navigators.
1.8
Chapter 1. Introduction
Mercator was extremely secretive about how he had produced his map-chart but this
only stimulated others to research their construction. The first to succeed, and publish (albeit
unwillingly), was a Cambridge professor of mathematics named Edward Wright (1558?
1615). Interestingly his 1599 publication is entitled The correction of certain errors in
navigation. He refers to certain errors in Mercators chart and explains how they could be
corrected by using a table of secants. More precisely he realized that if the parallel at latitude
has to move up by a factor proportional to sec then the net displacement of a parallel
from the equator would be given by summing all of the secants from the equator to the
parallel in question. In modern parlance this requires an integration of the secant function
(see equation 2.28). Of course Wright (like Mercator) was working in pre-calculus and prelogarithm days and integration of the secant function would not have been possible. Instead,
he presented his results in numerical tables of cumulative secants or meridional parts
derived by summation from secant tables evaluated at intervals of 1 minute of arc. In other
words he carried out a numerical integration with a very fine sub-division.. It is possible that
the errors he claims in Mercators chart were attributable to the fact that Mercators method
was equivalent to a much coarser approximation to the integration. Wrights tables certainly
allowed the construction of a very accurate chart based on a latitude and longitude values
taken from a rather fine globe modelled by his compatriot Emery Molyneux. For many
years thereafter the charts were widely described as Wright-Molyneux map projections but
the name of Mercator later became the standard appellation.
Wright gives a nice physical construction of the Mercator projection from a sphere.
Suppose a sphericall superficies with meridians, parallels, rumbes, and the whole hydrographical description drawne thereupon, to be inscribed into a concave cylinder,
their axes agreeing in one. Let this sphericall superficies swel like a bladder, (while
it is in blowing) equally always in every part thereof (that is, as much in longitude as
in latitude) till it apply, and join itself (round about and all alongst, also towards either pole) unto the concave superficies of the cylinder: each parallel on this sphericall
superficies increasing successively from the equinoctial [equator] towards either pole,
until it come to be of equal diameter with the cylinder, and consequently the meridians
still wideening themselves, til they become so far distant every where each from other
as they are at the equinoctial. Thus it may most easily be understood, how a sphericall
superficies may (by extension) be made cylindrical, . . .
Chapter 1. Introduction
1.9
called Henry Bond (16001678) stumbled on the numerical agreement between Wrights
tables and those for ln[tan()], as long as was identified with (/2 + /4). The mathematical proof of the equivalence immediately became noted as an important problem but it
was nearly thirty years before it was solved by James Gregory (16381675), Isaac Barrow
(16301677) and Edmond Halley (16561742) acting independently. These proofs eventually coalesced into direct integration of the secant function as presented in Chapter 2. The
modification of this integration for the ellipsoid (and NME) would be trivial. (But when?
By whom?)
Having given credit to Wright, Bond and others we must now remark that the first
to succeed, but not to publish, was almost certainly an English mathematician (and much
else besides) called Thomas Harriot (15601621). He left a large collection of unpublished manuscripts and many years late it became obvious that he had probably duplicated
Wrights calculation of meridional parts and moreover appears to have the link to the logarithmic tangent.
The transverse Mercator projection on the sphere was included in a set of seven new
projections published in 1772 by a continental (born in Alsace-Lorraine) mathematician and
cartographer, Johann Lambert (17281777). As we shall see in Chapter 3, the derivation
of this projection is a straightforward application of spherical trigonometry starting from
the normal Mercator result. Apparently Lambert even made some oblique references to the
transverse projection on the ellipsoid but it was Carl Friedrich Gauss (17771855) who first
constructed a conformal projection from the ellipsoid which preserves true scale on one
meridian, the projection we shall term TME. (This was in connection with the survey of
Hanover commenced in 1818).
Gausss method involved a double projection, from ellipsoid to sphere and then sphere
to plane. The modification of his work to construct a single equivalent projection was devel
oped only as late as 1912 by L Kruger.
For this reason the transverse Mercator projection
on the ellipsoid is often called the GaussKruger projection. This is the method we shall
examine.
The transverse Mercator projection was not much used until the middle of the twentieth
century when it was advocated for both the new British maps and the proposed world wide
system (UTM). In Britain the need for more precise series for the TME projections was
met by the papers of M Hotine, L P Lee and J C B Redfearn (see bibliography). The
last mentioned produced the most complete form and the solution is often referred to as the
Redfearn series. We shall derive these series in full detail.
1.10
Chapter 1. Introduction
The Mercator projection on the sphere (NMS) is defined as the single member of the class
which is such that an azimuth on the sphere and its corresponding grid bearing on the map
are equal. This property of conformality is then used to derive the projection formulae. The
scale is true only on the equator for basic NMS but we show how it may be modified to give
true scale on two parallels instead.
Chapter 3 discusses the transverse Mercator projection on the sphere (TMS). In this case
we are considering a projection onto a cylinder which is tangential to the sphere on a pair
of meridians which together form a great circle, such as 90 E and 90 W. These projections
are rather unusual when applied to the whole globe but in practice we intend to apply them
to a narrow strip on either side of the meridian of tangency which is then termed the central
meridian of the transverse projection. The crux is that by considering a large number of
such projection strips we can cover the whole sphere (except near the poles) with good accuracy. The derivation of the projection formulae is a straightforward exercise in spherical
trigonometry. An important new feature is that corresponding azimuths and grid bearings
are not equal (even though the transformation remains conformal) and we define their difference as the grid convergence. Finally we present low order series expansions for the
projection formulae.
Chapter 4 is the crunch. Our ultimate aim is to derive the projection equations for the
transverse Mercator projection on the ellipsoid (TME) in the form of series expansions.
The only satisfactory way of obtaining these results is by using a small amount of complex
variable theory. This method is complicated by both the geometrical problems of the ellipsoid and also by the fact that we need to carry the series to many terms in order to achieve the
required accuracy. Thus, for purely pedagogical reasons, in this chapter we use the complex
variable methods to derive the low order series solutions for TMS (derived in Chapter 3)
from the standard solution for NMS. That it works is encouragement for proceeding with
the major problem of constructing the TME projections from NME.
Chapter 5 derives the properties of the ellipse and ellipsoid that are required in later chapters. In particular we introduce (a) the principal curvatures in the meridian plane and its
principal normal plane, (b) the distinction between geocentric, geodetic and reduced latitudes, (c) the distance metric on the ellipsoid and (d) the series expansion which gives the
distance along the meridian as a function of latitude.
Chapter 6 derives the normal Mercator projection (NME) on the ellipsoid. The method
is a simple generalization of the methods used in Chapter 2 the only difference being in
the different form of the infinitesimal distance element on the ellipsoid. The results for
the projection equations are obtained in non-trivial closed forms. The inversion of these
formulae is not possible in closed form and we must revert to Taylor series expansions.
This chapter also contains a digression on double projections and includes a discussion
of the transformation of the ellipsoid to the sphere by means of the conformal latitude.
(The word double signifying that a second transformation from the sphere to the plane
is required to produce a map). As a corollary, the conformal latitude is used to provide a
second means of inverting the NME projection formulae.
Chapter 7 uses the techniques developed in Chapter 4 to derive the transverse Mercator
projection on the ellipsoid (TME) from that of NME. This derivation requires distinctly
Chapter 1. Introduction
1.11
heavy algebraic manipulation to achieve our main result, the Redfearn formulae for TME.
Chapter 8 returns to the definitions of point-scale factor and grid convergence presented
in Chapter 3 (for TMS) and derives the corresponding results for TME as series expansions.
Once again the algebra is fairly heavy.
Chapter 9 applies the general results of Chapters 7 and 8 for the TME projections to two
important cases, namely the Universal Transverse Mercator (UTM) and the National Grid
of Great Britain (NGGB). The former is actually a set of 60 TME projections each covering
6 degrees of longitude between the latitudes of 80 S and 84 N and the latter is a single
projection over approximately 10 degrees of longitude centred on 2 W and covering the
latitudes between 50 N and 60 N. We then discuss the variation of scale and grid convergence over the regions of the projection and also assess the accuracy of the TME formulae
by examining the terms of the series one by one. We find that for practical purposes some
terms may be dropped, as indeed they are in both the UTM and NGGB formulae. Finally
the projection formulae are rewritten in the completely different notation used in the OSGB
published formulae (see bibliography).
Finally, this work is continuing and two further chapters are intended:
Chapter 10 will cover the derivation of the oblique Mercator projections.
Chapter 11 will not be concerned with projections. It will discuss geodesics on the sphere
and ellipsoid and the problem involved in calculating arbitrary distances accurately. In
particular we will give the derivation of the Vincenty formulae for long geodesics on the
ellipsoid.
This article has tried to be as self-contained as possible and to this end there are seven
mathematical appendices. Many of these were developed for other uses so they are more
general in nature.
A Curvature in two and three dimensions.
B Inversion of series by Lagrange expansions.
C Plane Trigonometry.
D Spherical Trigonometry.
E Series expansions.
F Calculus of variations.
G Complex variable theory.
1.12
Chapter 1. Introduction
Chapter
2.1
Coordinates
The position of a point P on the sphere is denoted by an ordered pair (, ) of latitude, longitude values; the meridians ( constant), the equator ( = 0) and the small circles ( constant, nonzero) constitute the graticule on the sphere. The
figure shows a second point Q with coordinates
( + , + ), the meridians through P and
Q, arcs of parallels P M , KQ and the geodesic
(great circle) through the points P and Q. Such
geographical coordinates will be given in degrees
but in equations ALL angles will be in radians.
The unit mil, such that 6400mil=2 radians=360
is sometimes used for small angles, in particular
the grid convergence defined in Chapter 3.
Z
N
K
Q
M
Figure 2.1
(2.1)
2.2
If the distance of P from the axis (i.e. the radius of a parallel circle) is denoted by p(), so
that p() = a cos then the Cartesian coordinates are
X = p() cos = a cos cos ,
Y = p() sin = a cos sin ,
Z = a sin ,
(2.2)
Z
Z
= arctan
,
p
X2 + Y 2
= arctan
Y
X
.
(2.3)
For a point at a height h above the surface at P we simply replace a by a + h in the direct
transformations: the inverse relations for and are unchanged but they are supplemented
with the equation
h=
p
X 2 + Y 2 + Z 2 a.
(2.4)
In this chapter we consider a sphere with radius equal to the semi-major axis of the Airy
ellipsoid:
a = 6377563.396m 6378km 3963 miles .
This approximation value will suffice for the moment but more precise results will be needed
when we come to consider the large scale maps of TME. For the above radius the circumference of the equator (or any great circle) is approximately 40071km (24900 miles) and
the distance between pole and equator is one quarter of that value, 10018km (6225 miles).
The closeness of the latter value to 107 m reflects the original French definition of the metre
as 107 times the poleequator distance.
2.3
= 111.3km = 69.16miles,
00
= 31m
= 33.8yds,
= 1.312km = 0.815miles,
00
= 22m
1
1
= 23.9yds.
(2.5)
One minute of arc on the meridian (of a spherical Earth) was the original definition of the
nautical mile (nml). On the ellipsoid this definition of the nautical mile would depend
on latitude and the choice of ellipsoid so, to avoid discrepancies, the nautical mile is now
defined by international treaty as exactly 1852m, corresponding to 1.150779 miles. The
original definition remains a good rule of thumb for approximate calculations.
For two points in general position finding the distance is a non-trivial problem. There
are two important cases to consider: (a) given the geographic coordinates of two points
find the length of the geodesic between them and also the azimuths at the end points of the
geodesic joining them; (b) given a starting point and an initial azimuth find the coordinates
at a specified distance along the geodesic. These are the two principal geodetic problems
and their solution, for both sphere and ellipsoid, is presented in Chapter 11. In the present
chapter we consider only the problem of finding the distance between points on the surface
which are infinitesimally close. This will suffice for the calculation of scale factors.
Infinitesimal elements
In practical terms an element of area on the sphere can be said to be infinitesimal if, for
a given measurement accuracy, we cannot distinguish deviations from the plane. To be
explicit, consider the spherical element P M QK shown in Figure 2.1, and in enlarged form
in Figure 2.2a, where the solid lines P K, M Q, P Q are arcs of great circles, the solid lines
K
B
N
a
(a)
(b)
2 2
O
Figure 2.2
P M and KQ are arcs of parallel circles and the dashed lines are straight lines in three
dimensions. From Figure 2.2b, for (rad) 1 the arcchord difference is
1 3
a3
arc (AB) AB = a2a sin = a2a
+ =
+ O(a5 ). (2.6)
2
2 3! 8
24
Suppose the accuracy of measurement is 1m. Setting = we see that the difference
between the arc and chord P K will be less than 1m, and hence undetectable, if we take
2.4
(b)
a(cos)
M
+
Figure 2.3
We shall now prove that the small surface element P KQM may be well approximated
by a rectangular element. Figure 2.3a shows the planar trapezium which approximates the
surface element. Since P M = a cos we have
KQ P M =
d
(a cos ) = a sin .
d
(2.7)
1
P M KQ 1
= sin .
2
PK
2
(2.8)
QP
a cos
d
= cos
,
a
d
s2 = P Q2 = a2 2 + a2 cos2 2 .
(2.9)
(2.10)
(a cos )d,
ds2 =
dX 2 + dY 2 + dZ 2 ,
ds2 =
a2 d2 + a2 cos2 d2 .
(2.11)
(2.12)
2.5
2.2
The normal cylindrical projections of a sphere of radius a are defined on a cylinder of radius
a which is tangential to the sphere on the equator as shown in Figure (2.4). The axis of the
cylinder coincides with the polar diameter N S and the planes through this axis intersect
the sphere in its meridians and intersect the cylinder in its generators. The projection
=
y=af()
y
180 E
180 W
N
P(,)
=0
P(x,y)
O
x
S
a
x= a
Greenwich
meridian
x= a
(2.13)
y(, ) = a f (),
(2.14)
where and are in radians. With transformations of this form we see that the parallels on
the sphere ( constant) project into lines of constant y so that the orthogonal intersections
of meridians and parallels of the graticule on the sphere are transformed into orthogonal
intersections on the map; we shall see that this is not necessarily true for intersections at an
arbitrary angle. The spacing of the meridians on the projection is uniform but the spacing
of the parallels depends on the choice of the function f ().
1
2.6
Note that all normal cylindrical projections have singular points since the points N, S
at the poles transform into lines given by y = a f (/2). On the sphere meridians intersect
at the poles but on cylindrical projections meridians do not intersect. All other points of the
sphere are non-singular points. Of course there is nothing special about the poles; if we use
oblique or transverse projections the geographic poles are regular points and other points
become singularthe singularities at the poles are artifacts of the coordinate transformations. For example we shall find that the transverse Mercator projection has singular points
on the equator.
The equations (2.13, 2.14) define a projection to a super map of constant width equal to
the length of the equator, 2a or 40071km. Since the true length of a parallel is 2a cos ,
the scale (map length divided by true length) along a parallel is equal to sec , increasing
from 1 on the equator to infinity at the poles. Note that this statement about scale on a parallel applies to any normal cylindrical projection but the scale on the meridians, and at other
azimuths, will depend on f (). The actual printed projection in Figure 2.4 is about 8cm
wide on paper so the RF (representative factor) is 8cm/40071km which is approximately 1
to 500 million or 1:500M.
Angle transformations on normal cylindrical projections
In Figure 2.5 we compare the rectangular infinitesimal element P M QK on the sphere with
the corresponding element P 0 M 0 Q0 K 0 on the projection. The latter is also a rectangle but
without any approximation since the meridians map into lines of constant x and the parallels
map into lines of constant y. The angle K 0 P 0 Q0 is called the grid bearing corresponding
(a)
K sphere
s
a(cos)
Q +
(b)
K projection
s
P
x
Q y+y
y
M
x+x
Figure 2.5
to the azimuth on the globe. The geometry of these rectangular elements gives
(a)
tan =
a cos
a
and
so that
tan =
(b)
tan =
sec
tan .
f 0 ()
= 0
.
y
f ()
(2.15)
(2.16)
Note that = on the meridians (both zero) or on the parallels (both /2) but in general
6= unless f 0 () = sec . Therefore this is the condition for a normal cylindrical
projection to be conformal. It also provides the means of calculating f () for Mercator
projection.
2.7
(2.17)
QP
distance P 0 Q0 on projection
s0
= lim
,
QP s
distance P Q on sphere
(2.18)
QP
x2 + y 2
s02
=
lim
.
QP a2 2 + a2 cos2 2
s2
(2.19)
h = f 0 ().
(2.20)
Similarly, when P Q lies along a parallel of latitude, and y are zero and x = a. The
scale factor in this case is conventionally denoted by k. Therefore
parallel:
k = sec .
(2.21)
QP
which reduces to
x2 (1 + cot2 )
,
a2 cos2 2 (1 + cot2 )
sin
() = sec
,
sin
(2.22)
(2.23)
where we assume that has been found in terms of and from equation (2.16).
Area scale factor
The area scale is obtained by simply comparing the areas of the two rectangles P M QK
and P 0 M 0 Q0 K 0 . Denoting this scale factor by A and using equations (2.20) and (2.21).
A = lim
QP
x y
= sec f 0 () = hk.
(a cos )(a )
(2.24)
NB. All of these scale factors apply only to the normal cylindrical projections. They are
independent of , a reflection of the rotational symmetry.
2.8
2.3
x transformation
x-range
y transformation
y-range
f ()
f 0 ()
meridian scale (h)
parallel scale (k)
scale on equator
area scale (hk)
angles (eq. (2.16)
aspect ratio
equidistant
x = a
(a, a)
a
(a/2, a/2)
1
1
sec
1
sec
tan = sec tan
2
Figure 2.6
equal-area
x = a
(a, a)
a sin
(a, a)
sin
cos
cos
sec
1
1
tan = sec2 tan
Figure 2.7
Mercator
x = a
(a, a)
a ln[tan(/2+/4)]
(, )
ln[tan(/2 + /4)]
sec
sec
sec
1
sec2
tan = tan
0
Figure 2.8
Table 2.1
The three projections are shown in Figures (2.62.8) with every little smudge on the
maps indicating a small island or small lake! The maps all have the same x-range of
(a, a) (in metres) but varying y-ranges. They all have unit scale on the equator and the
RF on these pages is 12cm/40000km or about 1/300M.
On the equator = 0 so that (a) all the point scales (h, k, and ) are unity so
that the scale is isotropic and small elements retain their shape, (b) the area scale is unity
and (c) equation (2.16) shows that = so that all lines cross the equator on globe and
projection at the same angle. Thus on the equator the projections are perfectly behaved and
suitable for the accurate mapping of countries lying close to the equator.
Away from the equator the divergent parallel scale (sec as /2) produces
gross east-west stretching in high latitudes but the differing meridian stretching leads to
different shapes and areas. The acid test is to compare the projection with the outlines
shown on any toy or classroom globe. In particular look at the shape of Alaska and the area
of Greenland. The latter should be 1/8 that of South America and 1/13 that of Africa.
Finally, each of the projections is annotated on the right with a chequered column corresponding to 5 5 regions on the sphere. The width of these rectangles is the same on all
but their height depends on f ().
2.9
a/2
60
30
30
60
90
120
150
180
90
60
30
=0
y =0
30
a
a/2
x = a
60
x = a/2
x =0
x = a /2
90
x = a
by k = sec and the ruler distance must be divided by this factor and the RF to obtain
true distances. On a straight line on the projection with a constant bearing (not equal
to zero or /2), equation (2.23) shows that the scale is = sec sin () cosec with
() = arctan(tan cos ) obtained from (2.16). The dependence of this scale factor
means that there is no simple way of measuring along such a line with a ruler in order to
determine any kind of distance. A numerical integration of the scale factor along the line
on the projection is possible but it would give a distance on that line on the globe which
corresponds to constant; in general such a line is not a parallel, meridian, rhumb line or
great circle geodesic, so the utility of such an integration must be queried. The only simple
way to obtain the geodesic distance between general points on the projection is to transfer
their coordinates to the sphere and then use the standard geodesic formulae presented in
Chapter 11.
2.10
90
60
30
30
60
90
120
150
180
90
60
a/2
30
y=0
=0
a/2
30
60
90
a/2
x=0
a/2
2.11
90
60
30
30
60
90
120
150
180
85
2a
60
a
30
=0
y=0
30
a
60
2a
3a
a
85
a/2
x=0
a/2
2.12
of constant bearing on the map and vice-versa. For more on rhumb lines see Section 2.5.
The isotropy of scale implies that the shape of small elements is well rendered on the
projection; the property of orthomorphism. Consider a small square of side L at latitude
where the isotropic scale factor is sec . The projection will map this square into an exact
square of side L sec , preserving its shape, only if the variation of sec over the square can
be neglected at the precision of measurement in use. Larger shapes will be distorted if sec
varies appreciably over their size making the Mercator projection unsuitable for detailed
mapping of large countries.
The only true scale on the Mercator projection is attained at the equator and moreover
the area scale approaches unity too. Thus, since we still have conformality, it may well
be convenient to use Mercator for accurate large scale mapping near the equator. This is
discussed more quantitatively in the next section. Meanwhile we must accept the limitations
of using Mercator for small scale maps of the globecharts excepted!
Using ruler distance to deduce true distance is non-trivial away from the equator. On
parallels we can again divide the ruler distance by k = sec to find the true distance on the
parallel (which is not a geodesic difference). On a meridian where the scale is (h = sec )
we have to invert y = f () to find latitudes and hence an arc length: this is discussed in
the next section. On lines of constant bearing the scale factor is also sec . This can be
integrated but it gives a distance along a rhumb linesee Section 2.5
Conformal versus orthomorphic
The definition of conformal just given is that the projection preserves the angle between any two lines through any (non-singular) point of the map. This is an exact statement
referring to the behaviour of the transformations at a point. Many authors use the word
orthomorphic (right-shape) as a synonym for conformal but it is important to realise that
orthomorphism can only be satisfied approximately in a conformal projection because it is
a non-local property. No projection can be completely orthomorphic since even if small
shapes are preserved then one can find large shapes which are distorted. For this reason we
prefer to use the word conformal (exact) rather than orthomorphic (approximate).
2.4
f () =
sec d,
0
(2.26)
2.13
Now
cos = sin ( + /2)
= 2 sin (/2+/4) cos (/2+/4)
= 2 tan (/2+/4) cos2 (/2+/4)
so that
a
y = a f () =
2
Z
0
sec2 (/2 + /4)
.
d = a ln tan
+
tan (/2 + /4)
2
4
(2.27)
There is no need for modulus signs inside the logarithm. For /2 /2 the
argument of the tangent is in the interval [0, /2], therefore the argument of the logarithm
is in the range [0, ) and the logarithm itself varies from to .
Mercator parameter and isometric latitude
The function of which occurs in the expression for the y-coordinate in the Mercator projection is of importance in much that follows. We shall call it the Mercator parameter and
denote it by () and write the equations of the projection as
x(, ) = a,
with
y(, ) = a(),
() = ln tan
+
2
4
d
= sec .
d
(2.28)
The term Mercator parameter is not standard usage. It is usually called the isometric latitude. For example see Snyder(1987). But beware; other authors use the term for a different,
but related, function. Note, too, that the symbol is not universal: Lee(1946a) and Redfearn(1948) use , Maling(1992) uses q and so on.
The origin of the term isometric latitude relates to a re-parameterisation of the sphere
in such a way that the isometric latitude (replacing the geodetic latitude) and are involved
with equal weight in the metric formulae such as (2.12). If we write e for such an isometric
latitude then the metric ds2 = a2 d2 + a2 cos2 d2 becomes ds2 = a2 cos2 (de2 + d2 )
if we choose cos de = d. The coefficients in the metric are now both equal to a2 cos2 .
Since our is defined above (2.28) by cos d = d we see that the functions ()
e
and ()
must be identical. However the two functions are logically distinct and in an
elementary treatment we prefer to use different names but allow the same symbol .
Having urged care with notation we must flag a small problem in notation. When we
define the Mercator projection on the ellipsoid (NME in Chapter 6) we must define the
Mercator parameter in a slightly different way, but such that it reduces to (2.28) as the
eccentricity e tends to zero. It would therefore have been natural to define the Mercator
parameter on the sphere as 0 . We have not not done this, assuming that the correct interpretation will be obvious from the context.
2.14
1 + tan(/2)
cos(/2) + sin(/2)
=
1 tan(/2)
cos(/2) sin(/2)
(cos(/2) + sin(/2))2
1 + sin
=
= sec + tan .
2
2
cos
cos (/2) sin (/2)
(2.29)
Hence
() = ln[sec + tan ].
(2.30)
(2.31)
e = sec + tan .
After a little manipulation we find that
2 sinh = e e = 2 tan ,
and therefore
(a)
sinh = tan ,
(b)
cosh = sec ,
(c)
(2.33)
2.15
2.5
We have defined the rhumb line or loxodrome on the sphere as a curve which intersects
meridians at a constant azimuth and is therefore projected to a straight line on the Mercator projection with constant bearing such that = . Special cases are (a) the equator and
any parallel on which = = /2, (b) any meridian, on which = = 0. The utility
of a direct link between rhumb lines on the sphere and straight lines on Mercator charts is
obvious; it is discussed in many nautical publications and web-sites.
75
60
30
0
30
60
75
Figure 2.9: Rhumb line on the sphere and the Mercator projection
Figure 2.9 shows such a rhumb line crossing the equator at 30 W and maintaining a constant azimuth of 83 ; it spirals round the sphere from pole to pole whilst, on the projection,
it becomes a repeated straight line of infinite total length because the Mercator projection
extends to infinity along the y-axis. However, it is easy to show that corresponding rhumb
(a)
s
(b)
a(cos)
Figure 2.10:
line on the sphere has finite length. In Figure 2.10 we show an enlarged view of corresponding elements, P Q on the sphere and P 0 Q0 on the projection. The left hand triangle shows
that cos = a d/ds. Since is a constant this integrates directly to give
s12 = a sec (2 1 ),
1 6= 2 .
(2.36)
2.16
a rhumb line we need know only the constant azimuth and the change of latitude. Note that
the above result fails on the equator or any parallel since 1 = 2 . In this case we have
s12 = a cos (2 1 )
1 = 2 = .
It is very straightforward to plot the rhumb line on chart and on the sphere. On the chart it is
trivial for the equation of a line of gradient cot through the point (x1 , y1 ) on the projection
is simply
y y1 = cot (x x1 ).
(2.37)
The corresponding equation on the sphere follows immediately by using the Mercator projection formulae giving
a( 1 ) = a cot ( 1 ),
which, on substitution for from (2.28) , becomes
h
i
h
i
1
() = 1 + tan ln tan
+
ln tan
+
.
2
4
2
4
(2.38)
(2.39)
To plot the rhumb line on the surface of the sphere we choose a parallel () and then evaluate
at the point where the rhumb line crosses the parallel by using the above result.
Mercator sailing
The above equations solve the basic problem of Mercator sailing. That is, given a starting
point (1 , 1 ) and a destination (2 , 2 ) we have to find the azimuth of the rhumb line
and the sailing distance. The azimuth follows from (2.38) as long as we have tables of the
Mercator parameter or a calculator and the formula (2.28). Notice that in using (2.38) the
longitude difference must be expressed in radians by a preliminary calculation. Once we
have the azimuth then the sailing distance follows from (2.36) with the latitude in radians.
Of course, even before the advent of technology, there was no need for these calculation
as long as one was sailing short legs, no more than a few degrees, and as long as one
possessed a Mercator chart marked up with scales of latitude and longitude against which
we can plot start and final positions. The azimuth is trivial because it equals the bearing
marked on the chart as the line joining start and finish points.
The distance is a little more tricky because the chart scales vary non-linearly, but isotropically, as sec . To find the true distance ds, corresponding to ds0 we must first compensate
by dividing the ruler length of ds0 by sec . This will also be the case for finite distances
as long as the variation in sec along the sailing course can be neglected. Now rather than
working out ds0 / sec we can simply measure ds0 on a stretched rule, i.e. one which has
already been stretched by a factor of sec and there is one on the map already, namely
the latitude scale on the vertical edges of the chart. So all one need do is use dividers to
transfer ds0 to the latitude scale and if this is marked with minutes then we have our final
results in traditional nautical miles (on a sphere) without further calculation. For longer
voyages, more than a few hundred miles, reading the azimuth from the chart is still trivial
but the distance must be calculated using (2.36) because the latitude scale on the sides is too
non-linear to use as a measure over such distances.
2.17
Meridian parts
Equations (2.36) and (2.38) were the basic equations that navigators could use in the early
seventeenth century once they were furnished with Wrights table of Meridian parts. The
first edition was based on the division of the meridian into 540 parts of 100 interval but
he soon improved this to 5400 parts at 10 interval. For each part he calculated how much
it would be stretched by the Mercator scale factor and he then added these successively to
obtain his cumulative secants. Nautical texts often write the total of (stretched) Meridian
parts from the equator to a latitude as
MP() =
sec i
at intervals of 10 .
(2.40)
Bearing in mind that there are 3437.75 minutes of arc in one radian (Equation 2.1) we see
that the above sum is proportional to a good approximation to the integral of sec :
Z
MP() 3437.75
sec d = 3437.75 ().
(2.41)
0
Thus the meridional parts are simply proportional to the ordinate of the Mercator projection.
Now if we also express (2 1 ) in minutes of arc in equation (2.38) we obtain
tan =
(2 1 )0
MP(2 ) MP(1 )
(2.42)
where the MP values are obtained from Wrights tables. Thus, as long as we know latitude
and longitude values of start and finish, we can work out the required sailing course even if
we lack a chart. The sailing distance then follows immediately from equation (2.36) as
s12 = sec (2 1 )0 ,
(2.43)
where we have absorbed the factor of a so that if latitudes are expressed in minutes of arc
the result is in (traditional) nautical miles for a sphere.
2.6
We have already observed that since = for the Mercator projection the scale factor from
equation (2.23) is isotropic and equal to sec . It is customary to use k for the common value
of an isotopic scale factors therefore, in terms of the geographical coordinates, we have2
k(, ) = sec ,
(2.44)
(2.45)
Note the order. For all functions of two variables defined on the sphere we prefer to write , the absciisalike variable, first. This is unconventional but it makes subsequent chapters more logical.
2.18
Note that we should really distinguish these two scale factors, by using distinct functional
y), because replacing and in equation (2.44) by x and y
names such as k(, ) and k(x,
would give k(x, y) = sec x, which is not equation (2.45). This kind of distinction is rarely,
if ever, made in the literature of cartography and we shall not do so here; context will make
things clear.
Away from the equator the scale factor is k = sec > 1 and an infinitesimal ruler
distance on the map will therefore exaggerate distances on the sphere by this amount. It
is only on a parallel, where is constant, that we can measure large distances by dividing
the ruler distance by the RF and then dividing again by sec . Unlike the other projections
we have considerd, for Mercator we can also easily measure distances on a rhumb line by
using (2.36) but geodesic distances must be measured by the methods of Chapter 11.
Our main concern in this article is to investigate the application of the transverse Mercator projections to large scale accurate mapping so it is interesting to see how the normal
Mercator fares in this respect when we are close to the equator where k 1. It is necessary
to define the word accurate: we shall use it to mean that the scale variation is within 0.04%
of some specified value, corresponding to 4 parts in 10000. We shall call this the zone of
accuracy.
For normal Mercator the scale varies between k = 1 at the equator and 1.0004 as
latitudes 1 given by sec 1 = 1.0004, or cos 1 0.9996, corresponding to 1 = 1.62 .
Thus the zone of accuracy for the Mercator projection is a strip of about 3.24 width centred
on the equatorthis corresponds to 360km or 200 miles. The projection is certainly suitable
for accurate mapping that narrow band of latitudes; in the next section we shall see how the
projection may be modified to give a wider zone of accuracy.
We could argue that that the absolute value of the scale is not relevantonly the variation of scale over the mapped region is of interest. Consider, for example, the band of
latitudes starting at 10 N, where k = 1.015: we find that the scale has increased by 0.04%
when we reach 10.12 N so that the width of the zone of accuracy starting at 10 N is only
70 .7. This is a tiny strip indeed. By 30 N we have k = 1.155 and the zone of accuracy is
only 20 .4. Clearly such narrow zones are unsuitable and the Mercator projection must be
limited to a narrow equatorial zone for accurate mapping.
2.7
We have just seen that Mercator is accurate only within a narrow band centred on the equator. We now show how the width of this zone may be enlarged by making the scale on the
equator less than 1, but greater than 0.9996 so that we are still with the 0.04% tolerance. We
simply modify the projection cylinder by reducing its radius (Figure 2.11). If the cylinder
intersects the sphere in the parallels at latitudes 1 then the radius of the cylinder must be
a cos 1 . We then demand that the scale be true on the parallels 1 and we can achieve
this, and retain conformality, by multiplying both equations of the transformation by the
2.19
y=af()
P(x,y)
true scale
af(1)
x
1
true scale
y=af(-1)
x=k0 a
x=k0 a
+
,
y = ak0 ln tan
2
4
k0 = cos 1
(2.46)
The definition of the scale factor in equation (2.19) shows that it remains isotropic with a
value
k(, ) = k0 sec = cos 1 sec .
(2.47)
The scale factor is now k = cos 1 < 1 on the equator and increases with latitude. If
the lowest acceptable scale factor is k = 0.9996 then we must have cos 1 = 0.9996
corresponding to a value of 1 = 1.62 . Likewise, if the largest acceptable scale factor is
attained at = 2 then we must have
1.0004 = k(, 2 ) = k0 sec 2 = cos 1 sec 2 = 0.9996 sec 2 .
(2.48)
This equation gives cos 2 0.9992 and 2 = 2.29 , so that the projection is now
reasonably accurate in a strip of total width 4.58 centred on the equator. This corresponds
to a north-south distance of about 512km or 284 miles.
Thus for the modified normal Mercator projection the scale on the equator is not unity;
the parallels 1 , on which the scale is true, are called the standard parallels of the modified
projection. If we are willing to accept less accuracy then we can take the standard parallels
at higher latitudes. For example if we take 1 = 40 then the scale at the equator is
k = 0.76 and the latitudes at which k = 1.24 are 52 . Between these latitudes the
projection is accurate to within 24%. Similar considerations apply to all normal cylindrical
projections.
2.20
2.8
Mercator parameter
(2.49)
1
1+ sin
() = ln tan
= ln
,
+
2 4
2
1 sin
x
,
k0 a
(x, y) = 2 tan
(y) =
y
exp
,
k0 a
2
y
= sin1 tanh
,
k0 a
Conformality
Scale factors
y
k0 a
= ,
k(, ) = k0 sec ,
y
k(x, y) = k0 cosh
.
k0 a
(2.50)
(2.51)
(2.52)
(2.53)
(2.54)
(2.55)
(2.56)
Chapter
3.1
In Chapter 2 we constructed the normal Mercator projection (NMS). The strength of NMS
is its conformality, preserving local angles exactly and preserving shapes in small regions
(orthomorphism). Furthermore, meridians project to grid lines and conformality implies
that lines of constant azimuth project to constant grid bearings, thereby guaranteeing the
continuing usefulness of NMS as an aid to navigation.
As a topographic map of the globe, NMS has shortcomings in that the projection greatly
distorts shapes as one approaches the polesbecause of the rapid change of scale with
latitude. However, the (unmodified) NMS is exactly to scale on the equator and is fairly
accurate within a narrow strip of about three degrees centred on the equator (extending to
five degrees for modified NMS). It is this accuracy near the equator that we wish to exploit
by constructing a projection which takes a complete meridian great circle as a kind of
equator and uses NMS on its side to achieve a conformal and accurate projection within
a narrow band adjoining the chosen meridian. This is the transverse Mercator projection
(TMS) first demonstrated by Johann Lambert in 1772.
The modified versions of the transverse Mercator projection on the ellipsoid (TME, see
Chapter 7) are of great importance. They are used for map projections of countries which
have a predominantly north-south orientation, such as Great Britain. More importantly they
provide a systematic framework for covering the the whole of the globe with conformal
and accurate maps. The UTM (Universal Transverse Mercator) series covers the the globe
between the latitudes of 84 N and 80 S with 60 accurate projections of width 6 in longitude centred on meridians at 3 , 9 , 15 , . . . (The polar regions are always mapped with
projections centred at the poles).
3.2
Now for NMS the equator has unit scale because we project onto a cylinder tangential
to the sphere at the equator. Therefore, for TMS we seek a projection onto a cylinder which
is tangential to the sphere on some chosen meridian or strictly, a pair of meridians such as
the great circle formed by meridians at Greenwich and 180 E: the geometry is shown in
Figure 3.1a. This will guarantee that the scale is unity on the meridian: the problem is to
find the functions x(, ) and y(, ) such that the projection is also conformal.
'
$
&
!
"
%
#
Figure 3.1
The solution is remarkably simple. We first introduce a new graticule which is simply
the normal graticule of Figure 3.1a rotated so that its equator coincides with the chosen
great circle as in Figure 3.1b. Let 0 and 0 be the coordinates of P with respect to the new
graticule: they are the angles P CM 0 and OCM 0 on the figure. Note that, if 0 is measured
positive from M 0 to the rotated pole at E, then the sense in which 0 is defined on the
rotated graticule is opposite to the sense of in the standard graticule of Figure 3.1a. In
Figure 3.1b we have also shown x0 and y 0 axes which are related to the rotated graticule in
the same way that the axes were assigned for the normal NMS projection in Figure 2.4 so,
bearing in mind the sense of 0 , the equations (2.28) for NMS with respect to the rotated
graticule are
x0 = a0 ,
y 0 = a(0 ) = a ln tan 0 /2 + /4 .
(3.1)
Now the relation between the actual TMS axes and the primed axes is simply x = y 0 and
y = x0 , so that we immediately have the projection formulae with respect to the angles
(0 , 0 , ) of the rotated graticule:
x = a(0 ) = a ln tan 0 /2 + /4 ,
y = a0 .
(3.2)
It will prove more useful to use one of the alternative forms of the Mercator parameter, that
in equation (2.31), giving
a
1 + sin 0
x = ln
,
y = a0 .
(3.3)
2
1 sin 0
All that remains is to derive the relation between (0 , 0 ) and (, ) by a straightforward
exercise in applying spherical trigonometry to the triangle N M 0 P defined by the (true)
meridians through the origin and an arbitrary point P and by the great circle W M 0 P E
3.3
(Figure 3.2a). The right-hand figure shows a similar spherical triangle in standard notation
for which the sine and cosine rules (Appendix D) are (for a unit sphere)
sin A
sin B
sin C
=
=
,
sin a
sin b
sin c
(3.4)
(3.5)
(3.6)
(3.7)
=0
Figure 3.2
,
C ,
2
a 0 ,
b ,
c 0 ,
2
2
the first two terms of the sine rule and the first two cosine rules give
A ,
(3.9)
0
(3.8)
(3.10)
(3.11)
Note the simple expression for sin 0 in terms of and ; this explains why we chose the
alternative form of the NMS transformations in equation (3.3). To obtain the expression for
0 we eliminate cos 0 from the last two of these equations. On simplification we find
tan 0 = sec tan .
(3.12)
Therefore our final expressions for TMS centred on the Greenwich meridian are
a
1 + sin cos
x(, ) = ln
2
1 sin cos
y(, )
(3.13)
3.4
(0,-105)
(0,-120) (0,-135)
(0,+/-180)
(0,130.4)
(0,105.4)
(0,98)
(1.1,98.)
(7.3,93.0)
(15,-90)
(15,90)
(8.0,90)
(8.7,85.6)
y 0
(0,0)
(0,45) (0,60)
(0,75)
(0,82)
(-15,75)
-1
(-30,75)
(-45,75)
(-15,-90) (-30,-90) (-45,-90)
(-45,90) (-30,90)
(-15,90)
-2
-3
(0,-105)
-2
(0,-120)
(0,-120)(0,-135)
-1
3.2
3.5
(0,-175) (0,180)(0,170)
(0,180)
(0,155)
(0,110)
(0,60.4)
(0,35.4)
(0,28)
(1.1,28.)
(7.3,23.0)
(15,-160)
(15,20)
(8.0,20)
(8.7,15.6)
y 0
(0,-70)
(0,12)
(-15,5)
-1
(-30,5)
(-45,5)
(-15,-160) (-30,-160) (-45,-160)
(-45,20) (-30,20)
(-15,20)
-2
-3
(0,-175)
-2
(0,170) (0,155)
-1
sponding to the inverted meridian of 180 . On the other hand, the x-axis is now of infinite
length since E and W , the poles of the rotated graticule, are projected to infinity: we have
arbitrarily truncated the x-width at x = 2.5 in the figures.
Since all other meridians pass through N and S on the y-axis they are in general complicated curves running from top to bottom of the map. The exceptions are the meridians
at 90 which, since they are also great circles through the E and W poles, must extend
to infinity as horizontal lines on the map. The true parallels, except the equator, map into
a set of closed curves around the poles N and S; because of conformality these parallels
are orthogonal to the meridians at intersections, some of which are annotated with (latitude,longitude) geographical coordinates.
The equator itself appears on the map in three horizontal lines; the front equator lies
3.6
(0,45)
(0,-30)
(0,-79.6)
(0,-106.6)
(0,-112)
(1.1,-112)
(7.3,-117)
(15,60)
(15,-120)
(8.0,-120)
(8.7,-124.4)
(0,150)
y 0
(0,-165)(0,-150)
(0,-135)
(0,-128)
(-15,-135)
-1
(-30,-135)
(-45,-135)
(-15,60)
(-30,60) (-45,60)
-2
-3
(0,45)
-2
(0,30) (0,15)
-1
on the x-axis, extends to infinity and reappears as the back equator at the top and bottom
of the projection (corresponding to the generator along which the projection cylinder is cut
open). Note that longitude increases to the right on the x-axis and to the left at top and
bottom.
The grid lines which are parallel to the x-axis correspond to great circles which are
meridians of the rotated graticule. On the other hand the grid lines parallel to the yaxis have no readily defined precursors on the sphere. It is important to note that they are
parallel to the true meridians only where they cross the equator and nowhere else: the angle
between a (north-south) grid line and a curved meridian at a general point is called the grid
convergence, discussed in Section 3.6. We have also added some geographical coordinates
for grid intersections on the truncated sides and on the equator.
We shall give the formulae for the scale variation of this projection in Section 3.4.
3.7
Here we point out that the projection is quite faithful close to the central meridian but
there are gross distortions as we move away from the central meridianlook at Africa on
the projection centred on the Americas. This is now of no concern for we shall use the
projection centred on some 0 only near that meridian.
We know from Section 2.5 that a rhumb line, crossing meridians at a constant angle,
will reach both poles when extended on the sphere. Therefore, except for the zero azimuth
case, a rhumb line cannot be a straight line on the TMS projection. For short distances TM
maps are still suitable for navigation, particularly because conformality still guarantees a
simple relation (not equalitymore anon) between azimuth and grid bearing.
Finally note that we need not have shown the inverted sections of these maps for, by
suitable choices of central meridians, what is inverted on one map will be the right way
up on some other. Observe that in plotting these maps we cannot get the inverted sections,
where either y > a/2 or y < a/2, by using equation (3.13) with the arctan function in
its principal interval, (/2, /2). However arctan is a multivalued function, arbitrary to
within an additive factor of N . To plot the figures we used:
a
|| > /2, > 0
y=
+ a arctan [sec tan ]
for
.
(3.14)
0
|| < /2,
a
|| > /2, < 0
3.3
(3.15)
(3.16)
Eliminating gives
sec2 = sin2 coth2 (x/a) = 1 + cos2 tan2 (y/a),
tan2 coth2 (x/a) 1 = sec2 (y/a),
tan = sinh(x/a) sec(y/a),
(3.17)
(3.18)
(3.19)
(3.20)
(3.21)
where we take principal values in [/2, /2] to correspond to the front of the sphere.
3.8
The meridian distance m() a is the distance on the sphere measured along the
central meridian from the origin on the equator to a point at latitude . On the sphere
!"
!
m() = a.
(3.22)
The footpoint of P 0 (x, y) on the projection is that point P10 on the central meridian
of the projection which has the same ordinate as P 0 . The coordinates of the footpoint
are P10 (0, y).
The footpoint latitude, 1 , is the latitude of that point P1 on the central meridian
of the sphere which projects into the footpoint P10 (0, y). It is not the latitude of the
point P which is the inverse of P 0 .
=0
Figure 3.6
(3.23)
This is obvious because, by construction, the scale of the projection is true on the central
meridian so that y = a1 and hence 1 = y/a. In terms of the meridian distance function,
defined in equation (3.22) we see that 1 satisfies
m(1 ) = y.
(3.24)
We now take this equation as the definition of the footpoint latitude since we will find that it
continues to hold on the ellipsoid where m() is a non-trivial function. For future reference
we write equations (3.20) and (3.21) as
(x, y) = arctan [sinh(x/a) sec 1 ] ,
(x, y) = arcsin [sech(x/a) sin 1 ]
(3.25)
m(1 ) = y.
(3.26)
3.9
3.4
Because of the way in which the TMS was constructed, by applying NMS to a rotated
graticule, we know that the scale factor for TMS is isotropic and, in terms of the rotated
latitude 0 , its value is k = sec 0 . Using (3.9) and (3.15) we find the scale factor in terms
of either geographical or projection coordinates:
k(, ) =
1
(1 sin2 cos2 )1/2
(3.27)
k(x, y) = cosh(x/a).
(3.28)
[See comments after equation 2.45]. Note that the scale is a complicated function of the
geographical coordinates but is simply a function of the x-coordinate on the projection.
Both forms show that the scale is unity on the central meridian. ( = 0 or x = 0).
3.5
To investigate the relation between azimuths on the sphere and grid bearings on the projection we consider the relation of the infinitesimal elements shown in Figure 3.7. Now strictly,
an infinitesimal element on the projection would be a quadrilateral but we have drawn it as
curvilinear quadrilateral to emphasize the fact that the meridian M P , the parallel P N and
the displacement P Q will in general project to curved lines on the map. The relevant angles
must be defined with respect to the tangents of these lines at P 0 . The angles of concern are
(a)
(b)
Figure 3.7
, the angle between the meridian at P on the sphere and a general displacement P Q.
0 , the angle between the projected meridian and the projected displacement P 0 Q0 .
, the grid bearing, is the angle between the projected displacement and the y-axis.
, the angle between the projected meridian and the y-axis: this angle is the (grid)
convergence of the projection at P 0 ; in the next section we show how it is calculated.
Clearly 0 = + .
3.10
(3.29)
or, in words:
AZIMUTH
This equation is to be used in both directions. If we are given an azimuth at some point
on the sphere then the corresponding bearing on the map (chart) can be calculated from
= (, ). Likewise, given a bearing on the chart at (x, y) we find the azimuth
at the corresponding point on the sphere from = + (x, y). Clearly we need to find
expressions for the convergence in terms of both geographic and projection coordinates.
Although the convergence can take a wide range of values on small scale TMS projections (such as Figure 3.3), remember that the projection will be applied only in the region
very close to the central meridian where the the non-central meridian lines make very small
angles with the y-axis. For example, over Great Britain the convergence of the OSGB maps
is never greater than 5 .
3.6
The figure shows a section of the 45 E meridian between the equator and the north pole of the
/2
TMS projection of Figure 3.3. Since TMS is con
formal the angle between this projected meridian
(3.30)
so that in the quadrant shown in the figure, where x < 0 when y > 0, the convergence
is positive. (Thus a positive convergence is to the west of grid north).
3.11
Now the increments in x(, ) and y(, ) for arbitrary changes in and are
x
x
x =
+
,
y =
+
(3.31)
,
(3.32)
(3.34)
The partial derivatives must evaluated from equations (3.13); to put them in simpler forms
we use equation (3.9) and some equivalent forms
sin 0 = sin cos ,
(3.35)
(3.36)
Therefore
1 + sin cos
a
x = ln
,
2
1 sin cos
(3.37)
x
= a sec2 0 cos cos ,
y
= a sec2 0 sin sin cos ,
(3.38)
x
= a sec2 0 sin sin ,
y
= a sec2 0 cos .
(3.39)
(3.40)
This result can be written in terms of x and y by using equation (3.18) giving
(x, y) = arctan tanh(x/a) tan(y/a) .
(3.41)
It will prove useful to write this result in terms of x and the footpoint latitude as
(x, y) = arctan tanh(x/a) tan 1 ,
m(1 ) = y.
(3.42)
3.12
3.7
So far we have claimed, fairly, that TMS is conformal with an isotropic scale factor by virtue
of the method we used to derive the projection, viz NMS on its side. It is instructive to ask
how we may decide that an arbitrary projection from the sphere satisfies these conditions.
To this end consider Figure 3.7 where the azimuth angle of the displacement P Q on the
sphere is given by
cos
tan = lim
= lim R ,
(3.43)
QP
QP
where (for future developments) it is convenient to set R = cos . Now consider the grid
bearing of the corresponding displacement P 0 Q0 for an arbitrary projection. Using equations (3.31, 3.32), with the constraint implied by the above equation, we have
x + x
x tan + x R
x
tan = lim
= lim
=
.
(3.44)
y
y + y
y tan + y R
R
= tan
We already know tan from equation (3.33), therefore we can calculate 0 , the angle between the projected meridian and parallel, as
tan + tan
tan 0 = tan( + ) =
(3.45)
1 tan tan
=
=
y (x tan + x R ) x (y tan + y R )
y (y tan + y R ) + x (x tan + x R )
R (x2
(x y x y ) tan
.
+ y2 ) + (x x + y y ) tan
(3.46)
(3.47)
This is an identity which must hold for all values of , therefore the coefficient of tan and
the constant term must both vanish. This gives two conditions:
x x + y y = 0
R (x2 + y2 ) = (x y x y ).
Using (3.48), the second of these equations can be written as
R y2 1 + x2 /y2 = x y 1 + x2 /y2 ,
(3.48)
(3.49)
(3.50)
so that we must have x = R y . If we then substitute this back into (3.48) we obtain
y = R x . Thus, restoring R, the following conditions are necessary (and trivially
sufficient) for a conformal transformation from the sphere to the plane.
CAUCHY RIEMANN
x = cos y ,
y = cos x
(3.51)
It is trivial to check that these CauchyRiemann conditions are satisfied for both NMS and
TMS: in the first case have (from 2.28) x =a, y =a sec and x =y =0; TMS follows
immediately from equations (3.38, 3.39).
3.13
QP
= lim
QP
s02
x2 + y 2
=
lim
QP a2 2 + a2 cos2 2
s2
E2 + 2F + G2
.
a2 2 + a2 cos2 2
(3.52)
where
E(, ) = x2 + y2 ,
F (, )
G(, ) = x2 + y2 .
= x x + y y ,
(3.53)
An isotropic scale factor must be independent of the azimuth ; in other words (from 3.43)
it must be independent of the ratio of /. This is always the case when the Cauchy
Riemann equations (3.51) are satisfied,
for then we must have F = 0 and G = cos2 E.
1
k(, ) = a
x2 + y2 =
1
a cos
x2 + y2 .
(3.54)
Therefore ALL conformal transformations have isotropic scale factors. It is a simple exercise
to show that the above equation reduces to sec for NMS and to sec 0 for TMS; for the
latter use equations (3.353.39) to confirm the results of Section 3.4.
3.8
In the calculations for the transformations on the ellipsoid we shall have to resort to series
solutions. In this section we will derive the corresponding series for TMS. For the direct
transformations we hold constant and expand in terms of and for the inverse transformations we hold y constant and expand in terms of x/a. Typically the half-width (at the
equator) of a transverse projection is about 3 on the sphere and about 330km on the projection so that < .05 (radians) and x/a < 0.05. We shall drop terms involving fifth or
higher powers of these small parameters.
The coefficients of the direct series involve trigonometric functions of , which is not
generally a small term: for example tan is about 1.7 at 60 N. Likewise, the coefficients of
the inverse series will be functions of the footpoint latitude 1 which again is not generally
small. It is convenient to introduce the following compact notation for the trigonometric
functions of and 1 :
s = sin
c = cos
t = tan
s1 = sin 1
c1 = cos 1
t1 = tan 1
(3.55)
m(1 ) = y,
(3.56)
3.14
(3.57)
The argument of the arctan is not small but, using (E.16), we have
5 4
1 2
t sec = t 1 + + +
2
24
1 2
5 4
= t + z,
with z = t
+ + 1.
2
24
Using (E.9) with b = t we have
y(, ) = a arctan(t + z)
= a arctan(t) + at
1 2
5
+ 4
2
24
1
+ at2
1 + t2
1 2
+
2
2
(t)
.
(1 + t2 )2
asc 2 asc3 4
+
5 t2 + .
2
24
(3.58)
c
+
1
1
a 6 a3
3 a
1 2 x 3
1 x
1 1
= c1
+ c1
c
+
a
6 3 1
a
1 x (1 + 2t21 ) x 3
=
+ .
where 1 = y/a.
(3.59)
c1 a
6c1
a
3.15
1 x 2
5 x 4
s1 sech(x/a) = s1 1
+
+
2 a
24 a
1 x 2
5 x 4
= s1 + z
with z = s1
+
+
2 a
24 a
Using (E.8) with b = s1 we obtain
1 x 2
5 x 4
+
2 a
24 a
2
1 x 2
s1
1
2
s
+
+
+
2 (1 s21 )3/2 1
2 a
x 4
t 1 x 2 t1
+ (5 + 3t21 )
+ .
where 1 = y/a.
= 1
2 a
24
a
s1
(x, y) = arcsin s1 +
(1 s21 )1/2
(3.60)
(3.61)
1 x 2
1 x 4
+
+ .
2! a
4! a
(3.62)
(3.63)
3.16
Equation (3.41) is
(x, y)= arctan tanh(x/a) tan(y/a) = arctan t1 tanh(x/a) .
Expanding tanh(x/a) for small x with (E.23) and again using (E.20) for arctan gives
x
3
x 1 x 3
(x, y) = t1
+ (1/3)t31
+
a 3 a
a
x 1 t x 3
1
= t1
2
+
(3.64)
a
3 c1 a
3.9
Modified TMS
In Section 2.7 we showed how the NMS was modified to obtain greater accuracy over wider
areas by reducing the scale factor on the equator. We do the same for the TMS, reducing the
scale on the central meridian by simply multiplying the transformation formulae in equations (3.13) by a factor of k0 . The corresponding equations for the inverses, scale factors
and convergence are easily deduced: they are listed below along with the corresponding
series solutions. We continue to use the abbreviations for the trig functions of and 1
(equations 3.55, 3.56)
Direct transformations
1
1 + sin cos
1
x(, ) = k0 a ln
= k0 a c + c3 3 (1 t2 ) +
(3.65)
2
1 sin cos
6
sc 2 sc3 4
y(, ) = k0 a arctan [sec tan ] = k0 m()+k0 a
+
5t2 +
2
24
(3.66)
Note that on the central meridian given by = 0 we have y(, 0) = k0 m() = k0 a.
Therefore for the modified projection we must define the footpoint latitude by
1 =
Inverse transformations
x
y
(x, y) = arctan sinh
sec
k0 a
k0 a
y
.
k0 a
1
=
c1
(3.67)
x
k0 a
(1 + 2t21 )
6c1
x
k0 a
3
+
(3.68)
x
y
(x, y) = arcsin sech
sin
k0 a
k0 a
t1
= 1
2
x
k0 a
2
t1
+ (5+3t21 )
24
x
k0 a
4
+
(3.69)
3.17
Convergence
1
(, ) = arctan(tan sin )
= s + sc2 3 +
3
x 1 t x 3
x
y
1
(x, y) = arctan tanh
2
= t1
tan
+
k0 a
k0 a
a
3 c1 a
(3.70)
(3.71)
Scale factors
1 2 2
1 4 4
k0
2
= k0 1 + c + c (5 4t ) +
(3.72)
k(, ) =
2
24
(1 sin2 cos2 )1/2
"
#
x
1
x 2
1
x 4
k(x, y) = k0 cosh
= k0 1 +
+
+
(3.73)
k0 a
2! k0 a
4! k0 a
Consider the scale factor in terms of projection coordinates, that is k(x, y). If we choose
k0 = 0.9996 then we easily find that the scale is true when x/a = 0.0282 corresponding
to x = 180km (approximately). Once outside these lines the accuracy decreases as k
increases without limit. (The value of k reaches 1.0004 when x = 255km so that k increases
from 1 to 1.0004 in a distance of 75km. This is less than half of the distance over which the
scale changes from k = 0.9996 on the central meridian to k = 1 at x = 180km.)
Thus we see that the modified TMS is reasonably accurate over a width of approximately 510km. We shall see later that this includes most of the area covered by the British
grid.
3.18
Chapter
4.1
Introduction
(4.3)
In addition to these closed forms we also derived series expansions for TMS which neglect
terms of order 5 and higher:
1
x(, ) = ac + ac3 (1 t2 )3 + ,
6
y(, ) = a +
asc 2 asc3 4
+
5 t2 + .
2
24
(4.4)
(4.5)
The purpose of this chapter is to show how the above projection equations for TMS , both
closed forms and series, may be derived directly from the NMS projection equations. To
avoid confusion of the two sets of projection coordinates we shall reserve (x, y) for the TMS
projection and refer to the NMS projection by coordinates (, )note the order. For the
transformation from NMS to TMS we seek functions x(, ) and y(, ): for the inverse
transformations we seek two functions (x, y) and (x, y): the latter then gives (x, y) by
inverting one of (2.28, 2.32) or by a Taylor series.
4.2
Chapter 4. NMS to TMS by complex variables
(,)
(,)
Figure 4.1
Figure 4.1 summarizes the properties of the two projections. NMS is conformal, true
to scale on the equator (the -axis), with finite extent in and infinite extent in . TMS is
also conformal, true to scale on the central meridian (the y-axis), with finite extent in y and
infinite extent in x. A general point on the sphere P (, ) projects into points P 0 (, ) and
P 00 (x, y) for NMS, TMS respectively; a general point on the central meridian of the sphere,
taken as the Greenwich meridian for simplicity, projects into points K 0 (0, ) and K 00 (0, y).
We shall prove that the following conditions are sufficient to determine the functions x(, )
and y(, ) which define the transformation of NMS to TMS.
(4.6)
The scale on the y-axis of TMS is true so that the distance OK 00 on TMS is equal to
OK on the sphere; therefore y = m() = a.
The functions x(, ) and y(, ) must describe a conformal transformation: any
two lines through P 0 project into lines intersecting at the same angle at P 00 .
The first of these conditions is trivial but the others require further discussion.
The meridian distance as a function of
We need to consider the meridian distance as a function M of as well as the usual function
m of ; we equate M () and m() or, more strictly, we set
M = m () .
(4.7)
Thus the scale condition on the y-axis may be written as
y(0, ) = M .
(4.8)
On the sphere, where m() = a(), we can use equation (2.32a) to obtain an explicit
expression for M ():
M () = a arctan sinh() .
(4.9)
4.3
In our subsequent calculations we shall need the first four derivatives of M () with respect
to . These are straightforward enough to obtain as functions of from the last equation
but, at the end of the day, it will prove more useful to express the derivatives in terms of .
For example we have
M 0 ()
dM ()
dm() d
=
= a cos ,
d
d d
(4.10)
(4.11)
which is the definition of () in Section 2.4. Proceeding in this way we can construct the
first four derivatives of M () with respect to but with the results expressed as functions
of (using the compact notation for sin etc. defined in Section 3.8).
M 0 = a cos
d(a cos ) d
M 00 =
d
d
M 000 = a(c2 s2 )c
= ac3 (1 t2 ),
= asc,
ac,
asc3 (5 t2 ).
(4.12)
Figure 4.2
In Section 3.3, where we discussed the inverse TMS transformations, we introduced the
footpoint and the footpoint latitude. In considering the inverse transformations from TMS
to NMS it is useful to introduce the footpoint parameter 1 . All of these parameters are
indicated in Figure 4.2 which shows the (x, y) plane of TMS and the central meridians only
of NMS and the sphere. Given a point P 00 (x, y) the footpoint in the TMS plane is K 00 (0, y)
and the footpoint latitude 1 at K on the sphere is such that m(1 ) = y. The footpoint
parameter in the NMS plane is defined as the point K 0 (0, 1 ) such that
M (1 ) = y = m(1 ).
(4.13)
4.4
(a)
(b)
Figure 4.3
x = y ,
y = x
(4.14)
Satisfying the CauchyRiemann conditions and fitting the scale condition on the y-axis,
namely y(0, ) = M (), is a non-trivial problem. It becomes much more tractable when
we use complex numbers to effect the transformation. The basic idea is to associate with
the point P 0 (, ) of NMS a complex number = + i. Form a new complex number
by constructing a function z(); the real and imaginary parts of z are then used to define the
coordinates of a point P 00 (x, y) in TMS. Clearly this construction defines two real functions
x(, ) and y(, ). The theory of complex numbers tells us that if the function z()
is differentiable (or analytic, an equivalent term) then x and y must satisfy the Cauchy
Riemann conditions and define a conformal transformation.
There is a (very) concise introduction to complex functions in Appendix G. Here we
consider a trivial example of a conformal transformation. Consider z() = + 2 : this is
differentiable, giving z 0 () = 1 + 2 2 . Using i2 = 1 we have
z() = + 2 = (+i) + (+i)2 = +i + [2 +2i 2 ]
= + 2 2 + i( + 2).
(4.15)
Taking the real and imaginary parts defines the functions x(, ) = + 2 2 and
y(, ) = + 2 which satisfy the CauchyRiemann equations x = y = 1 + 2 and
y = x = 2. On the other hand this example does not satisfy the boundary conditions
x(0, )=0 and and y(0, )=M () of equations (4.6), (4.8) and (4.9).
4.5
Figure 4.4
4.2
The stages of the transformation are summarized in the above figure. We proceed anticlockwise from the sphere and derive the series solutions for TMS that were presented
in Chapter 3, albeit at a very low order which which would be inappropriate for accurate
mapping. The same steps will be used when we come to the ellipsoid projections.
Start at a general point on the sphere with coordinates P (, ) and map to the NMS
plane at P 0 (, ) with the Mercator parameter for the sphere.
Associate this point with the complex number = + i in the complex -plane.
Use a differentiable function z() to construct a conformal map from the complex
-plane to the complex z-plane.
Expand z() in a Taylor series about a point 0 on the imaginary axis (=0).
(4.16)
(4.17)
Demand that the scale be true on the central meridian in the z-plane.
y(0, ) = M ().
(4.18)
The result is a pair of series for x and y agreeing with those derived in Section 3.8.
Invert the Taylor series to find (z) and use the real and imaginary parts to find the
series for (x, y) and (x, y).
4.6
z0 = z(0 ) = iy0 = iM (0 ) = iM0
(4.20)
=+
Figure 4.5
It is instructive to consider the derivatives of z() from first principles, (as in (G.27):
z 0 () = lim
z( + ) z()
.
(4.21)
Now because z() is analytic we know that this limit is independent of direction and we
choose to take it in the direction so that = 0 and = i. Therefore we have
1 d
0
z () =
z().
(4.22)
i d
Once again, equations (4.17) and (4.18) imply that z() reduces to iM () at an arbitrary
point on the imaginary axis. therefore
d
0
z (0 ) = i
iM ()
= M 0 (0 ),
d
0
d
z 00 (0 ) = i
= iM 00 (0 ),
M 0 ()
d
0
d
000
00
z (0 ) = i
iM () = M 000 (0 ),
d
0
d
z 0000 (0 ) = i
M 000 () = iM 0000 (0 ).
(4.23)
d
0
Finally, if we abbreviate
be written as
z = z0 +(0 )M00
M 0 (
0 )=
M 00 ,
M 00 (
00
0 )=M0
i
1
i
(0 )2 M000 (0 )3 M0000 + (0 )4 M00000 +
2!
3!
4!
(4.24)
4.7
z = x + iy = iM +M 0
i 2 00 1 3 000 i 4 0000
M M + M +
2!
3!
4!
(4.25)
=+
()
Figure 4.6
The real and imaginary parts of equation (4.25) give x and y as functions of and . The
derivatives of M are real so that the transformations from NMSTMS are
1
x(, ) = M 0 3 M 000 +
(4.26)
3!
1
1
y(, ) = M 2 M 00 + 4 M 0000 + .
(4.27)
2!
4!
On substituting for M and its derivatives using equations (4.12), we obtain the corresponding expressions in terms of and (with s = sin etc. ) which define the transformation
from sphere to TMS:
1 3
ac (1 t2 )3 + ,
3!
1
1
y(, ) = a + asc 2 + asc3 (5 t2 ) 4 + .
2!
4!
These results agree with the expansions obtained in equations (3.57, 3.58).
x(, ) = ac +
(4.28)
(4.29)
y = M 0
(4.30)
(4.31)
4.8
Note also that the equations (4.28, 4.29) satisfy the the CauchyRiemann equations (3.51)
which apply to the transformation from the sphere to the TMS plane. Explicitly
1 3
ac (1 t2 )2 ,
2!
1
y = cos x = asc + asc3 (5 t2 ) 3 + .
3!
x =
cos y = ac +
(4.32)
(4.33)
(4.34)
where
b2 = i
M000
,
M00
b3 =
M0000
,
M00
b4 = i
M00000
.
M00
(4.35)
The series (4.34) and (B.13) are identical if we replace z and in the latter by (z z0 )/M00
and 0 respectively. Using (B.14) we can immediately find the inverse series of (4.34)
as
0 =
z z0
p 2 z z0 2 p 3 z z0 3 p 4 z z 0 4
+
M00
2! M00
3! M00
4! M00
(4.36)
p2 = b2
(4.37)
(4.38)
4.9
When we derived the inverse series in Section 3.8 we expanded and in power series in
x keeping y constant. The corresponding approach for the complex planes is indicated in
Figure 4.7. We start with a given point P 00 (x, y) on the projection corresponding to the point
z = x+iy in the complex z-plane. Let z0 = iy be the point K 00 in the z-plane corresponding
to the footpoint of P 00 . Clearly this point projects back to the central meridian of the -plane
=+
Figure 4.7
+
(4.39)
M1
2! M10
3! M10
4! M10
where the p coefficients are also evaluated at 1 using equations (4.37) with M00 M10
etc. Taking the real and imaginary parts (noting that p3 is real whilst p2 and p4 are pure
imaginary) we find
3
1
x
x
+
(4.40)
(x, y) = 0 p3
M1 3!
M10
2
4
x
x
1
1
(x, y) = 1 Im (p2 )
Im (p4 )
+ .
2!
M10
4!
M10
(4.41)
Now substitute for the p-coefficients from equations (4.38): these coefficients must now be
evaluated at the footpoint latitude 1 = (1 ) corresponding to the footpoint parameter 1 .
Since M10 = ac1 (equation 4.12) and we find
x 3
1 x
1 1
1 + 2t21
+ ,
(4.42)
(x, y) =
c1 a 3! c1
a
(x, y) = 1
x 4
1 t1 x 2
1 t1
+
5+6t21
+ .
2! c1 a
4! c1
a
(4.43)
The series for is in agreement with equation (3.59) but we must now derive the series for
from that for .
4.10
1
( 1 )2 sin 1 cos 1 + .
2!
(4.47)
+
s1 c1
2!
2! a
c1
which simplifies to
x 4
t1 x 2 t1
2
(x, y) = 1
5 + 3t1
+
+ ,
2 a
24
a
where m(1 ) = a1 = y, in agreement with equation (3.60).
4.3
(4.48)
Another way of deriving the inverse series is to take the development given in the first part
of Section 4.2 and run it backwards from the z-plane to the -plane. That is we assume the
existence of an analytic function (z) such that (a) the central meridian maps from x = 0
to = 0 and (b) on the -axis the scale is true. Therefore
(z) = (x, y) + i(x, y),
(4.49)
(0, y) = 0,
(4.50)
(0, y) = M (y),
(4.51)
4.11
where M (y) is an inverse to M () in the sense that M M (y) = y and M M () = .
The Taylor series analogous to (4.24) is then an expansion of (z) about a point on the
z0 = iy0 on the y-axis of the z-plane:
i
1
i
4
0000
(zz0 )2 M 000 (zz0 )3 M 000
0 + (zz0 ) M 0 + ,
2!
3!
4!
(4.52)
where 0 = i0 = iM (y0 ) and the derivatives of M are with respect to y at y0 .
(z) = 0 +(zz0 )M 00
(4.53)
= M (M ()) = M (y),
(4.54)
to give
dy
= M 0 (),
d
d
= M 0 (y).
dy
(4.55)
d
1
= 0
,
dy
M ()
(4.56)
and in general
d( )
1
d( )
= 0
.
dy
M () d
(4.57)
It is now straightforward to calculate all the derivatives in equation (4.52). For compactness
we suppress the argument in M () and all its derivatives, M 0 (), M 00 () etc. Comparing
the results with the p-coefficients in equation (4.37) we find
M 0 (y) =
1
M0
1 d h 0 i
1 d
1
M 00
M (y) = 0
M (y) = 0
=
M d
M d M 0
(M 0 )3
1 d h 00 i
1 d
M 00
M 000
(M 00 )2
000
M (y) = 0
M (y) = 0
=
+
3
M d
M d (M 0 )3
(M 0 )4
(M 0 )5
1 d h 000 i
M 0000
M 00 M 000
(M 00 )3
M 0000 (y) = 0
M (y) =
+
10
15
M d
(M 0 )5
(M 0 )6
(M 0 )7
00
1
,
M0
ip2
,
M 02
p3
,
M 03
ip4
.
M 04
(4.58)
Substituting these derivatives (evaluated at y0 for M and at 0 for M ) into equation (4.52)
clearly gives a complex series identical to that of (4.36) and the same results follow. We
choose not to follow this method since the calculation of the derivatives to eighth order for
the ellipsoid becomes very intricate.
4.12
4.4
Finally, we present a closed analytic expression whose real and imaginary parts give the
TMS transformations of equation (4.2). This section does not readily generalise to the
transformations on the ellipsoid: it is included as an interesting digression.
Finding a conformal transformation which satisfies given conditions is not always simple: there is no general technique and we have to rely mainly on our experienceand the
fact that there are books which give exhaustive lists of the conformal transformations which
have been studied in the last two hundred years. The required transformation is
z=
ia
2ia arccot [exp(i)]
2
(4.59)
We now verify that this transformation has the required properties. Substituting z = x + iy
and = + i in the above gives
x + iy
y
ix
cot
+
= cot
+
= exp(i + ).
(4.60)
2ia
4
4 2a 2a
The real and imaginary parts of the cotangent function are given in Appendix G, equay
x
tion (G.17). We substitute A =
4 2a and B = 2a in that identity but to clarify the
algebra we temporarily replace x/a and y/a by x and y respectively. The result is
cos y i sinh x
= e (cos i sin ).
cosh x sin y
Taking the real and imaginary parts gives
cos y
= e cos p
cosh x sin y
sinh x
= e sin q.
cosh x sin y
(4.61)
(4.62)
(4.63)
(4.64)
(4.65)
(4.66)
Eliminate y from (4.65) and (4.66) using cos2 y + sin2 y = 1; eliminate x from (4.64) and
(4.66) using cosh2 x sinh2 x = 1. On simplification we find
2q
2e sin
=
= sin sech ,
p2 + q 2 + 1
e2 + 1
p2 + q 2 1
e2 1
tan y =
=
= sec sinh .
2p
2e cos
tanh x =
(4.67)
(4.68)
This is the final result for the real and imaginary parts of the transformation from the complex -plane to the complex z-plane.
4.13
The final step is to transform from to in (4.67, 4.68) . To do this we use the inverse
expressions from equation (2.32), that is
sinh = tan ,
cosh = sec ,
(4.69)
(4.70)
(4.71)
The second equation agrees directly with the y transformation equation (4.2) after restoring
y y/a. Equation (4.3) for the x transformation follows since
1 + sin cos
1 + tanh x
cosh x + sinh x
2ex
=
=
= x = e2x e2x/a .
1 sin cos
1 tanh x
cosh x sinh x
2e
(4.72)
cos sech ,
2
(4.73)
(4.74)
sech2 x = 1 sin2 sech2 = sech2 cosh2 sin2 ,
sec2 y = 1 + sec2 sinh2 = sec2 cosh2 sin2 .
(4.75)
(4.76)
4.14
Chapter
5.1
We now model the Earth as an ellipsoid of revolution for which the Cartesian coordinates
with respect to its centre satisfy
X2 Y 2 Z2
+ 2 + 2 = 1.
(5.1)
a2
a
b
The definition of longitude is exactly the same as on the sphere. The geodetic latitude ,
which we will simply call the latitude, is the angle at which the normal at P intersects the
equatorial plane (Z = 0). The crucial new feature is that the normal does not pass through
Figure 5.1
the centre of the ellipsoid (except when P is on the equator and at the poles). The line
joining P to the centre defines the geocentric latitude c . We introduce the notation p()
for the distance P N of a point P from the central axis and we also set () for the length
CP of the normal at P to its intersection with the axis. Therefore
p
p() = () cos = X 2 + Y 2 .
(5.2)
5.2
5.2
Instead of using (a, b) as the basic parameters of the ellipse we can use either (a, e) where
e is the eccentricity, or (a, f ) where f is the flattening. These parameters are defined and
related by
ab
,
e2 = 2f f 2 .
(5.3)
a
For numerical examples we use the values for the Airy (1830) ellipsoid which is used for
the OSGB maps:
b2 = a2 (1 e2 ),
f=
a = 6377563.396m,
e = 0.0816733724,
b = 6356256.910m,
e2 = 0.00667053982,
f = 0.0033408505,
1
= 299.3249753.
f
(5.4)
The flattening of the Earth is small. For example, in the figures on the previous page the
difference between a sphere of radius a and the ellipsoid should be about the width of one of
the lines in the figure. Thus the ellipses shown here, and elsewhere, are greatly exaggerated.
Other parameters used to describe an ellipse
There are several other small parameters which arise naturally in the study of the properties
of the ellipse. Two which we shall need are: (a) the second eccentricity, e0 , and (b) the
parameter n (sometimes e1 ). They are defined by
a2 b2
e2
ab
=
,
n = e1 =
.
(5.5)
2
2
b
1e
a+b
There are many possible relations between all these parameters. For example we will need
the following results:
1+n
2 1/2
a=b 1e
=b
(5.6)
1n
= b 1 + 2n + 2n2 + 2n3 + ,
(5.7)
e02 =
2
1n 2
4n
b
e =1
=1
=
a
1+n
(1 + n)2
2
5.3
(5.8)
(5.9)
The equation of the cross-section ellipse follows from (5.1) and (5.2):
p2 Z 2
+ 2 = 1.
a2
b
(5.10)
5.3
(5.11)
Since the normal and tangent are perpendicular the product of their gradients is 1 and
therefore the gradient of the normal is
1
dZ
Za2
Z
tan =
= 2 =
.
(5.12)
dp
pb
p(1 e2 )
Substituting Z = p(1 e2 ) tan in (5.10) gives
p2 [1 + (1 e2 ) tan2 ] = a2 .
(5.13)
P M = Z() =
a cos
[1 e2 sin2 ]1/2
a(1 e2 ) sin
[1 e2 sin2 ]1/2
(5.14)
(5.15)
a
[1
(5.16)
e2 sin2 ]1/2
in terms of which
p() = () cos ,
(5.17)
Z() = (1 e2 ) () sin .
The triangle OCE
Later we shall require the sides of the triangle OCE
defined by the normal and its intercepts on the axes.
OE = OM EM = p Z cot
= cos (1 e2 ) cos
= e2 cos .
=()
CE = e2
OC = e2 sin .
(5.18)
Figure 5.2
(5.19)
The sides of this small triangle are all of order ae2 ; for example at latitude 45 the sides
OE and OC are about 30km and CE is about 42.5km.
5.4
(5.22)
Therefore
c has a turning point, clearly a maximum, when the right hand side vanishes
at tan = 1/ 1 e2 . Using the value of e for the Airy ellipsoid (equation 5.4) shows that
the maximum difference occurs at 45 .095, for which c 44 .904 and the latitude
difference c 11.50 . (Note that e2 0.00667 is the radian measure of 22.90 ).
A comment on other latitudes
In addition to the geodetic latitude and geocentric latitude c , we have already discussed
the isometric latitude (Section 2.4) and we shall meet three further latitude definitions:
the reduced (or parametric) latitude U (Section 5.5), the rectifying latitude (Section 5.9)
and the conformal latitude (Section 6.5). With the exception of the isometric latitude all
of these latitudes coincide with the geodetic and geocentric latitudes at the poles and on the
equator and the maximum deviations from are no more than a few minutes of arc. The
isometric latitude agrees with the others at the equator only but diverges to infinity at the
poles: it is a radically different in character.
5.4
Using (5.17) and (5.18) the Cartesian coordinates of a point on the surface are
X() = p() cos = () cos cos ,
(5.23)
(5.24)
(a) = arctan
,
(b) = arctan
.
X
(1 e2 ) X 2 + Y 2
(5.25)
(5.26)
5.5
=()
Figure 5.3
For the inverse relations dividing equation (5.28) by (5.27) gives explicitly, as in equation (5.26a). To find and h we can eliminate from (5.27) and (5.28) and rewrite equation (5.29) for Z to give
p
X 2 + Y 2 = () + h cos ,
(5.30)
2
Z + e () sin = () + h sin .
(5.31)
Dividing these equations gives an implicit equation for :
"
#
Z + e2 () sin
p
= arctan
.
X2 + Y 2
(5.32)
There is no closed solution to this equation but we can develop a numerical solution by
considering the following fixed point iteration:
"
#
Z + e2 (n ) sin n
p
n+1 = g(n ) = arctan
,
n = 0, 1, 2 . . . .
X2 + Y 2
(5.33)
Now in most applications we will have h a so that a suitable starting approximation is
the value of obtained by using the h = 0 solution, equation (5.26b):
Z
0 = arctan
.
(5.34)
(1 e2 ) X 2 + Y 2
If the iteration scheme converges so that n+1 and n in (5.33) then must
be the required solution of equation (5.32). The condition for convergence of this fixed
point iteration is that |g 0 ()| < 1: this is true here since g 0 () = O(e2 ). Once we have
found it is trivial to deduce h from equation (5.30):
p
h = sec X 2 + Y 2 ().
(5.35)
Comment: The formulae of this section are presented without proof in the OSGB web
publication mentioned in the Bibliography. We do not use them elsewhere but the method
of iteration to a fixed point is important and we will apply it in several places.
5.6
5.5
There is another important and obvious parameterisation of the ellipse. Construct the auxiliary circle of the
ellipse: it is concentric and touches the ellipse at the ends
of its major axis so that the radius is equal to a. Take
a point P on the ellipse and project its ordinate until it
meets the auxiliary circle at P 0 and let angle P 0 OA be U .
The angle U is called the reduced latitude (or parametric
latitude) of the point P on the ellipse.
The points P and P 0 clearly have the same abscissa,
p = a cos U . Substituting this abscissa into the equation
of the ellipse (5.10) we have
p
Z = b 1 p2/a2 = b sin U.
Figure 5.4
(5.36)
Z = b sin U,
(5.37)
constitutes the required parametric representation of the ellipse. It is clear that that the
ellipse is related to the auxiliary circle by scaling in the Z direction by a factor of b/a.
Relations between the reduced and geodetic latitudes
Comparing the parameterisations of p and Z in equations (5.17, 5.18) and (5.37) gives
p = () cos
= a cos U,
(5.38)
but it is more useful to divide the expressions for Z and p to find (using b = a 1 e2 )
tan U = 1 e2 tan
(5.39)
It will also be useful to have an expression for in terms of U . Using (5.38) and (5.16)
e2 2
e2 cos2
2
cos
=
1
a2
1 e2 sin2
1 e2
=
,
1 e2 sin2
1 e2 cos2 U = 1
(5.40)
so that
1/2
a
a
2
2
=
=
1
e
cos
U
1/2
1 e2
1 e2 sin2
(5.41)
5.7
(5.42)
(5.43)
occurred when tan = 1/ 1 e2 we deduce that the maximum value of U will occur
when tan = 1/ 4 1 e2 . This corresponds to 45 .048 for which the corresponding
value of U is 44 .952 so that the maximum difference is U 50 .7.
$
#
$
!
"
!
"
&
'
%
#
%
&'
%
&'
5.6
We now investigate the properties of the two dimensional curves formed by the intersection
of some, but not all, planes with the surface of the ellipsoid: we use the mathematical results
established in Appendix A. In particular we investigate two special families of planes.
The first family (S) has the normal at P as a common axis and the intersections of its
planes with the surface are called the normal sections at P . One member of the family is
Figure 5.5
the meridian plane (black) containing P and the symmetry axis of the ellipsoid. Another
important member of the family is the plane at right angles to the meridian plane: it is called
the prime vertical plane (shaded grey). Other members of the family are labelled by the
angle between a specific plane and the meridian plane.
The second family of planes (T) has as its axis the tangent to the parallel circle at P : we
are interested in just two of its planes. One (black) is the plane of the parallel: its section
on the surface is the parallel circle. The other is that which contains the normal at P : this is
the prime vertical plane (grey), the only plane common to both families.
5.8
1
1 e2
=
.
(5.44)
a [1 e2 cos2 U ]3/2
It will be more useful to work with the meridian radius of curvature defined by = 1/
and expressed as a function of . Using equation (5.41) we have
() =
a(1 e2 )
3/2
1 e2 sin2
(5.45)
3
1 e2 .
2
a
(5.46)
Furthermore, we have
1 e2
=
.
(5.47)
1 e2 sin2
Since (a) the denominator is less than or equal to 1 and (b) the numerator is less than or
equal to the denominator, we have
(1 e2 ) .
(5.48)
(5.49)
We identify the prime vertical plane and parallel plane of the family T with the normal and
oblique planes of the theorem. Now the parallel plane intersects the surface in a parallel
circle so we know that its radius of curvature is simply N P = p() in Figure 5.2. But this
is just () cos and therefore
Rprime vertical = sec Rparallel = sec p() = ().
(5.50)
5.9
Thus we have the important result that the distance CP = () may be identified as the
radius of curvature of the normal section made by the prime vertical plane. The point C
where the normal meets the axis is the centre of curvature of this section.
Radius of curvature along a general azimuth
Returning to S, the family of planes on the normal, we now know the curvature of two
of the normal sections: 1 on the meridian plane and 1 on the prime vertical. Now
consider the curvature, K(), of the section made by that plane of the family at an angle ,
measured clockwise from the meridian plane. Clearly the symmetry of the ellipsoid about
any meridian plane implies that K() = K() so that K() is a symmetric function of
and it must therefore have a turning point at = 0. Therefore the curvature of the meridian
section must be either a minimum or maximum and it is therefore one of the principal
curvatures at P see Appendix A.
In the appendix we proved that the planes containing the principal curvatures are orthogonal. Therefore the curvature of the normal section made by the prime vertical plane
must be the other principal curvature. Furthermore, equation (5.47) gives and
therefore 1 1 so that the curvature in the meridian section is the maximum normal
section curvature at any point. Introducing the radius of curvature on the general section by
R() = 1/K(), we use Eulers formula, equation (A.36), to deduce that
1
1
1
= cos2 + sin2 .
R()
(5.51)
1 e2 sin2
=
.
1 e2
1=
e2 cos2
.
1 e2
(5.53)
We shall frequently require the derivatives of the curvatures and their quotient. It is straightforward to show that
d
= ( 1) tan ,
d
d
( 1)
=3
tan ,
d
d
= 2( 1) tan .
d
(5.54)
1 d2
( 1)
1
2
=
+
2
5
+
3
tan2 . (5.55)
d2
5.10
Finally we note that the cross-section coordinates (5.17, 5.18) and their derivatives are
Z() = (1 e2 ) () sin ,
p() = () cos ,
dp
= sin ,
d
(5.56)
dZ
= cos .
d
(5.57)
Spherical limit
We shall refer to the limit e 0 as the spherical limit. Clearly in this limit
a,
5.7
a,
0 , 0 , 0 0.
1,
(5.58)
(5.59)
we have
dX = p cos d p sin d,
where DOT
dY = p sin d + p cos d,
dZ = Z d.
d
d
(5.60)
(5.61)
2
ds = d + cos d .
(5.62)
(5.63)
(5.64)
5.11
The infinitesimal element on the sphere was discussed in Section 2.1 and the same considerations apply on the ellipse. In particular the infinitesimal element may be approximated
by a planar rectangular quadrilateral to which we can apply plane trigonometry. Equations (5.63) and (5.64) show that the sides are equal to on the meridians and cos
Figure 5.6
on a parallel. The metric (5.62) may be viewed as the application of Pythagoras to the
infinitesimal element. The azimuth of an arbitrary displacement is calculated from
tan =
cos d
.
(5.65)
smeridian =
()d = a(1 e )
dsmeridian =
1
5.8
3/2 . (5.67)
1 e sin2
2
On the sphere we defined the meridian distance m() as the distance along a meridian from
the equator to a point at latitude ; this was trivially m() = a. On the ellipsoid we use
the same notation but the definition follows from equation (5.67):
Z
Z
Z
d
2
m() =
dsmeridian =
()d = a(1 e )
(5.68)
3/2 .
0
0
0
1 e2 sin2
5.12
m() = a(1 e )
1 + b2 e2 s2 + b4 e4 s4 + b6 e6 s6 + b8 e8 s8 + d,
(5.69)
b4 =
15
,
8
b6 =
35
,
16
b8 =
315
.
128
(5.70)
Using the trigonometric identities (C.32) to (C.38) we can express the sin2 , . . . sin8 in
terms cos 2, . . . cos 8. Collecting terms with the same cosine factors and integrating will
then give a series starting with a term and followed by terms in sin 2, . . . sin 8. The
result is:
m() = A0 + A2 sin 2 + A4 sin 4 + A6 sin 6 + A8 sin 8 + ,
(5.71)
2
8
16
128
4 64 256 161024
a(1e2 )
b2 e2 b4 e4 15b6 e6 7b8 e8
3e2 3e4 45e6
420e8
A2 =
=a
2
2
2
32
16
8
32 1024 161024
15e4 45e6
525e8
a(1e2 ) b4 e4 3b6 e6 7b8 e8
+
+
=a
+
+
,
A4 =
4
8
16
32
256 1024 161024
a(1e2 )
b6 e6 b8 e8
35e6
175e8
A6 =
=a
,
6
32
16
3072 121024
a(1e2 ) b8 e8
315e8
A8 =
=a
.
(5.72)
8
128
1281024
If we use the numerical values for the Airy ellipsoid (5.4), then, in metres,
m() = 633691460915979859 sin 2+16711 sin 40022 sin 6+000003 sin 8
(5.73)
The first four terms have been rounded to the nearest millimetre whilst the last term shows
that the O(e8 ) terms give rise to sub-millimetre corrections. We shall therefore drop O(e8 )
terms in expressions for the meridian distance .
The distance from equator to pole is defined by
1
mp = m(/2) = A0 = 10001126081 metres.
2
(5.74)
5.13
a2 =
15
,
8
a3 =
35
.
16
(5.76)
Apart from the overall constant multiplier the integral becomes, to O(n3 ),
Z z
dz
1
1
1
2 2
3
2 2
3 3
1 + a1 n + a1 n+a1 a2 n
z+
+ a2 n z + 2 + a3 n z + 3
z
z
z
1 2iz
z
1
1
a2 n2 2 1
a3 n3 3 1
2 2
3
=
1 + a1 n ln z+ a1 n+a1 a2 n
z
+
z 2 +
z 3
2i
z
2
z
3
z
1
Now the contributions
from the lower limit all vanish and at the upper limit we have ln z =
ln exp(2i) = 2i and z z 1 = 2i sin 2 etc. Therefore the final result is
m() = B0 + B2 sin 2 + B4 sin 4 + B6 sin 6 + ,
(5.77)
= b
n+ n + n ,
B2 = b(1 n)(1 + n)
2
16
2
2
16
2
15 2 15 3
2 15n
B4 = b(1 n)(1 + n)
= b
n + n ,
16
16
16
3
35 3
35n
B6 = b(1 n)(1 + n)2
= b
n .
(5.78)
48
48
5.14
It is straightforward to show that these coefficients are exactly the same as those obtained
in Method I. Simply substitute in the An with expressions for e2 etc. obtained from (5.9).
(Ignoring terms of O(e8 )). Therefore we can write
1
1
mp = m(/2) = A0 = B0
2
2
(5.79)
(5.80)
5.9
When we derived the inverse series for TMS in Chapter 3 we expressed the coefficients
in terms of the footpoint latitude 1 which was defined by msph (1 ) = a1 = y for a
given point (x, y) on the projection. Trivially, 1 = y/a. We will have to do the same for
the ellipsoid but now we are faced with inverting the series (5.71). There are two methods
to choose from. The first is a numerical solution by a fixed point iteration, the second is to
apply the Lagrange method to a series for a function closely related to the meridian distance,
the rectifying latitude.
Inverse meridian distance by numerical methods
To solve m() = y when m() is given by (5.71) or (5.77) we consider the iteration
n+1 = g(n ) = n
(m(n ) y)
,
a
n = 0, 1, 2 . . . ,
(5.81)
where the initial value is that for the spherical approximation: namely 0 = y/a. Now if
the iteration scheme converges so that n+1 and n then (5.81) becomes
(m( ) y)
(5.82)
= g( ) =
a
so that m( ) y = 0 and is the required solution for the footpoint 1 . Note that since
g 0 () 1 B0 /a = O(e2 ) < 1 the iteration will indeed converge.
5.15
m()
m()
m()
=
=
2 mp
A0
B0
(5.83)
where mp = m(/2) is the distance from equator to pole given by (5.71) or (5.77): the
constants A0 and B0 are given by (5.78) or (5.77) respectively. The rectifying latitude may
be used to construct projections from the ellipsoid to the sphere which preserve the meridian
distance but here we use it simply as a parameter which facilitates the series inversion. We
choose to use the expansion in n, equation (5.77).
() =
1
m() = + b2 sin 2 + b4 sin 4 + b6 sin 6 +
B0
(5.84)
where b2 = B2 /B0 , b4 = B4 /B0 etc. The coefficients are to be calculated to O(n3 ) from
equations (5.78):
5
5
set = n + n2 + n3 ,
4
4
B01 = b1 (1 + )1 = b1 (1 + 2 3 + )
1
1
= b1 1 n n2 + n3 +
4
4
3 2 21 3
3
9
1 3
b2 = bB0
n + n + n +
= n + n3 + ,
2
2
16
2
16
15
15
15 2
b4 = bB01
n2 + n3 +
=
n +
16
16
16
35
35
b6 = bB01
n3 +
= n3 + .
(5.85)
48
48
Finally, we invert equation (5.84) by the Lagrange expansion of Appendix B, Section B.5:
= + D2 sin 2 + D4 sin 4 + D6 sin 6 + ,
where, to O(n3 ),
1
D2 = b2 b2 b4 + b32
2
D4 = b4 + b22
3
D6 = b6 + 3b2 b4 b32
2
3
27
n n3 ,
2
32
21 2
=
n ,
16
151 3
=
n .
96
m
,
B0
(5.86)
(5.87)
This completes the derivation of a series for the inverse meridian distance. Note that the
definition of the rectifying latitude is such that it is zero on the equator and /2 at the pole.
From the first term of the above series we see that the maximum of must occur when
45 and it will be of a magnitude given by D2 , approximately 90 .
5.16
5.10
Ellipsoid: summary
p2
Z2
+
= 1.
a2
b2
(5.88)
Parameters
b2 = a2 (1 e2 ),
e02 =
a2 b2
e2
=
,
2
b
1 e2
f=
ab
,
a
e2 = 2f f 2 ,
e1 = n =
(5.89)
ab
.
a+b
(5.90)
Airy ellipsoid
a = 6377563.396m,
e = 0.081673372415,
b = 6356256.910m,
e2 = 0.006670539762,
f = 0.03340850522,
1
= 299.3249753.
f
(5.91)
Cartesian coordinates
X() = p() cos = () cos cos ,
Y () = p() sin = () cos sin ,
Z() = (1 e2 ) () sin .
(5.92)
Coordinate derivatives
dp
= sin ,
d
dZ
= cos .
d
(5.93)
,
[1 e2 sin2 ]1/2
() =
3
2
1
e
,
a2
() =
1 e2 sin2
=
.
1 e2
(5.94)
Curvature derivatives
d
= ( 1) tan ,
d
d
( 1)
=3
tan ,
d
d
= 2( 1) tan .
d
(5.95)
Metric
ds2 = 2 d2 + 2 cos2 d2 .
(5.96)
/continued
5.17
Meridian distance
m() = A0 + A2 sin 2 + A4 sin 4 + A6 sin 6 + ,
(5.97)
(5.98)
Rectifying latitude
() =
m()
.
2 mp
(5.99)
(5.100)
In the above series the An , Bn and Dn coefficients are given by (5.72), (5.78) and (5.87).
5.18
Chapter
6.1
(6.1)
y(, ) = a f ().
(6.2)
The geometry of the projection may be illustrated by Figure 2.4 with only one change,
the normal at P no longer passes through the centre in general. The meridians on the
ellipsoid are projected into lines parallel to the y-axis on the projection and parallel circles
are projected into lines parallel to the x-axis. The meridian spacing is equal but the spacing
6.2
of the projected parallels will of course depend on the nature of the function f (). Note
that the orthogonal intersections of meridians and parallels on the graticule are transformed
into orthogonal intersections on the map but this is not necessarily true for intersections at
an arbitrary angle.
The essential difference is that the metric on the ellipsoid is given by (5.62):
ds2 = 2 d2 + 2 cos2 d2 ,
(6.3)
where () and () are defined by equations (5.16) (5.45) respectively. The corresponding
infinitesimal elements on the ellipsoid and the plane of projection are shown in Figure 6.1.
Figure 6.1
tan =
cos
,
and
so that
tan =
(b)
tan =
x
a
=
,
y
a f 0 ()
sec
tan .
f 0 ()
(6.4)
(6.5)
() sec
.
()
(6.6)
QP
s02
x2 + y 2
=
lim
.
QP 2 2 + 2 cos2 2
s2
(6.7)
On P Q and P 0 Q0 equations (6.4) give = (/) cot cos and y= cot x. Therefore we have
x2 (1 + cot2 )
2 = lim 2
.
(6.8)
QP cos2 2 (cot2 + 1)
6.3
a sec sin
,
() sin
(6.9)
where we assume that has been found in terms of and from equation (6.5). Clearly if
the projection is conformal, with = , we have an isotropic scale factor with
k() =
a sec
.
()
(6.10)
This scale factor differs from that on the sphere by the factor of a/, a difference of O(e2 ).
There is no simple expression for k as a function of y. Given y we must first use one of the
methods of inversion discussed in Section 6.3 to find from y and then we apply (6.10).
6.2
The Mercator parameter (or isometric latitude) for NME will be be denoted by () so that
the equations of the projection are still written as
x(, ) = a,
y(, ) = a().
(6.11)
Warning. We use the same notation for the Mercator parameter on both the sphere and the ellipsoid although they are
of course different functions. From this point will always
denote the ellipsoidal form which is derived below.
The condition that NME be conformal follows from equation (6.6) when f () ().
d
() sec
=
.
d
()
(6.12)
Substituting for and from equations (5.16) and (5.45), splitting into partial fractions and
noting that the first term simply gives the same integral as (2.26), we have
Z
1
(1 e2 )
() =
d
2 sin2
cos
e
0
Z
=
0
1
e2 cos
cos
2
1
1
+
1 + e sin 1 e sin
e
1 + e sin
= ln tan
+
ln
.
2
4
2
1 e sin
d
(6.13)
The first term is just the NMS Mercator parameter for the sphere so we see that the parameters on sphere and ellipse differ by terms of O(e2 ). Note that still diverges to at
= /2. Obviously we recover the parameter for the spherical case when e 0.
6.4
Alternative forms
As in Section 2.4 we can rewrite the Mercator parameter for the ellipsoid in many different
forms. The most useful are (a) a simple rearrangement of (6.13), (b) likewise but with
tan(/2 + /4) replaced by cot(/4 /2) and (c) a replacement of the tangent term as
in equation (2.31):
"
() =
ln tan
"
() = ln tan
1
() = ln
2
6.3
+
2
4
4
2
1 + sin
1 sin
1 e sin
1 + e sin
e/2 #
1 + e sin
1 e sin
e/2 #
1 e sin
1 + e sin
(6.14)
(6.15)
e
.
(6.16)
Inverting the equations x = a and y = a to find and is trivial but finding a value of
from = y/a is anything but trivial. There is no way in which we can manipulate any
of the expressions for () to give () in a closed form. We have to resort to one of the
following methods.
Construct an iterative scheme from which we can find a numerical value for for any
given value of = y/a.
Use a Taylor series expansion of (). In particular we shall need an expansion about
the footpoint parameter 1 when we come to the transverse Mercator on the ellipse
(TME) in the next chapter.
It is clear that the first method is problematic. The range of is finite, [/2, /2],
but we have as /2; there is obviously no hope of obtaining a series
of the form = + terms of order e, e2 , . . . which would be amenable to inversion by
a Lagrange expansion. However, in Sections 6.46.6, we shall see that this is possible for
another parameter which is closely related to . In the remainder of this section we consider
the numerical solution and the Taylor series.
6.5
1 e sin e/2
= 2 arctan exp()
,
(6.17)
2
1 + e sin
from which we construct an iteration
n+1
"
#
1 e sin n e/2
= 2 arctan exp()
,
2
1 + e sin n
(6.18)
where n = 0, 1, 2 . . . For the initial value 0 we use the spherical approximation. Setting
e = 0 in equation (6.17) we have
0 = 2 arctan[exp()].
(6.19)
2
The choice of (6.15) rather than (6.14) introduces the factor of exp() which clearly
facilitates the convergence when is large and positive. If the iteration does converge to a
value , then n and n+1 in (6.18) so that
"
#
e/2
e
sin
= 2 arctan exp()
.
(6.20)
2
1 + e sin
Thus is the required solution of (6.17).
(6.22)
c = cos
t = tan
= () = / .
(6.23)
6.6
Note that an expression for d/d evaluated at the footpoint latitude 1 is exactly what we
shall need in the application of this series.
We now construct expressions for all the derivatives in the Taylor series as functions
of . We need the derivative of () given in equation (5.54) as 0 = (2 2)t where
t = tan (and dt/d = 1 + t2 ).
d
= c,
d
d
d
d 0
d2
=
[c] =
[c]
= c s (c)
2
d
d
d
d
= [(2 2)tc s] (c) = c2 t 3 2 + 2 ,
d
d 2
d3
=
c t 3 2 + 2
3
d
d
d
= (c) 2cst + c2 (1 + t2 ) 3 2 + 2 + c2 t [6 + 2] (2 2)t
= c3 3 (3 + 15t2 ) + 2 (2 18t2 ) + (4t2 ) ,
d
d 3 3
d4
=
c (3 + 15t2 ) + 2 (2 18t2 ) + (4t2 )
4
d
d
d
= c4 t 4 (57105t2 ) + 3 (68+180t2 ) + 2 (1684t2 ) + (8t2 )
(6.24)
These derivatives must be evaluated at 1 and substituted into the Taylor series which we
now write as
1 = (1 ) 1 c1 +
(1 )2
(1 )3
(1 )4
1 c21 t1 D2 +
1 c31 D3 +
1 c41 t1 D4 ,
2!
3!
4!
(6.25)
D4 =
(6.26)
5 t21 .
(6.27)
6.7
6.4
Projections can be defined from the ellipsoid to any surface of reasonable shape, not just to
the plane. Here we shall consider only the case of projections from the ellipsoid to a sphere
of radius R. Such a projection can then be followed by one of the many known projections
from the sphere to the plane to produce a double projection with desirable properties.
If we denote the latitude and longitude coordinate on the sphere by (, ) then the
projection to the sphere is defined by two (well-behaved) functions, (, ) and (, ),
where and are the usual geodetic longitude and latitude on the ellipsoid. To complete
the double projection we define coordinates (x, y) on the plane by specifying a further two
functions x(, ) and y(, ).
The basic properties of such a projection from the ellipsoid to the sphere can be investigated by comparing the infinitesimal elements shown in Figure 6.2.
Figure 6.2
tan =
cos
and
(b)
tan 0 =
R cos
.
R
(6.28)
cos
tan .
0 () cos
(6.29)
It is straightforward to calculate the scale factor for any azimuth but we shall consider only
the scale factor on the meridian of the sphere; denoting this by h(), as in Section 2.2, we
6.8
have
h() =
R
R0 ()
=
.
(6.30)
Note that we have not specified how the radius R is to be chosen. There are many possible
choices which, whilst not affecting the angle transformations, will of course influence the
scale factor. All of the following have been used for R.
The meridian radius of curvature, , at a latitude where we seek the best fit.
The Gaussian radius of curvature, , at a latitude where we seek the best fit.
R such that the ellipsoid and sphere have the same volume.
R such that the ellipsoid and sphere have the same surface area.
R such that the ellipsoid and sphere have the same equatorpole distance.
sec 0 () =
which integrates to
Z
()
Z
sec d =
(6.31)
() sec
d.
()
(6.32)
The integral on the left is the same as that for the Mercator parameter on the sphere, equation (2.26), whilst the integral on the right is that which gives the Mercator parameter on
the ellipse, equation (6.13). Therefore
"
#
()
1 e sin e/2
ln tan
+
= ln tan
+
(),
(6.33)
2
4
2
4
1 + e sin
"
() = 2 arctan tan
+
2
4
1 e sin
1 + e sin
e/2 #
.
2
(6.34)
6.9
This rather complicated transformation (along with = ) has been constructed to guarantee a conformal projection from ellipsoid to sphere. But note that equations (6.30) and (6.31)
clearly show that the scale cannot be uniform on any meridian of the sphere. Therefore following this projection with TMS from sphere to plane will produce a conformal projection
of the ellipsoid to the plane but with a non-uniform scale on the central meridian. When we
meet TME (next Chapter) we shall find that it is a conformal projection of the ellipsoid to
the plane with a uniform scale on the central meridian.
A rectifying projection from ellipse to the sphere
As a second example of a projection to the sphere consider that defined by setting =(),
where () is the rectifying latitude defined in Section 5.9, equation (5.83):
() = () =
m()
2 mp
(6.35)
where mp = m(/2) is the meridian distance between equator and pole on the ellipsoid.
The corresponding distance on the sphere is (/2)R and the two are clearly equal if we set
R = 2mp /. In terms of the notation developed in Section 5.8 we have R = A0 = B0 so
that we have a > R > b as expected.
The scale factor on the meridian is then
h=
R0 ()
R m0 ()
=
= 1,
2 mp
(6.36)
since m0 () = from equation (5.68). Thus this projection to the sphere conserves the
scale factor and total length on every meridian but on the other hand it is not a conformal
projection since () and () are different functions.
6.5
Equation (6.34) shows that and are equal at the equator and at the pole and elsewhere
they differ by O(e) terms, (actually by O(e2 ) as we shall see), so there is every reason to
expect them to be related by a series which can be inverted by the Lagrange method. In this
section we will find coefficients such that
() = + b2 sin 2 + b4 sin 4 + b6 sin 6 + b8 sin 8 + ,
(6.37)
(6.38)
To develop the direct series for () we first introduce two new (small) parameters, A
and , such that
1 e sin
1 + e sin
e/2
= exp A = 1 + .
(6.39)
6.10
Taking logarithms
A =
e
ln
2
1 e sin
1 + e sin
e
= ln
2
1 es
1 + es
.
Using the series (E.12) we obtain A and its powers to order O(e8 ):
1 2 2 1 4 4 1 6 6
2
A = e s 1 + e s + e s + e s ,
3
5
7
2
23
A2 = e4 s2 1 + e2 s2 + e4 s4 ,
3
45
A3 = e6 s3 1 + e2 s2 ,
A4 =
Therefore
e8 s4 .
1 2
1
1
A + A3 + A4 +
2!
3!
4!
= p2 e2 + p4 e4 + p6 e6 + p8 e8 + ,
= exp A 1 = A +
(6.40)
p2 = s
p4 =
(6.41)
(6.42)
+
= cos = .
2
2
1+b
sec (/2 + /4)
2
2
2
2
With this notation equation (6.34) becomes
+ = arctan[b(1 + )] = arctan(b + b).
2
4
For the inverse tangent we use the series (E.9) with z = b. Therefore
(6.43)
(6.44)
b 1
(b)2
2b
+ = + +
(6.45)
2
4
2
4
1! 1 + b2
2! (1 + b2 )2
(b)3
2
8b2
(b)4
24b
48b3
+
+
+ ,
+
3!
(1 + b2 )2 (1 + b2 )3
4!
(1 + b2 )3 (1 + b2 )4
6.11
Substitute first for the powers of [b/(1 + b2 )] from equation (6.43) and then substitute for
the remaining powers of b from equation (6.42). The result is
c
c
c
= + c (1 + s) 2 + (1 + 3s + 2s2 ) 3 (s + 2s2 + s3 ) 4 + .
2
6
4
(6.46)
Substitute for the powers of in terms of the p-coefficients (equation 6.40) and then substitute for the p-coefficients using equation 6.41. This gives a series of the form
= + q2 e2 + q4 e4 + q6 e6 + q8 e8 + ,
(6.47)
,
2
24
32
5760
5e4 7e6
697e8
b4 =
+
+
+ ,
48
80
11520
13e6
461e8
,
b6 =
480
13440
1237e8
b8 =
+ .
161280
(6.49)
b2 =
(6.50)
(6.51)
6.12
where
1
d2 = b2 b2 b4 + b32
2
4
d4 = b4 + b22 2b2 b6 + 4b22 b4 b42
3
3 3
d6 = b6 + 3b2 b4 b2
2
8
d8 = b8 + 2b24 + 4b2 b6 8b22 b4 + b42
3
6.6
e2 5e4
e6
13e8
+
+
+
+ ,
2
24
12
360
7e4 29e6
811e8
=
+
+
+ ,
48
240
11520
7e6
81e8
=
+
+ ,
120 1120
4279e8
=
+ .
(6.52)
161280
=
Returning to equation (6.33) we see that the relation between the Mercator parameter and
the conformal latitude is
h
i
h
i
() = ln tan
+
= ln cot
.
(6.53)
2
4
4
2
Inverting the second of these gives
() =
2 arctan [exp()] .
2
(6.54)
This result, alongwith the last section, provides another inversion of the Mercator parameter.
Given = y/a we use equation (6.54) to calculate and then use the series (6.51) to find .
This solves the problem of the inverse transformation and also allows us to find the scale
factor for any y using (6.10).
6.7
NME can be modified exactly as NMS to provide a slightly wider domain near the equator
in which the scale is accurate to within a given tolerance. For a tolerance of 1 in 2500 the
range of latitude will differ from that for NMS (Section 2.7) by terms of order e2 . The
technique is the same, we simply introduce a factor of k0 into the transformations.
Direct transformation
"
x = k0 a,
y = k0 a(),
() = ln tan
+
2
4
1 e sin
1 + e sin
e/2 #
(6.55)
6.13
Inverse transformation
= x/k0 a,
= () with = y/k0 a
(6.56)
1 e sin n e/2
n+1 = 2 arctan exp()
,
2
1 + e sin n
(6.57)
2 arctan [exp().]
2
(6.58)
where the coefficients are given in equation (6.52) and (the conformal latitude) is defined
in terms of by
() =
2 arctan [exp()] .
2
(6.59)
k0 a sec
()
(6.60)
Scale factor
k() =
To find the scale for a given y on the projection we use the above results to first find when
given = y/k0 a.
6.14
Chapter
TME is derived as a series by a complex transformation from the NME projection. The method parallels that used in Chapter 4 for the derivation of the TMS
series from NMS.
Figure 7.1
7.1
Introduction
In Chapter 6 we derived the NME projection: it can be considered as a conformal transformation from a point P (, ) on the ellipsoid to a point on the complex -plane defined by
= + i where () is the Mercator parameter for the ellipsoid given in (6.14).
Let (x, y) be the coordinates of the required TME projection and let z = x + iy be a
general point on the associated complex plane. The aim of this chapter is to find a conformal
transformation
z() x(, ) + iy(, ),
(7.1)
such that (a) the central meridians, = 0 and x = 0, map into each other, and (b) the scale
is true on the y-axis. Our method parallels that of Chapter 4 with the appropriate definitions
of the Mercator parameter and meridian distance for the ellipsoid. The resulting series are
those first given by Lee (to sixth order) and Redfearn (to eighth order)see bibliography.
7.2
(7.2)
where the A-coefficients are given in equations (5.72). In considering the transformation
from the complex -plane to the complex z-plane it is useful to express the meridian distance
as a function of and write it as M (), where
M () = m().
(7.3)
There is no closed expression for M () analogous to (4.6); this is of no import since we
only need its derivatives. (See next page).
Footpoint latitude and parameter
Given a point P 00 with projection coordinates (x, y) then the projection coordinates of the
footpoint are (0, y). The definition of the footpoint latitude 1 and the footpoint parameter 1 are unchanged from those of Sections 3.3 and 4.1: they are the solutions of
m(1 ) = y,
M (1 ) = y.
(7.4)
We shall need to calculate the footpoint latitude (but not the footpoint parameter) for a
given y. One method of finding the solution of m()=y is to use the fixed point iteration
given in equation (5.81),
(m n ) y
n+1 = g(n ) = n
,
n = 0, 1, 2 . . . ,
(7.5)
a
starting with the spherical approximation 0 = y/a.
Alternatively, we can use the series given in (5.86) with m() = y,
= + d2 sin 2 + d4 sin 4 + d6 sin 6 + ,
y
,
B0
(7.6)
(7.7)
We shall not need this explicit form but we shall require require its derivatives. From (6.12)
d
()
=
,
d
() cos
d
() cos
=
.
d
()
(7.8)
7.3
We shall also need (), the inverse of the Mercator parameter, as a fourth order Taylor
series about the footpoint parameter 1 . This is given in equation (6.25).
1 = (1 ) 1 c1 +
(1 )2
(1 )3
(1 )4
1 c21 t1 D2 +
1 c31 D3 +
1 c41 t1 D4
2!
3!
4!
(7.9)
where the D-coefficients are given in (6.26). The 1 suffix of course denotes a term calculated at the footpoint latitude.
The derivatives of the meridian distance
We shall need the derivative of the M () as functions of . From (5.68) we have
dm()
= (),
d
(7.10)
dM (()) d
dm() cos
dM ()
=
=
= () cos .
d
d
d
d
(7.11)
Proceeding in this way we can construct all the derivatives of M () with respect to but
with the results expressed as functions of . Denoting the n-th derivative of M with respect
to by M (n) the exact results for the first six derivatives are given below. We use the usual
compact notation for sin etc. and also make frequent use of the derivatives of () and
() given in equation (5.54):
d
= ( 1) tan ,
d
d
= 2( 1) tan .
d
M (1) =
dM
=
d
c.
M (2) =
d2 M
=
d 2
d (1) d
d
d
c
M
=
(c)
= [( 1)tc s]
d
d
d
d
= sc.
M (3) =
c
d3 M
= ( 1)tsc + (c2 s2 )
= c3 t2
3
d
(7.12)
(7.13)
(7.14)
(7.15)
c3 W3 .
/cont.
7.4
M (4) =
h
i c
d4 M
3
2
2
3
2
=
(1)tc
3c
s
(t
)
+
c
2(1)t
2t(1+t
)
d 4
= sc3 4 2 + t2
(7.16)
M (5) =
d5 M
=
d 5
=
1ex
M
(6)
d6 M
=
=
d 6
sc3 W4 .
h
( 1)tsc3 + (c4 3s2 c2 ) (4 2 + t2 )
i c
+ sc3 (8 + 1)(2 + 2)t 2t(1 + t2 )
c5 4 3 (1 6t2 ) + 2 (1 + 8t2 ) 2t2 + t4
(7.17)
c W5 .
i
5
0 c
( 1)tc + (5sc ) W5 + c W5
h
= sc5 (4 1) 4 3 (1 6t2 ) + 2 (1 + 8t2 ) 2t2 + t4
+ t1 12 2 (1 6t2 ) + 2(1 + 8t2 ) 2t2 (2 + 2) t
i
24 3 8 2 + 2 2t(1 + t2 ) + 4t3 (1 + t2 )
= sc5 8 4 (11 24t2 ) 28 3 (1 6t2 ) + 2 (1 32t2 ) 2t2 + t4
h
sc5 W6
(7.18)
We shall find that the derivatives M (7) and M (8) multiply 7 and 8 terms respectively
and we shall later justify the neglect of terms of order e2 7 and e2 8 . Accordingly, we
evaluate these derivatives in the spherical limit in which e 0 and 1, (except that the
overall multiplicative factors of are not set equal to a for the sake of visual conformity
with the lower order derivatives, not to improve accuracy). Noting that 0 = 0 = 0 in this
limit we find
h
i c
d7 M
6
2 4
5
0
M (7) =
=
5s
c
W
|
+
sc
W
6
6 =1
=1
d 7
7
2
2
4
= c 1 5t
61 58t + t + t 116t + 4t3 (1 + t2 )
= c7 61 479t2 + 179t4 t6 .
c7 W7 .
M (8) =
d8 M
=
d 8
(7.19)
i
0 c
7c6 sW7 c7 W7
1
7
2
4
6
3
5
2
sc 7 61479t +179t t 958t+716t 6t (1+t )
t
7
2
4
6
sc 1385 3111t + 543t t .
sc7 W8 .
(7.20)
7.5
Summary of derivatives
M (1) =
M (2) = sc
M (3) = c3 W3
W3 () = t2
M (4) =
sc3 W4
W4 () = 4 2 + t2
M (5) =
c5 W5
M (6) = sc5 W6
M (7) = c7 W7
M (8) =
sc7 W8
(7.21)
The bar on W7 and W8 denotes that the term is evaluated in the spherical limit. This notation
will be standard from here on. Later we will need the expressions for W3 , . . . W6 in the
spherical approximation: setting = 1 gives
W3 () W3 () = 1 t2 ,
W4 () W4 () = 5 t2 ,
W5 () W5 () = 5 18t2 + t4 ,
W6 () W6 () = 61 58t2 + t4 .
7.2
(7.22)
z = z0 +(0 )M0
(n)
(7.24)
7.6
For the direct series we start from a given (arbitrary) point P 0 at = + i and choose
0 = i with the same ordinate in the -plane. Therefore in the Taylor series (7.23) we set
=+
()
Figure 7.2
(7.25)
The real and imaginary parts of equation (7.25) give x and y as functions of and :
1
1
1 3 (3)
M + 5 M (5) 7 M (7) +
3!
5!
7!
1
1
1
1
y(, ) = M 2 M (2) + 4 M (4) 6 M (6) + 8 M (8) + .
2!
4!
6!
8!
x(, ) = M (1)
(7.26)
(7.27)
Writing M and its derivatives as functions of from (7.3) and (7.21) gives the Redfearn
formulae for the direct transformation as power series in (radians):
x(, ) = c +
3 c3
5 c5
7 c7
W3 +
W5 +
W7 ,
3!
5!
7!
y(, ) = m() +
2 sc 4 sc3
6 sc5
8 sc7
+
W4 +
W6 +
W8 .
2
4!
6!
8!
(7.28)
(7.29)
Since all the coefficients on the right hand sides are now expressed in terms of , we have
replaced x(, ) and y(, ) on the left hand side by x(, ) and y(, ) respectively.
Conformality and the CauchyRiemann equations
The conformality of the above transformations may be confirmed by evaluating the Cauchy
Riemann equations (4.14)
1 2 (3)
1
1
M + 4 M (5) 6 M (7) + ,
2!
4!
6!
1
1
1
x = y = M (2) 3 M (4) + 5 M (6) 7 M (8) + .
3!
5!
7!
x = y = M (1)
(7.30)
(7.31)
7.7
z z0
b2
b3
b8
= (0 ) + (0 )2 + (0 )3 + + (0 )8 +
0 c0
2!
3!
8!
(7.32)
b2 =
iM0
0 c0
b3 =
M0
0 c0
b4 =
iM0
0 c0
b5 =
M0
0 c0
b6 =
iM0
0 c0
b7 =
M0
0 c0
b8 =
iM0
0 c0
is0
(3)
c20 W3 (0 )
(4)
= is0 c20 W4 (0 )
(5)
c40 W5 (0 )
(6)
= is0 c40 W6 (0 )
(7)
c60 W7 (0 )
(8)
= is0 c60 W8 (0 )
(7.33)
where the functions on the right hand sides are evaluated at 0 such that 0 = (0 ).
The Lagrange inversion of an eighth order series is developed in Appendix B, Sections B.6B.8. If we identify the series (7.32) with (B.23) by replacing (z z0 )/0 c0 and
( 0 ) by w and z respectively we can use (B.24) to deduce that the inverse of (7.32) is
p 2 z z0 2 p 3 z z0 3
p 8 z z0 8
z z0
0 =
,
(7.34)
0 c0
2! 0 c0
3! 0 c0
8! 0 c0
where the p-coefficients are given by equations (B.25) and (B.30). We shall actually need
these coefficients at the footpoint latitude 1 and we choose to write them as
p2 = ic1 t1 ,
p3 =
c21 V3
p4 = ic31 t1 V4
p5 =
c41 V5
p6 = ic51 t1 V6
p7 =
c61 V7
p8 = ic71 t1 V8
V3 = 1 + 2t21 ,
V4 = 412 91 6t21 ,
V5 = 413 (1 6t21 ) 12 (9 68t21 ) 721 t21 24t41 ,
V6 = 814 (1124t21 )8413 (38t21 )+22512 (14t21 )+6001 t21 +120t41 ,
V7 = 61 + 662t21 + 1320t41 + 720t61 ,
V8 = 1385 7266t21 10920t41 5040t61 .
(7.35)
7.8
=+
Figure 7.3
For the inverse we start from an arbitrary point with projection coordinates P 00 (x, y) and
move to the footpoint at K 00 (0, y) so that we set zz0 = (x+iy)iy = x in equation (7.34).
We must then set 0 = i1 where 1 is the footpoint parameter such that M (1 ) = y. Let
1 be the corresponding footpoint latitude such that m(1 ) = y. Therefore (7.34) becomes
x
p2 x 2 p3 x 3
p8 x 8
+ i i1 =
,
(7.36)
1 c1
2! 1 c1
3! 1 c1
8! 1 c1
where the p coefficients at the footpoint latitude 1 have already been given in (7.35). Taking the real and imaginary parts we have
(x, y) =
x3
x5
x7
x
V7 ,
3
5
1 c1 3! 13 c1
5! 15 c1
7! 17 c1
1 =
where m(1 ) = y,
x2 t1
x4 t1
x6 t1
x8 t1
V8 .
4
6
2! 12 c1 4! 14 c1
6! 16 c1
8! 18 c1
(7.37)
(7.38)
We see that the spherical approximation has been used only in the last term of each series.
Therefore it is equivalent to neglecting terms of order e2 (x/a)7 and e2 (x/a)8 . This will be
justified when we look at the typical magnitude of such terms.
Finally, we note for future reference the spherical limits of the terms V3 , . . . V6 : since
1 1 as e 1 we have
V3 V 3 =
1 + 2t21 ,
V4 V4 = 5 6t21 ,
V5 V5 = 5 28t21 24t41 ,
V6 V 6 =
61 + 180t21 + 120t41 .
(7.39)
(1 )2
(1 )3
(1 )4
1 c21 t1 D2 +
1 c31 D3 +
1 c41 t1 D4 ,
2!
3!
4!
(7.40)
7.9
where the D-coefficients are given in (6.26). All that remains is to substitute for ( 1 )
using (7.38). It is convenient to use a temporary abbreviation, setting x
e = x/1 .
t1 1 2
1
1
1
x
e + V4 x
e 4 + V6 x
e 6 + V8 x
e8 ,
c1 2!
4!
6!
8!
t21 1 4
2
2
1 2 8
6
8
( 1 )2 =
,
x
e
+
V
x
e
+
V
x
e
+
V
x
e
4
6
2!4!
2!6!
4!4! 4
c21 2!2!
t3
1
3
( 1 )3 = 13
x
e6 +
V4 x
e8 ,
2!2!4!
c1 2!2!2!
1
x
t41
8
where x
e= .
( 1 )4 =
x
e
1
c41 2!2!2!2!
1
(7.41)
1 =
where we use the spherical approximation in evaluating the eighth order term. Substituting for the D-coefficients from equations (6.26, 6.27) and the V -coefficients from equations (7.35,7.39) our final result for is
(x, y) =
x2 1 t1 x4 1 t1
x6 1 t1
x8 1 t1
U 8,
4
6
212
4!14
6!16
8!18
(7.43)
where
U4 = 412 91 (1 t21 ) 12t21 ,
U6 = 814 (11 24t21 ) 1213 (21 71t21 ) + 1512 (15 98t21 + 15t41 )
+ 1801 (5t21 3t41 ) + 360t41
U 8 = 1385 3633t21 4095t41 1575t61 .
(7.44)
61 + 90t21 + 45t41 ,
(7.45)
7.10
7.3
Direct
x(, ) = c +
3 c3
5 c5
7 c7
W3 +
W5 +
W7
3!
5!
7!
y(, ) = m() +
Inverse
(x, y) =
(7.46)
2 sc 4 sc3
6 sc5
8 sc7
+
W4 +
W6 +
W8 (7.47)
2
4!
6!
8!
x5
x7
x3
x
V
V7
3
5
1 c1 3! 13 c1
5! 15 c1
7! 17 c1
(7.48)
(x, y) = 1
x2 1 t1 x4 1 t1
x6 1 t1
x8 1 t1
U8
4
6
212
4!14
6!16
8!18
(7.49)
(x, y) = 1
x4 t1
x6 t1
x8 t1
x2 t1
V8
4
6
2! 12 c1 4! 14 c1
6! 16 c1
8! 18 c1
(7.50)
where
W 3 = t2
W4 = 4 2 + t2
W5 = 4 3 (1 6t2 ) + 2 (1 + 8t2 ) 2t2 + t4
W6 = 8 4 (1124t2 ) 28 3 (16t2 ) + 2 (132t2 ) 2t2 + t4
W7 = 61 479t2 + 179t4 t6 + O(e2 )
W8 = 1385 3111t2 + 543t4 t6 + O(e2 )
(7.51)
V3 = 1 + 2t21
V4 = 412 91 6t21
V5 = 413 (1 6t21 ) 12 (9 68t21 ) 721 t21 24t41
V6 = 814 (1124t21 )8413 (38t21 )+22512 (14t21 )+6001 t21 +120t41
V7 = 61 + 662t21 + 1320t41 + 720t61
V8 = 1385 7266t21 10920t41 5040t61
(7.52)
(7.53)
and (a) m() is the meridian distance which may be calculated from the series (5.71)
or (5.77); (b) s, c, t denote sin , cos , tan ; (c) the functions () and () are defined in equation (5.52, 5.53); (d) is measured in radians from the central meridian; (e) the
7.11
subscripted terms in the inverse series are to be evaluated at the footpoint latitude 1 such
that m(1 ) = y.
Comments
1. The series given on the previous page are in full accordance with those printed in Redfearns article in the (Empire) Survey Review. They differ in format since
Redfearn
writes the series in terms of () and () rather than () = / and (). In
the original paper Redfearns formulae a few 1 subscripts are omitted in the inverse
series.
2. We shall defer an analysis of the accuracy of the above series until we present the
results for scale and convergence in the next chapter. Then in Chapter 9 we discuss
the results in the context of two important applications, UTM and NGGB.
3. The TME projection may be modified so that the scale on the central meridian is less
than unity, equal to or approximately equal to 0.9996 for UTM and NGGB respectively. Unit scale is achieved on two lines but for TME these are neither meridians
nor grid lines. This modification reduces the range of scale variation.
4. For computational purposes the series should be written in a nested form. For example equation (7.46) can be written as
!)#
"
(
W
W
W
5
3
e2
e2 7
e 1+
e2
e = c.
+
+
x =
,
(7.54)
3!
5!
7!
7.12
Chapter
8.1
x = y ,
x = y .
(8.1)
x = y ,
y = x .
(8.2)
It is straightforward, and useful, to show that the CauchyRiemann equations for the inverse
transformation follow from (8.1). Consider the identities
x = x (x, y), (x, y) ,
y = y (x, y), (x, y) .
(8.3)
Differentiate both these identities by x and then both by y, giving four equations:
1 = x x + x x ,
0 = x y + x y ,
0 = y x + y x ,
1 = y y + y y .
(8.4)
Eliminating x then x from the left pair and then y and y from the right pair:
y = J x ,
y = J x ,
x = J y ,
x = J y ,
(8.5)
where
J = x y y x = x2 + y2 = |x + iy |2 = |z 0 ()|2
=
1
1
1
=
= 2
.
| 0 (z)|2
|x + ix |2
x + x2
(8.6)
8.2
Now neither z 0 () nor 0 (z) can vanish, for otherwise the transformation and its inverse
would degenerate into z=constant or =constant respectively (so that every P (, ) would
map to the same point). Therefore J must be non-zero and bounded and equations (8.5) can
be used to show that the equations (8.2) follow from (8.1).
8.2
Azimuths and grid bearings in TME
Following Section 3.5 we consider the geometry of infinitesimal elements on the ellipsoid
and corresponding elements on the NME and TME complex planes. Both steps of the
projection are conformal and the azimuth on the ellipsoid is equal to the angle between
projected meridian and the projected displacement on both NME and TME planes. The
(b)
(a)
(c)
Figure 8.1
grid bearing on the TME projection is defined as the angle between grid north (on which
x = 0) and the tangent to the projected displacement P 00 Q00 . The TME projection convergence, , is defined as the angle between TME grid north and the projected meridianand
because of conformality we see that the displacents P 0 M 0 , P 0 Q0 and P 0 N 0 are all rotated
by to the corresponding displacements P 00 M 00 , P 00 Q00 and P 00 N 00 . With these definitions
we see that just as in Section 3.5 we have
=+
(8.7)
This equation can be used to find bearings from azimuths or vice versa if we know (, )
and (x, y) respectively. That is =+(x, y) or =(, ).
8.3
We have defined grid (or projection) convergence, , as the angle between true north and
grid north on the projection at P 00 , where true north is defined by the tangent to the projected
meridian at P 00 . Since = 0 on the projected meridian, and remembering that we are now
considering x and y as functions of and , we have
x + x
x
dx
y
tan =
=
lim
=
=
,
(8.8)
dy =0
y + y =0
y
x
where we use the Cauchy-Riemann conditions (8.1) in the last step. This last form is to be
8.3
preferred since the Redfearn formulae (7.46, 7.47) give x and y as power series in with
coefficients as functions of . Therefore we have
y
(, ) = arctan
.
(8.9)
x
The calculation of the convergence in projection coordinates follows immediately if we use
equations (8.5) to set y = J x and x = J y :
x
(x, y) = arctan
,
(8.10)
x
where and are power series in x (7.50, 7.48) with coefficients evaluated at the footpoint 1 .
8.4
If we work from first principles (compare equation 2.19) the square of the point scale factor
is given by the ratio of the metric distances on the TME projection and the ellipsoid:
2 =
x2 + y 2
.
2 2 + 2 cos2 2
(8.11)
From the derivative of the Mercator parameter on the ellipsoid, equation (6.12), we have
that = cos . Therefore we can write the above as
2 =
a2
x2 + y 2
2
kNME
m2
2 cos2 a2 ( 2 + 2 )
(8.12)
The first of these factors is just the square of the scale for the transformation from the ellipsoid to NMEsee equation (6.10). The second factor, which we have defined as m2 ,
is the square of the scale factor between the NME and TME projections (or the magnification of the conformal transformation between the NME, TME complex planes). Using the
CauchyRiemann equations (8.1) we can rewrite the numerator as follows:
x2 + y 2 = (x + x )2 + (y + y )2
=
=
(x2 + y2 ) 2 + 2(x x
(x2 + y2 )( 2 + 2 ).
+ y y ) +
(8.13)
(x2
y2 )2
(8.14)
(8.15)
where we use (8.9) in the second step. We have also replaced by k, the usual notation for
an isotropic scale factor, since x and y are independent of (depending only on position).
Using equations (8.6, 8.10) this result can be transformed to
2
1/2
1
1
1
k(x, y) =
x + x2
=
,
(8.16)
() cos
() cos x sec (x, y)
where it is assumed that we have first calculated (x, y).
8.4
r
Therefore
y
= arctan
,
m = {x2 + y2 }1/2 ,
(8.18)
x
reproducing the previous definitions in equations (8.9) and (8.15), the latter after multiplication by the scale factor kNME for the transformation from the ellipsoid to NME and after
dividing by a factor of a to allow for the fact that we chose (Chapter 4) to label the NME
complex plane by (, ) rather than (a, a).
8.5
The expressions for the convergence and scale factors in equations (8.9), (8.10), (8.15) and
(8.16) depend on the partial derivatives x , y , x , x of the Redfearn series (7.467.50).
e = c and x
Setting
e = x/1 we have
x
1 e2
1 e4
1 e6
e = c
x =
(8.19)
=
c 1 + W3 + W5 +
W7 , ,
2
24
720
1 e2
1 e4
1 e6
y
e
y =
= s 1 + W4 +
W6 +
W8 ,
(8.20)
6
120
5040
1
1 2
1 4
1 6
x
x =
=
1 x
e V3 x
e V5
x
e V7 ,
x
e=
(8.21)
x
1 c1
2
24
720
1
t1 x
e
1 2
1 4
1
6
x =
=
e V4 +
x
e V6 +
x
e V8 .
(8.22)
1+ x
x
1 c1
6
120
5040
Note that, apart from the overall multiplicative terms, the series for x is obtained from that
ex
for x by the replacements
e and Wn Vn (n odd); the series for x is obtained
ex
from that for y by the replacements
e and Wn Vn (n even).
NB. In constructing the Redfearn series we discarded terms of order 9 and (x/a)9 so we
must discard terms of order 8 and (x/a)8 (and higher order) in the derivatives and in any
e7 ,
expressions obtained by manipulation of the above series. Moreover the coefficients of
7
8
8
e
x
e , and x
e terms of the Redfearn series were evaluated in the spherical approximation
(e = 0, = 1); therefore, for consistency, we must use the spherical approximation in
e6 , x
e7 and x
coefficients of terms of the order
e 6,
e 7 wherever they arise in the manipulation
of the series for the derivatives.
8.5
.
=1
+
x
2
24
4
720
24
8
(8.23)
Note that the W3 and W5 which the inversion casts into the last term have been replaced by
their spherical limits. The product of the above with (8.20) gives
h
i
y
e 1 + a2
e2 + a4
e4 + a6
e6 ,
= t
(8.24)
x
where the coefficients (and required spherical limits) are calculated using (7.21, 7.22)
i
1h
1
W4 W3
1 + t2
= 2 2 + t2 ,
a2 =
a2 =
6
2
3
3
W6
W3 W4 W32 W5
,
120
12
4
24
i
1h 4
=
(1124t2 ) 3 (1136t2 )+ 2 (24t2 )4t2 +2t4 ,
15
a4 =
a4 =
1
2+4t2 +2t4
15
W8 W3 W6 W3 W4 W4 W5 W3 W3 W5 W7
a6 =
5040
240
24
144
8
24
720
i
h
1
17 + 51t2 + 51t4 + 17t6 .
=
315
(8.25)
r4 =
r4 =
2
15
V8
V3 V6 V3 V4 V4 V5 V3 V3 V5 V7
17
r6 =
+
+
+
+
+
+
=
.
5040 240
24
144
8
24
720
315
Note the absence of terms in t2 , t4 or t6 in the coefficients r2 , r4 or r6 .
(8.27)
8.6
8.6
(8.28)
1
1
1
tan3 + tan5 tan7 + .
3
5
7
(8.29)
h
i
e5 t5 1 + 5a2
e2
tan5 =
e 6 t6 [ 1 ]
tan6 =
e 7 t7 [ 1 ]
tan7 =
(8.30)
so that
e +
= t
e5 t
e7 t
e3 t
3a2 t2 +
5a4 5a2 t2 +t4 +
7a6 7 a4 +a22 t2 +7a2 t4 t6 .
3
5
7
(8.31)
e = c.
(8.32)
H3 = 2 2 ,
H5 = 4 (11 24t2 ) 3 (11 36t2 ) + 2 (2 14t2 ) + t2 ,
H 7 = 17 26t2 + 2t4 .
8.7
(8.33)
(8.34)
Comparing this equation with (8.28) we see from equation (8.31) that
=x
e t1 +
x
x
x
e 3 t1
e 5 t1
e 7 t1
3r2 t21 +
5r4 5r2 t21 +t41 +
7r6 7 r4 +r22 t21 +7r2 t41 t61 .
3
5
7
(8.35)
8.7
x
e=
x
,
1
(8.36)
where
K3 = 212 31 t21 ,
K5 = 14 (11 24t21 ) 313 (8 23t21 ) + 512 (3 14t21 ) + 301 t21 + 3t41 ,
K 7 = 17 77t21 105t41 45t61 .
(8.37)
Setting 1 = 1 /1 gives the Redfearn series except that he writes rather than a more
correct 1 etc.
8.8
x sec (, )
,
cos
(8.38)
where x is given in equation (8.19) and we evaluate sec by using the binomial series (E.28) and the expressions for tann given in (8.30):
x
1 e2
1 e4
1 e6
=1+
W3 +
W5 +
W7 ,
(8.39)
c
2
24
720
1/2
1
1
1
sec = 1 + tan2
= 1 + tan2 tan4 +
tan6 +
2
8
16
1 e6
1 e2 2 1 e4
16a4 t2 + 8a22 t2 8a2 t4 + t6
=1+
t + 8a2 t2 t4 +
2
8
16
(8.40)
1 e2
1 e4
k(, ) = 1 +
W 3 + t2 +
W5 + 6t2 W3 + 24a2 t2 3t4
2
24
e6
+
W7 +15t2 W5 +360a2 t2 W3 45t4 W3 +720a4 t2 +360a22 t2 360a2 t4 +45t6
720
(8.41)
Finally, using the W -coefficients from (7.21 7.22) and the a-coefficients from (8.25),
1 e2
1 e4
1 e6
k(, ) = 1 +
H2 +
H4 +
H 6,
2
24
720
e = c,
(8.42)
H2 =
H4 = 4 3 (1 6t2 ) + 2 (1 + 24t2 ) 4t2
H 6 = 61 148t2 + 16t4 .
(8.43)
8.8
8.9
1
1
.
cos x sec
(8.44)
Ignoring the factor of cos for the moment we follow the same steps as in the previous
section but calculate sec using the expression for tan given in equation (8.34):
1
1 2
1 4
1 6
x =
1 x
e V3 x
e V5
x
e V7 , ,
(8.45)
1 c1
2
24
720
1/2
1
1
1
sec = 1 + tan2
= 1 + tan2 tan4 +
tan6 +
2
8
16
1 6
1 2 2 1 4
e t1 + x
e 8r2 t21 t41 + x
e 16r4 t21 + 8r22 t21 8r2 t41 + t61 ,
=1+ x
2
8
16
(8.46)
1 h
1 2
1 4
1+ x
e V3 + t21 + x
e V5 6t21 V3 + 24r2 t21 3t41
1 c1
2
24
i
x
e6
+
V7 15t21 V5 360r2 t21 V3 +45t41 V3 +720r4 t21 +360r22 t21 360r2 t41 +45t61
720
1
1 2
1 4
1 6
=
1+ x
e p2 + x
e p4 +
x
e p6 ,
(8.47)
1 c1
2
24
720
x sec =
where the coefficients and their spherical limits are evaluated using (7.35, 7.39) and (8.27):
p2 = 1 t21
p2 = 1 t21
p4 = 5 + 14t21 + 9t41
(8.48)
f () = c,
f 0 () = 0 cs,
f 00 () = 00 c2 0 sc,
f 000 () = as,
and
1
1
cos = 1 c1 +[10 c1 1 s1 ](1 )+ [100 c1 210 s1 1 c1 ](1 )2 + [as1 ](1)3 .
2!
3!
(8.50)
8.9
(8.51)
(8.53)
q2 = t21
q 2 = t21 ,
q 4 = 8t21 3t41 ,
(8.54)
(8.55)
g2 = p2 + q2 = 1
g 2 = 1
g4 = 5
g 6 = p6 + q 6 + 15(p2 q 4 + q 2 p4 ) = 61.
(8.56)
The actual scale factor (8.44) is the inverse of (8.55); using (E.32) we have,
1 2
1 4
1 6
k(x, y) = 1 + x
e K2 + x
e K4 +
x
e K 6,
2
24
720
x
e=
x
1
(8.57)
where
K2 = g2
= 1 ,
K4 = g4 + 6g22
K 6 = g 6 + 30g 2 g 4 90g 32 = 1.
(8.58)
Setting 1 = 1 /1 gives the Redfearn series except that in the sixth order term he has a
denominator of 13 31 as against 16 here. Since we assume the spherical limit for this term
both could be replaced by a6 and there is no inconsistency.
8.10
8.10
Direct series
As for NMS, TMS and NME simply multiply (7.46) and (7.47) by a factor of k0 .
"
#
e3
e5
e7
e + W3 + W5 + W7 ,
e = c
x(, ) = k0
3!
5!
7!
"
#
e2 t
e4 t
e6 t
e8 t
y(, ) = k0 m() +
+
W4 +
W6 +
W8 .
2
4!
6!
8!
(8.59)
(8.60)
Inverse series:
Set xx/k0 and replace x
e = x/1 by x
b = x/k0 1 . The footpoint latitude, 1 , must be
found from (7.4) or (7.5) with yy/k0 : it is the solution of m(1 ) = y/k0 . Equations 7.48,
7.50 and 7.49 become
(x, y) =
x
b
x
b3
x
b5
x
b7
V3
V5
V7 ,
c1 3! c1
5! c1
7! c1
x
b=
(x, y) = 1
x
b 2 t1 x
b 4 t1
x
b 6 t1
x
b 8 t1
V4
V6
V8 .
2 c1
4! c1
6! c1
8! c1
(x, y) = 1
x
b 2 1 t1 x
b 4 1 t1
x
b 6 1 t1
x
b 8 1 t1
U4
U6
U 8,
2
4!
6!
8!
x
k0 1
(8.61)
y
k0
(8.62)
m(1 ) =
(8.63)
(8.64)
2
24
720
e 3 t H3 + 1
e 5 t H5 + 1
e7 t H 7 ,
e + 1
(, ) = t
3
15
315
1 4
1 6
1 2
k(x, y) = k0 1 + x
b K2 + x
b K4 +
x
b K6 ,
2
24
720
1 3
1 5
1 7
(x, y) = x
b t1 + x
b t1 K 3 + x
b t1 K 5 +
x
b t1 K 7 ,
3
15
315
(8.65)
m(1 ) =
x
b=
y
, (8.66)
k0
x
. (8.67)
k0 1
Note 1: both series for the scale factor give a value of k0 on the central meridian.
Note 2: as usual c= cos , t= tan , =()/() from (5.53) and the 1 subscript
denotes a function evaluated at the footpoint latitude. For convenience all of the required
coefficients are collected on the following page.
8.11
All coefficients
x(, )
W 3 = t2
W5 = 4 3 (1 6t2 ) + 2 (1 + 8t2 ) 2t2 + t4
W7 = 61 479t2 + 179t4 t6 + O(e2 )
y(, )
W4 = 4 2 + t2
W6 = 8 4 (1124t2 ) 28 3 (16t2 ) + 2 (132t2 ) 2t2 + t4
W8 = 1385 3111t2 + 543t4 t6 + O(e2 )
(x, y)
V3 = 1 + 2t21
V5 = 413 (1 6t21 ) 12 (9 68t21 ) 721 t21 24t41
V7 = 61 + 662t21 + 1320t41 + 720t61
(x, y)
V4 = 412 91 6t21
V6 = 814 (1124t21 )8413 (38t21 )+22512 (14t21 )+6001 t21 +120t41
V8 = 1385 7266t21 10920t41 5040t61
(x, y)
k(, )
H2 =
H4 = 4 3 (1 6t2 ) + 2 (1 + 24t2 ) 4t2
H 6 = 61 148t2 + 16t4
(, )
H3 = 2 2
H5 = 4 (11 24t2 ) 3 (11 36t2 ) + 2 (2 14t2 ) + t2
H 7 = 17 26t2 + 2t4
k(x, y)
K2 = 1
K4 = 413 (1 6t21 ) 312 (1 16t21 ) 241 t21
K6 = 1
(x, y)
K3 = 212 31 t21
K5 = 14 (11 24t21 ) 313 (8 23t21 ) + 512 (3 14t21 ) + 301 t21 + 3t41
K 7 = 17 77t21 105t41 45t61
(8.68)
8.12
Chapter
9.1
So far we have rather casually mixed the terms projection coordinates and grid coordinates. We must now be a little more precise for the two are logically distinct and rarely
equal. In fact we must consider four reference systems:
We will discuss each of these points in turn with examples relating to the National Grid of
Great Britain (NGGB) and the Universal Transverse Mercator (UTM).
Geographic coordinates
The basic data of any survey are the latitude and longitude of all features of interest. It is
important to remember that in each survey the (, ) coordinates are defined with respect
to one particular ellipsoid of revolution. For example the NGGB projection is based on the
Airy 1830 ellipsoid and UTM is based on the International 1924 (aka Hayford 1909) ellipsoid. Remember that there is no such thing as an absolute value of latitude or longitude
the same location has different geographic coordinates with respect to different ellipsoids.
9.2
Projection coordinates
We consider here only the modified Transverse Mercator projection coordinates derived
from the latitude and longitude values and the parameters of the chosen ellipsoid by the
Redfearn formula for x(, ) and y(, ) in equations (8.59, 8.60). These formulae give
Cartesian coordinates with a projection origin where the central meridian crosses the equator. Our convention (in agreement with Snyder) is that the positive y-axis is directed northwards along the projection of the central meridian and the positive x-axis is eastwards along
the projection of the equator. (Beware of alternative conventions. For example the original
paper of Redfearn, and most continental texts, have x and y interchanged).
In the introduction we called the Cartesian representation defined by equations (8.59)
and (8.60) a super-map. From (8.59) we see that if we restrict to be less than say 5 ,
or about 0.1 radians, then the x-range of this super-map will of the order 106 metres and
from (8.60) we see that the distance between the equator and pole on the central meridian
is k0 m(/2) 107 metres. Actual printed maps are then obtained by scaling the supermap by the representation factor (RF). For example the Landranger maps produced by the
Ordnance Survey of Great Britain (OSGB) have an RF of 1:50000 so that a sheet which
measure 80cm square represents an area approximately 40km square on the super-map or
on the ground. The caveat approximately is necessary since we know that the scale in the
TME projection is the complicated function given by equation (8.64), approximately unity
when close to the central meridian.
Note that the maps produced in this way need not be embellished by lines of any kind.
There is no obligation to show the x- and y-axes or the origin of the projection coordinates.
There is also no obligation to show any lines marking constant values of x and y. There is
no obligation to show lines (curves) of constant latitude or longitude.
9.3
coordinates referred to a true origin on the central meridian will also take positive and
negative values. This is unsatisfactory for many practical applications and we therefore
introduce a false origin of the grid at a position such that the grid coordinates relative to
that point are positive throughout the region of interest. These positive coordinates, rounded
to the nearest metre, are called Eastings (E) and Northings (N). The choice of false origin
for UTM and NGGB is described in the next sections.
9.2
!
Figure 9.1
E0 = 500000
(
0
(N)
N0 =
10000000 (S)
(9.1)
(9.2)
where x, y are calculated (to the nearest metre) using the modified TME series (8.59, 8.60)
taking relative to 3 W (in radians) and k0 =0.9996. Note that (E0, N 0) are the coordinates of the true origin (at x = y = 0) relative to the false origin on each grid. The
9.4
figure shows the eastings and northings of several points on the perimeter of the zone: (a) at
= 0 on the equator we have E = 833960, N = 0; (b) at = 0 on the parallel at 84 N
we have E = 534993, N = 9329291. Thus the width of the UTM zone at the equator is
approximately (recall k0 6= 1) 668km and at latitude 84 N it is approximately 70km.
The inverse relations are the series (8.61, 8.63) with the replacement of x by E E0
and coefficients evaluated at a footpoint latitude such that
m(1 ) =
y
N N0
=
.
k0
k0
(9.3)
Each zone is split into twenty latitude sub-zones, nineteen of extent 8 starting from
80 S and one of 12 finishing at 84 N. Each of these twenty latitude bands is designated a zone letter from C to X, with I and O excepted (to avoid ambiguity with digits
1 and 0). My home (approximately 55 570 N, 3 110 W) is just within the northern limit
of zone 30U, centred on 3 W and running from 48 to 56 N. The extent of this subzone in projection coordinates is approximately 890km north-south whilst the width
varies from 447km to 373km as one goes north.
Within each of the sub-zones the 100km squares are labelled with row and column
letters; for example Edinburgh is in a 100km square labelled UG. Thus 30UUG fixes
my home to within 100km. (For a full description of the labelling of the 100km
squares see the DMA manual listed in the bibliography.)
Within a 100km square the Eastings and Northings range from 0 to 99999m so that
1m accuracy is given by two five digit numbers. Thus the full UTM grid reference of
my armchair is 30UUG 88533 00666.
Such precision is often superfluous and a pair of numbers with 4, 3, 2, 1 digits may
be used for accuracy to within 10m, 100m, 1km, 10km (of the left hand and bottom
edges of a box of that size). Thus
30UUG 8853 0066 fixes my home to within 10m
30UUG 885 006 fixes my home to within 100m
30UUG 88 00 fixes the centre of Edinburgh within 1km
30UUG 8 0 fixes Edinburgh within 10km
9.5
9.3
The British national grid (NGGB) is a grid overlain on the TME projection centred on
longitude 2 W and based on the Airy 1830 ellipsoid, for which the parameters are given in
Section 5.2. It is a modified projection with a value of k0 = 0.9996012717 on the central
meridian. This value arose when the 1936 re-survey was constrained to be as close as
possible to the previous survey at a number of selected points. For most practical purposes,
and in the remainder of this chapter, we will take k0 0.9996.
Clearly only a small area of the projection is
needed, just the small box shown in Figure 9.2,
so we take the true origin at latitude 0 = 49 ,
or 0.855rad; this parallel is slightly south of the
area mapped by the OSGB. The false origin is
then chosen west and north with E0=400000m
and N0=100000m so that mapped area is completely to its east and north and all E and N values
are positive. Recall that E0 and N0 give the position of the true origin relative to the false origin.
1.4
1.2
60N
1
30W
45N
0.8
false
true
9.6
12W
10W
8W
6W
4W
2W
2E
4E
6E
62N
1300
HL
HM
HN
HO
HP
JL
JM
HQ
HR
HS
HT
HU
JQ
JR
HV
HW
HX
HY
HZ
JV
JW
NA
NB
NC
ND
NE
OA
OB
NF
NG
NH
NJ
NK
OF
OG
NL
NM
NN
NO
NP
OL
OM
1200
60N
1100
1000
58N
800
56N
700
NQ
NR
NS
NT
NU
OQ
OR
NV
NW
NX
NY
NZ
OV
OW
SA
SB
SC
SD
SE
TA
TB
SF
SG
SH
SJ
SK
TF
TG
SL
SM
SN
SO
SP
TL
TM
SQ
SR
SS
ST
SU
TQ
TR
SV
SW
SX
SY
SZ
TV
TW
Latitude
900
600
500
54N
400
300
52N
200
100
FALSE
ORIGIN
50N
k=0.9996
k=1.0004 k=1
-100
k=1
TRUE
ORIGIN
k=1.0004
48N
-200
-100
100
200
300
400
500
Figure 9.3
600
700
800
9.7
on a slightly modified Airy ellipsoid centred on the meridian at 8 W; the grid having a
true origin on this meridian at latitude 53 300 N and a false origin such that E0=200000m
and N0=250000m. The scale modification is k0 = 1.000035. (See the Bibliography for a
reference to a discussion of this unusual factor).
Returning to NGGB we see that the equations for E and N must take into account the
shift along the meridian to the true origin as well as the translation to the false origin:
E = E0 + x(, )
= x + 400000,
(9.4)
(9.5)
where we have calculated m(0 ) with 0 = 49 = 0.8552 radians using (5.71). For x
and y we use the Redfearn series (8.59, 8.60) with the understanding that (in radians)
is measured from the central meridian at 2 W. The inverse series are again calculated
from (8.61, 8.63) with x replaced by E E0 and the footpoint calculated from
m(1 ) =
y
N N0
=
+ m(0 )
k0
k0
(9.6)
Using the above formulae I calculate that, to the nearest metre, I am sitting at E=325701,
N=673642. This is in the 100km square labelled NT so the alpha-numeric grid reference
is NT 25701 73642 to within 1m. As for UTM we can use shorter grid references such as
NT25707364, NT257736, NT2573, NT27 for accuracy to within 10m, 100m, 1km, 10km.
9.4
x
b=
x
.
k0 a
(9.8)
This scale factor is a function of x only and it is straightforward to calculate (using 9.7) the
values at which k = 1 and k = 1.0004. Taking k0 = 0.9996 and the radius of the sphere
as a = 6378km, (approximately equal to the semi-major axis of the Airy ellipsoid), we
find that these isoscale lines are at approximately x=180km and x = 255km respectively.
Figure 9.4 shows a plot of the scale factor as a function of x. We shall prove that the
differences between the scale factors of TMS and TME are so small that the corresponding
plot for TME at any fixed y value is indistinguishable from Figure 9.4.
9.8
1.0015
scale k
1.0010
1.00098
1.0007
1.0005
1.0000
0.9995
0
50
100
150
200
250
300
350
400
x=|E-E0| in kilometres
Figure 9.4
Scale variation in TME
The corresponding result for TME is given by equation (8.66):
1 2
1 4
1 6
k(x, y) = k0 1 + x
b K2 + x
b K4 +
x
b K6
2
24
720
x
b=
x
,
k0 1
(9.9)
(9.10)
where 1 and t1 are evaluated at the footpoint latitude such that m(1 ) = y/k0 and y is
related to the northing coordinate by equation (9.2) or (9.5) for UTM or NGGB respectively.
To simplify the comparison of the TMS and TME scale factors it is instructive, and
useful, to write K2 and K4 in terms of the following parameter of order e2 :
2 = 1 =
e2 cos2
1=
= e02 cos2 .
1 e2
(9.11)
The importance of this parameterisation is that it allows us to identify the way in which the
K-coefficients depend on the eccentricity of the ellipse. It was not introduced at an earlier
stage because the parameter is much better suited to the derivations of all the coefficients.
For the non-trivial coefficients in the TME scale factor we have
K2 = 1 + 12
K4 = 1 + 612 + 14 (9 24t21 ) + 16 (4 24t21 )
(9.12)
Clearly as e2 0 the TME scale factor reduces to the TMS scale factor. Their difference
is given by terms of order 2 x
b 2 and 2 x
b 4 : these are typically no more than 105 and 108
respectively so that Figure 9.4 is an adequate representation of the TME scale variation
9.9
with x for any value of 1 calculated for the footpoint value corresponding to y. (Note that
in working with the Redfearn series we must not assume that the t1 terms are small: the
footpoint latitude may be of the order 80 for which t1 = tan 1 5.7 so that t21 32).
Figure 9.4 is annotated with the scale factors for the extreme cases which may arise in
TME. For UTM the greatest extent in projection coordinate is on the equator where x
334km and k reaches its worst value of almost 1.001still perfectly acceptable in all practical applications. For NGGB the worst cases are the extremes of the East Anglia coast and
the Outer Hebrides where |x| is about 300km and k1.0007. The bulk of the NGGB grid is
within approximately 255km of the central meridian where |k 1|<0.0004. Note that a typical map sheet corresponds to a small section of the plot. For example on my local sheet of
the OSGB, bounded by E=316km and E=356km and close to the central meridian on which
E=400km, the scale varies from k=0.999686 on the west to k=0.999624 on the east.
Isoscale lines of TME
We will now investigate the lines on which k is a constant and show how little they differ
from the straight lines x=constant which we found for the TMS projection. To do this we
make a Lagrange inversion of the series (9.9) for the TME scale factor. Start by writing this
equation as
2
1 4 K4
1 6 K6
x
k
w
1 =x
b2 + x
b
+
x
b
x
b=
.
(9.13)
K2 k0
12 K2 360 K2
k0 1
From equations (B.10B.11), with z = x
b 2 and b4 = 0, we immediately obtain the inverse
!
2
5K
K
K
K
k
2
6
4
2
2
2
3
4
x
b =w
1 ,
(9.14)
w
w
w=
12K2
K2 k0
360K22
from which we can find x = k0 1 x
b , and hence (E E0), as a function of k for a given value
of N . We have used this result to calculate the eastings at which k = 1 and k = 1.0004 for
the UTM projection at four particular values of N for each of which we have first calculated
the footpoint from equation (9.3). From the fourth column we see that the deviation in
N
9000000
6000000
3000000
0
k=0.9996
k=1.0
k=1.0004
footpoint
180946
255887
180556
255336
180010
254564
179759
254208
81 .05846848
54 .14694587
27 .12209299
0 .00000000
eastings on the k = 1 isoscale is just approximately 1km over in a range of 9000km from
the equator to just over 80 N. The deviation of the k = 1.004 isoscale is still under 2km.
Thus the isoscale lines are essentially parallel to the central meridian.
9.10
footpoint
1200000
k=0.9996
k=1.0
k=1.0004
180369
255275
180302
255180
180231
255080
180157
254976
60 .68382350
800000
57 .09116995
400000
53 .49645197
49 .89956809
Clearly the lines on which k = 1 and k = 1.0004, which are shown in Figure 9.3, are
indistinguishable from straight lines parallel to the central meridian. For the k = 1 locus
the change in eastings is only 212m over the full northings extent of 1200km.
To end this discussion of the scale factor in TME note that much of mainland Britain is
just within the scale variation of 0.9996 to 1.0004. This may be why these numbers were
adopted as suitable criteria for mapping in the first place. This is of course speculation and
we would be interested to hear of any evidence either way.
9.5
When we first discussed convergence in Section 3.6 we observed that on any particular
meridian (on the TMS projection) the convergence, defined as the angle between grid north
and true north (the tangent to the meridian), must increase from zero at the equator to
at the pole. The same must be true for TME although we require only values up to the
northerly limit of UTM or NGGB. The following table shows the convergence for these two
projections for several latitude valueson a bounding meridian for UTM (at 0 +3 ) and on
the extreme meridian intersecting the land area covered by NGGB.
Convergence along a projected meridian
UTM at 0 +3
NGGB at 7 W
84 N
2 590 100 W
60 N
4 190 5800 E
80 N
2 570 1600 W
58 N
4 140 3600 E
60 N
2 350 5500 W
56 N
4 80 5500 E
40 N
1 550 4600 W
54 N
4 20 5500 E
20 N
1 10 3700 W
52 N
3 560 3800 E
50 N
3 500 300 E
0 N
For UTM the values clearly approach the limiting value of 3 which would be attained
at the pole. For NGGB the specified meridian is 5 west of the central meridian and the
convergence clearly approaches this value at the northern extremity of the grid.
9.11
The convergence varies only slightly over any one of the map sheets. For example exact
calculations give the convergence at the corners of the NGGB Edinburgh sheet bounded by
E=316km, E=356km, N=650km and N=690km as
Boundary of Edinburgh sheet
E
NW
SW
316
316
00
00
1 7 14.94 E
690
1 6 20.85 E
650
00
356
690
NE
00
356
650
SE
35 13.82 E
34 45.48 E
Note that the variation of convergence from top to bottom is much less than the variation
from east to west. Similar figures are given at the corners of every map in the OSGB series:
the values at other points may be approximated by interpolation.
It is important to observe that convergence () values are not vanishingly small and
they must be taken into account in relating an azimuth () to a grid bearing () by the
relation = + discussed in Section 8.2. This correction is important in high accuracy
applications.
9.6
One obvious test of the accuracy of the TME transformations is to start from given geographical coordinates (, ), transform to projection coordinates with the direct Redfearn
series (8.59, 8.60) and then reverse the transformation with the inverse series (8.61, 8.63) We
should then be back where we started. Before doing so we must decide on our standard of
accuracy. We shall work to within 1mm in the projection coordinates and to within 0.000100
in geographical coordinates. These accuracies are approximately equivalent, for we see
from (2.5) that 0.000100 is equivalent to 3mm along the meridian and less than 2mm along
a parallel for examples calculated within the region of the NGGB (at latitudes from 50 N
60 N). To cope with rounding errors we compute to two extra places of decimals. These
accuracies are purely to assess the mathematical consistency of our transformations for, in
practice, no survey claims accuracies better than 10cm. The following example shows that
this first test is satisfied with flying colours. Note that we use eastings and nothings rather
than x and y, the projection coordinates. For NGGB they are related by (9.4, 9.5) and the
footpoint is to be calculated from (9.6).
Redfearndirect
Redfearn-Inverse
Lat
Lon
E
N
00
0
00
52 39 27.2531 1 43 4.5177
651409.903
313177.270
E
N
Lat
Lon
651409.903
313177.270 52 390 27.253100 1 430 4.517700
Another test is to assess the outcome of small changes in the inputs to the direct and
inverse series. Sticking to the same coordinates as above we perturb the geographical coordinates by 0.000100 in latitude, longitude separately and together: for the inverse we perturb
the projection coordinates by 0.001mm. The results are shown overleaf.
9.12
NGGB-direct
Lat + 0.000100
Lon + 0.000100
Both together
Lat
52 39 27.253100
52 390 27.253200
52 390 27.253100
52 390 27.253200
Lon
1 43 4.517700
1 430 4.517700
1 430 4.517800
1 430 4.517800
E
651409.903
651409.903
651409.905
651409.905
N
313177.270
313177.273
313177.270
313177.274
NGGB-Inverse
E + 1mm
N + 1mm
Both together
E
651409.903
651409.904
651409.903
651409.904
N
313177.270
313177.270
313177.271
313177.271
Lat
52 39 27.253100
52 390 27.253100
52 390 27.253100
52 390 27.253100
Lon
1 43 4.517700
1 430 4.517800
1 430 4.517700
1 430 4.517800
Thus we see that 0.000100 changes induce a maximum change on the projection of no more
than 4mm: for the inverse transformation 1mm changes in the projection coordinates change
the geographical coordinates by no more than 0.000100 .
Now when Redfearn published his series he was simply extending the series that had
been published earlier by Lee who had discarded terms smaller that 4 e2 , 5 e2 , (x/a)4 e2
and (x/a)5 e2 in the series for x, y, and respectively. Redfearn observed that the coefficients in Lees series were increasing rapidly, particularly at larger latitudes where t and t1
are not small, and consequently it seemed possible that some omitted terms might actually
be larger than the smallest ones retained. This in fact proved to be the case for two of the
terms omitted by Lee. Redfearns analysis to higher order makes clear which terms can be
safely omitted, as is done in many published expressions for the transformations.
To compare the size of the terms we again introduce the parameter 2 , which is O(e2 ),
defined in equation (9.11). The transformation equations are those of Section 8.10 but we
now write them in terms of eastings and northings using equations (9.49.6) and we also
replace by 0 .
"
#
e3
e5
e7
e + W3 + W5 + W7 ,
E E0 = k0
(9.15)
3!
5!
7!
e2 t
e4 t
e6 t
e8 t
N N 0 + k0 m(0 ) = k0 m() +
+
W4 +
W6 +
W8 ,
2
4!
6!
8!
x
b3
x
b5
x
b7
x
b
(E, N ) =
V3
V5
V7 ,
c1 3! c1
5! c1
7! c1
(E, N ) = 1
x
b 2 1 t1 x
b 4 1 t1
x
b 6 1 t1
x
b 8 1 t1
U4
U6
U 8,
2
4!
6!
8!
(9.16)
(9.17)
(9.18)
where
e = ( 0 )c ,
x
b=
E E0
,
k0 1
m(1 ) =
N N0
+ m(0 )
k0
(9.19)
and the coefficients, given in equations (7.517.53), are now written in terms of as shown
overleaf.
9.13
The significance of the vertical rules in the rewritten coefficients will be discussed shortly.
W 3 = 1 t2 + 2
W4 = 5 t2 + 9 2 + 4 4
W5 = 5 18t2 + t4 + 2 (14 58t2 ) + 4 (13 64t2 ) + 6 (4 24t2 )
W6 = 61 58t2 + t4 + 2 (270 330t2 ) + 4 (445 680t2 )
+ 6 (324 600t2 ) + 8 (88 192t2 )
W7 = 61 479t2 + 179t4 t6
W8 = 1385 3111t2 + 543t4 t6
(9.20)
V3 = 1 + 2t21 + 12
V5 = 5 28t21 24t41 12 (6 + 8t21 ) + 14 (3 4t21 ) + 16 (4 24t21 )
V7 = 61 + 662t21 + 1320t41 + 720t61
(9.21)
(9.22)
9.14
In: = 58 = 7
0
2
1
e 104482.705
e3
e5
164.425
e7
4.9E-06
e6
e8
Out: E =104647.323m
8
-0.199
0.389
0.003
In: = 58 = 7
0
2
0
e 901166.420
e2
e4
6.0E-06 4.3E-09
Out: N =912106.244m
8
10935.050
4.753
0.032
2.8E-05
-7.1E-13
-1.6E-05
14
x
b 0 58 50 53.772800
x
b2
50 54.571800
0.804200
0.002500
x
b6
0.002700
1.0E-0500
x
b8
1.1E-0500
x
b
-8.9E-0700
-2.1E-0800 -9.4E-1200 2.3E-1400
Out: =7 00 0.000000
14
16
18
x
b 1 7 00 39.421300
x
b3
39.571100
0.012000
x
b5
0.162600
-3.4E-0500
00
x
b
-1.8E-0800 -2.6E-1000
0.0008
Lee would no doubt have had this in mind when he dropped the sixth order terms in his
calculations. To be exact, of the required terms he dropped the term in 6 in the direct
series for y, the term in (x/a)6 in the inverse series for and the term in (x/a)7 in the
inverse series for ; at the same time he included the term in (x/a)5 e2 in the inverse series
for although we now see that it is negligible.
The approximations used in Snyder are similar but he retains the terms 12 (46252t21 ) in
U6 ; we have shown this is negligible for NGGB but it is required at the higher latitudes that
we meet in UTM. (Comment: Snyder uses e02 explicitly in some coefficients. To compare
with the expressions above first set e02 = 2 sec2 = 2 [1 + t2 ]).
9.15
9.7
We have stressed that although no penalty is incurred by using the full Redfearn series
as summarised in Section 8.10, the coefficients truncated at the vertical rules in equations (9.209.22) will be perfectly adequate. Dropping the small terms equations (9.15
9.19) become
"
#
e3
e5
e + WT + WT ,
E(, ) = E0 + k0
(9.23)
3! 3
5! 5
"
N (, ) = N 0 + k0 [m() m(0 )] + k0
(E, N ) =
#
e2 t
e4 t
e6 t
+
W4T +
W6T ,
2
4!
6!
x
b
x
b3 T
x
b5 T
x
b7 T
V3
V5
V ,
c1 3! c1
5! c1
7! c1 7
(E, N ) = 1
(9.24)
(9.25)
x
b 2 1 t1 x
b 4 1 t1 T x
b 6 1 t1 T
U4
U6 ,
2
4!
6!
(9.26)
where
e = ( 0 )c,
x
b=
E E0
,
k0 1
m(1 ) =
N N0
+ m(0 ),
k0
(9.27)
and, as usual, c= cos , t= tan , =()/() from (5.53) and the 1 subscript denotes
a function evaluated at the footpoint latitude.
The truncated coefficients follow from equations (9.209.22)
W3T = 1 t2 + 2
W4T = 5 t2 + 9 2
W5T = 5 18t2 + t4 + 2 (14 58t2 )
W6T = 61 58t2 + t4
(9.28)
V3T = 1 + 2t21 + 12
V5T = 5 28t21 24t41
T
(9.29)
e2 cos2
1=
= e02 cos2 .
1 e2
(9.30)
(9.31)
9.16
9.8
The published form of the series used by the OSGB uses a different notation for the truncated series of the previous section. First use equations (5.78, 5.80) to set
M = k0 [m() m(0 )]
(9.32)
5
5
21
= k0 b 1+n+ n2 + n3 (0 ) 3n+3n2 + n3 sin(0 ) cos(+0 )
4
4
8
35 3
15 2 15 3
n + n sin 2(0 ) cos 2(+0 )
n sin 3(0 ) cos 3(+0 )
+
8
8
24
where n is defined in equation (5.5) as (a b)/(a + b).
Equations (9.24), (9.23), (9.26) and (9.25) may be written as
N = I + II( 0 )2 + III( 0 )4 + IIIA( 0 )6
3
(9.33)
E = E0 + IV( 0 ) + V( 0 ) + VI( 0 ) ,
2
(9.34)
(9.35)
(9.36)
y
N N0
=
+ m(0 )
k0
k0
(9.37)
X=
t1
2e
1 e1
1
e1 c1
II =
esc
2
III =
esc3 W4T
4!
V=
ec3 W3T
3!
VI =
ec5 W5T
5!
VIII =
t1 U4T
4!e
1 e13
IX =
t1 U6T
6!e
1 e15
XI =
V3T
3!e
13 c1
XII =
V5T
5!e
15 c1
IIIA =
XIIA =
esc5 W6T
6!
V7T
7!e
17 c1
and, as usual, c= cos , t= tan , =()/() from (5.53) and the 1 subscript denotes
a function evaluated at the footpoint latitude. The OSGB publication exhibits the above
formulae without the tildes on and because the definitions of and in that publication
already absorb a factor of k0 .
Chapter
10
IN PREPARATION
10.2
11
Chapter
This chapter is a not directly concerned with the Mercator (or any) transformation, rather
it tackles the problem of finding geodesics, curves of shortest length, between points on a
sphere and on an ellipsoid. To be precise we solve the following two problems on both the
sphere and the ellipsoid.
1. The direct problem.
Define a geodesic by an initial azimuth 1 at a starting point P1 with geographic
coordinates (1 , 1 ). Find the coordinates (2 , 2 ) of a point P2 at a distance s
measured along the geodesic; find also the azimuth 2 of the geodesic at P2 .
2. The inverse problem.
Given the coordinates (1 , 1 ) and (2 , 2 ) of two points P1 and P2 find s, the geodesic distance between them, and also 1 , 2 , the azimuths of the geodesic at the
points.
Note that we need only relative longitudes so we calculate 1 in the direct
problem; for the inverse problem only is given.
For the case of the sphere we shall see that we can derive the results (a) by spherical
trigonometry or (b) by integrating in closed form simple first order differential equations
describing the geodesics. For the ellipsoid there is no equivalent to spherical trigonometry and we have no option but to integrate the differential equations; moreover, since the
equations cannot be integrated in closed form we must resort to series expansions.
The results for the sphere have been long established and approximations to the ellipsoid
problems have been presented by many authors over the last 150 years. The results we
present here were given in the Survey Review by Vincenty. (See bibliography)
11.2
Appendix
A.1
Planar curves
A straight line has zero curvature. The curvature, , of a general curve in the plane is
defined as the rate of change of the direction of its tangent with respect to the distance
travelled along the line:
d
=
.
(A.1)
ds
If we are given the equation of the curve as y = f (x) with respect to Cartesian axes then it
is natural to choose the x-axis as the reference for the direction of the tangent.
dy
= y 0 (x),
dx
cos =
dx
ds
(A.2)
d
d[y 0 (x)]
dx
=
= y 00 (x)
ds
ds
ds
= y 00 (x) cos .
(A.3)
Figure A.1
p
2
02
Now sec = 1 + tan = 1 + y so we obtain d/ds and
=
y 00 (x)
.
[1 + y 02 ]3/2
DASH
d
dx
(A.4)
The choice of sign is a matter for convention in every case. We shall illustrate this point
immediately.
A.2
y 0 (x) =
x
,
a2 x2
y 00 (x) =
a2
.
(a2 x2 )3/2
(A.5)
Substituting these values in equation (A.4) we see that the curvature of the upper semicircle
is = (1/a) whilst for the lower semicircle = (1/a). Now it is conventional to
define the curvature of a circle to be positive so we must choose the negative sign in the
definition for the case of the upper semicircle and the positive sign for the lower; with these
choices of sign we have a constant curvature = 1/a. Therefore the curvature of a circle
is the inverse of the radius and vice-versa.
The osculating circle and the radius of curvature
The particular circle which touches the curve at P (Figure A.1) and also shares the same
curvature at that point is called the osculating circle (osculating=kissing) or the circle of
curvature. The radius of this circle defines R, the radius of curvature of the curve at that
point. Clearly
R=
1
.
(A.6)
dy
y
= ,
dx
x
d y du
du x dx
1
x
y y
x
= =
R
[x 2 + y 2 ]3/2
x
y y
x
x 3
DOT
d
.
du
(A.7)
A.3
Set DASH d
ds
DASH
d
ds
(A.8)
Set DOT d
d
dy
Since is the angle between the tangent and the x-axis we have tan =
= y/
x
dx
x
y y
x
Differentiating with respect to gives sec2 =
.
2
x
x 2 + y 2
But we also have sec2 = 1 + tan2 =
x 2
2
Therefore we must have x
y x
y = x + y 2 so that (A.7) becomes
1
1
= =
2
R
[x + y 2 ]1/2
DOT
d
d
(A.9)
Figure A.2
y = b sin U,
y = b cos U,
y = b sin U,
(A.11)
giving
+ab
1
1 e2
= 2
=
.
a [1 e2 cos2 U ]3/2
[a (a2 b2 ) cos2 U ]3/2
(A.12)
A.4
A.2
First consider two neighbouring points, P (s) and Q(s+s), on a curve parameterised by its
arc length s (Figure A.3a). The chord length between these points is given by s2 = r r.
The tangent vector at P is the limit of r/s and therefore has the properties
t = r0 ,
t t0 = 0,
tt=1
(A.13)
Figure A.3
Consider the tangents at neighbouring points (Figure A.3b); t(s) and t(s + s) = t + t are
compared in the third figure by bringing them together at some point B. Tangent vectors
are unit vectors so that |t| = |t + t| = 1; therefore in the limit of s 0 we see that t is
in the direction of n, a unit vector normal to t and in the osculating plane defined by the
two vectors t(s) and t(s + s). Furthermore, if the angle between the unit tangent vectors
is then as s 0 we must have |t| = . Consequently
t
d
=
n = n.
s0 s
ds
t0 = lim
(A.14)
The vector n is called the principal normal to the curve and the curvature , is defined as
for planar curves. We can invert this relation and write
n=
1 0
1
t = r00 .
(A.15)
Note that n is defined to be a unit vector; on the other hand t0 is not a unit vector, its length
is equal to the curvature . Since t t0 = 0 we must have n t = 0.
Now introduce the unit binormal vector b, defined so that it
forms a right-handed orthogonal triad with t and n. Therefore
b = t n,
(A.16)
tn=0
nb=0
b t = 0,
t n = b,
n b = t,
b t = n.
(A.17)
Figure A.4
A.5
(A.18)
Now since t0 = n the first of these terms must vanish so we must have t b0 = 0.
Consequently b0 is perpendicular to both t and b and it is therefore a vector in the direction
of n. The vector b0 is not a unit vector and its magnitude is defined to be , the torsion of
the curve. We choose to set
b0 = n.
(A.19)
The torsion of the curve is a measure of the rate of rotation of the vectors b, and hence n,
about the tangent vector as s increases. The negative sign associates a right-handed rule
as part of the definition.
We have expressions for the derivatives of t and b in equations (A.15) and (A.19). We
now evaluate the derivative of n from b t:
n0 = b0 t + b t0 = n t + b (n) = b t.
(A.20)
This equation together with the derivatives of t and b constitute the set of Frenet-Serret
relations:
t0 = n,
n0 = b t,
b0 = n.
(A.21)
These equations show that the form of a curve in three dimensions is essentially determined
by the two functions (s) and (s) and an initial orthonormal triad.
A.3
Curvature of surfaces
At any point on a surface we can define the curvature of a line on the surface passing
through that point. Rather than build up a large part of differential geometry we shall give
elementary proofs of two important results that we need.
First consider those curves which are defined by the intersection of a plane with the surface. The most important case is a plane which contains the normal N at the point P ; such a
plane defines a normal section of the surface. We shall consider all the normal sections at
a given point and investigate the curvature at P of their intersections. The principal result is
that the maximum and minimum curvatures arise on two normal sections at right-angles to
each other; these are the principal curvatures which we will denote by 1 and 2 . Eulers
formula gives the curvature on any other normal section in terms of the principal curvatures.
The other main result that we need is Meusniers theorem. This relates the curvature
on a normal section to the curvature of the sections made by planes oblique to the chosen
normal plane, i.e. sharing the same tangent at P . It is convenient to prove this theorem first.
A.6
A.4
Meusniers theorem
Without loss of generality we choose axes with the origin O at an arbitrary point on a surface
and such that the xy-plane is tangential at the point. Consider the family of planes which
contain the tangent directed along the x-axis. Each plane intersects the the surface in a plane
curve; let g(x) be the curve on the normal plane and w(x) on an oblique plane inclined at
an angle to the normal plane. If , R denote the curvature and radius of curvature at the
origin of g(x) and w(x) on the normal and oblique planes repectively, then
oblique = sec normal ,
(A.24)
(A.25)
Ax2
Cw2
+Bxw sin +
sin2
2
2
It is clear from this equation that for small x and arbitrary we must have w(x) = O(x2 ). (For suppose on
the contrary that w(x) = O(x), then the LHS of the previous equation would be O(x) and the RHS would be
O(x2 )).
1
w(x) = sec Ax2 + O(x3 ) .
(A.26)
2
Figure A.5
Equations (A.26) and (A.24) give Meusniers theorem
w00
oblique =
= A sec = sec normal .
3/2
0
2
[1 + (w ) ]
x=0
(A.22)
(A.27)
A.7
A.5
(A.29)
Figure A.6
Principal Axes
We exploit the freedom to choose any pair of orthogonal lines as axes in the xy-plane. If
new x0 y 0 -axes are rotated from the original by an angle then we must set
x = cos x0 sin y 0 ,
(A.30)
(A.31)
y = sin x + cos y .
Abbreviate c = cos and s = sin and set x = cx0 sy 0 and y = sx0 + cy 0 in the equation
of the surface (A.28). In terms of these new coordinates
1
1
z = h(x0 , y 0 ) = A(cx0 sy 0 )2 + B(cx0 sy 0 )(sx0 +cy 0 ) + C(sx0 +cy 0 )2 + . (A.32)
2
2
Now the coefficient of the x0 y 0 cross term is equal to [Asc + B(c2 s2 ) + Csc] and this
will vanish if (A C) sin 2 = 2B cos 2 or tan 2 = 2B/(A C). There is always a
solution for and therefore we can always rotate axes so that the equation of the surface
may be taken without a cross term. This defines the principal axes in the tangent plane at
the given point.
A.8
() = 1 cos2 + 2 sin2
(A.35)
(A.36)
for the curvature of the normal section made by a plane making an angle with one of the
principal normal planes.
Without loss of generality let us take 1 > 2 , then we have
1 () = (1 2 ) sin2 0,
2
() 2 = (1 2 ) cos 0.
1 () 2 .
(A.37)
(A.38)
(A.39)
Thus we have proved that the curvatures of normal sections at a point are such that the
minimum and maximum values, the principal curvatures, are associated with orthogonal
planes and the curvature on any other plane is given by the Euler formula. Note that we have
not provided any way of calculating the curvature for an arbitrary surface for in general we
do not have equations for the surface in the form of (A.28). The general study requires the
machinery of differential geometry (see Bibliography) but for surfaces of revolution such
as the ellipsoid we shall find that this not required.
A.9
(A.40)
(A.41)
These definitions are useful in various waysfor example, when we seek to approximate
the surface of small part of the ellipsoid by a sphere.
A.10
Appendix
Lagrange expansions
B.1
Introduction
(B.1)
where both z and w are assumed to be small, less than 1, whilst the coefficients are of
of order unity. The series we shall meet in the cartographic applications will typically be
Taylor series truncated after a few terms. Now since z n z for z < 1 and n > 1 we must
have z w and we might expect it to be represented by a series of the form
z = b 1 w + b 2 w 2 + b 3 w 3 + b 4 w 4 + b5 w 5 + .
(B.2)
One way of finding the coefficients is to substitute the series for z into every term on the
right hand side of (B.1) and compare coefficients of wn on both sides. This is demonstrated
explicitly in the next section. Fortunately a more general method exists, namely the Lagrange expansions defined in Section B.3. This is essential for the inversion of the eighth
order series that we shall encounter.
A second category of problem is illustrated by a series of the form
w = z + c2 sin 2z + c4 sin 4z + c6 sin 6z + ,
(B.3)
where we might have z and w as O(1) whilst the coefficients cn are small. The method of
Lagrange expansions will show that there is an inverse given by
z = w + d2 sin 2w + d4 sin 4w + d6 sin 6w + .
(B.4)
B.2
B.2
The power series may be solved simply by back substitution, i.e. we substitute z from (B.2)
into the terms on the right hand side of (B.1) and compare coefficients of w. If we retain
only terms up to O(w4 ) we have
w = (b1 w+b2 w2 +b3 w3 +b4 w4 ) + a2 w2 (b1 +b2 w+b3 w2 + )2
+ a3 w3 (b1 +b2 w+ )3 + a4 w4 (b1 + )4 ,
= (b1 w+b2 w2 +b3 w3 +b4 w4 ) + a2 w2 (b21 +2b1 b2 w+b22 w2 +2b1 b3 w2 + )
+ a3 w3 (b31 +3b21 b2 w+ ) + a4 w4 b41 + O(w5 ).
Comparing coefficients gives
w1 :
1 = b1 ,
w2 :
0 = b2 + a2 b21 ,
w3 :
0 = b3 + 2a2 b1 b2 + a3 b31 ,
w4 :
b2 = a2 ,
b3 = a3 + 2a22 ,
b4 = a4 + 5a2 a3 5a32 .
(B.5)
B.3
Lagranges theorem
(B.6)
with |f (z)| |z| and w z. The theorem of Lagrange states that in a suitable domain the
solution of this equation is
X (1)k d (k1)
z=w+
[f (w)]k
k=1
k!
dw
The proof of this theorem will be given in the last section of this appendix.
(B.7)
B.3
B.4
(B.8)
a2 w2 + a3 w3 + a4 w4 ,
[f (w)]2 =
w4 (a22 ) + w5 (2a2 a3 ) + ,
[f (w)]3 =
w6 (a32 ) + .
a2 w2 + a3 w3 + a4 w4 ,
1
D[f (w)]2 = 2w3 (a22 ) + 5w4 (a2 a3 ) + ,
2!
1 2
D [f (w)]3 = 5w4 (a32 ) + .
3!
Substitute in Lagranges expansion:
z = w f (w) +
1
1
D[f (w)]2 D2 [f (w)]3 + .
2!
3!
(B.9)
B.4
w=z+
(B.10)
(B.11)
p3 = b3 3b22 ,
p4 = b4 10b2 b3 + 15b32 .
(B.12)
Alternative notation
For the applications to cartography it is convenient to use the following notation for the
direct and inverse series:
b 2 2 b 3 3 b4 4
+ + ,
2!
3!
4!
p2 2 p3 3 p4 4
= z z z z +
2!
3!
4!
z=+
B.5
(B.13)
(B.14)
(B.15)
(B.16)
where the coefficients bn are small enough for the condition |f (z)| z, w to be valid;
note that we are assuming that w and z are of order unity. For the applications we have in
mind we shall have |bn | = O(en ) where e is the eccentricity of the ellipsoid. In deriving
the inversion we shall truncate the infinite Lagrange expansion at terms of order e8 ; thus we
retain terms proportional to b2 , b4 , b22 , b6 , b2 b4 , b32 , b8 , b2 b6 , b22 b4 , b24 , b42 and drop higher
powers.
B.5
In the following steps we make use of several trigonometric identities from Appendix C.
f (w) = b2 sin 2w + b4 sin 4w + b6 sin 6w + b8 sin 8w,
[f (w)]2 = b22 sin2 2w+2b2 b4 sin 2w sin 4w+b24 sin2 4w+2b2 b6 sin 2w sin 6w +
1
= b22 (1 cos 4w) + b2 b4 (cos 2w cos 6w)
2
1
+ b24 (1 cos 8w) + b2 b6 (cos 4w cos 8w) +
2
[f (w)]3 = b32 sin3 2w + 3b22 b4 sin2 2w sin 4w +
1
3
= b32 (3 sin 2w sin 6w) + b22 b4 (2 sin 4w sin 8w) +
4
4
1
[f (w)]4 = b42 sin4 2w + = b42 (3 4 cos 4w + cos 8w) + .
8
Calculate the derivatives
f (w) = b2 sin 2w + b4 sin 4w + b6 sin 6w + b8 sin 8w,
1
D[f (w)]2 = b22 sin 4w + b2 b4 ( sin 2w + 3 sin 6w)
2!
+ 2b24 sin 8w + 2b2 b6 ( sin 4w + 2 sin 8w) +
1 2
1
D [f (w)]3 = b32 ( sin 2w + 3 sin 6w) + 4b22 b4 ( sin 4w + 2 sin 8w) +
3!
2
1 3
4
D [f (w)]4 = b42 ( sin 4w + 2 sin 8w) + .
4!
3
Finally, substituting into
z = w f (w) +
1
1
1
D[f (w)]2 D2 [f (w)]3 + D3 [f (w)]4 +
2!
3!
4!
(B.17)
where
1
d2 = b2 b2 b4 + b32 ,
2
4
d4 = b4 + b22 2b2 b6 + 4b22 b4 b42 ,
3
3 3
d6 = b6 + 3b2 b4 b2 ,
2
8
d8 = b8 + 2b24 + 4b2 b6 8b22 b4 + b42 .
3
(B.18)
B.6
B.6
(B.19)
retaining only the terms up to w8 in the series for z. This problem is a trivial generalisation
of the derivation for the fourth order series developed in Section B.4, only the algebra is
a little more involved. We set f (w) = a2 w2 + a3 w3 in the Lagrange expansion and
evaluate the powers of [f (w)]k ; recall that we need retain only those powers of w which
give terms no higher than w8 after differentiating k 1 times.
f (w) =
[f (w)]2 =
a2 w2 + a3 w3 + a4 w4 + a5 w5 + a6 w6 + a7 w7 + a8 w8 ,
w4 (a22 )
+ w5 (2a2 a3 )
+ w6 (2a2 a4 + a23 )
+ w7 (2a2 a5 + 2a3 a4 )
+ w8 (2a2 a6 + 2a3 a5 + a24 )
+ w9 (2a2 a7 + 2a3 a6 + 2a4 a5 ) + O(w10 ),
[f (w)]3 =
w6 (a32 )
+ w7 (3a22 a3 )
+ w8 (3a22 a4 + 3a2 a23 )
+ w9 (3a22 a5 + 6a2 a3 a4 + a33 )
+ w10 (3a22 a6 + 6a2 a3 a5 + 3a2 a24 + 3a23 a4 ) + O(w11 ),
[f (w)]4 =
w8 (a42 )
+ w9 (4a32 a3 )
+ w10 (4a32 a4 + 6a22 a23 )
+ w11 (4a32 a5 + 12a22 a3 a4 + 4a2 a33 ) + O(w12 ),
[f (w)]5 =
w10 (a52 )
+ w11 (5a42 a3 )
+ w12 (5a42 a4 + 10a32 a23 ) + O(w13 ),
[f (w)]6 =
w12 (a62 )
[f (w)]8 =
O(w16 ).
B.7
z = w f (w) +
B.8
(B.21)
z = ww2 [a2 ]
w3 a3 2a22
w4 a4 5a2 a3 +5a32
w5 a5 3(2a2 a4 +a23 )+21a22 a3 14a42
w6 a6 7(a2 a5 +a3 a4 )+28(a22 a4 +a2 a23 )84a32 a3 +42a52
w7 a7 4(2a2 a6 +2a3 a5 +a24 )+12(3a22 a5 +6a2 a3 a4 +a33 )
30(4a32 a4 +6a22 a23 )+330a42 a3 132a62
w8 a8 9(a2 a7 +a3 a6 +a4 a5 )+15(3a22 a6 +6a2 a3 a5 +3a2 a24 +3a23 a4 )
165(a32 a5 +3a22 a3 a4 +a2 a33 )+99(5a42 a4 +10a32 a23 )
1287a52 a3 +429a72
B.7
(B.22)
w=z+
(B.23)
(B.24)
where
p2 = [b2 ]
p3 = b3 3b22
p4 = b4 10b2 b3 +15b32
p5 = b5 (15b2 b4 +10b23 )+105b22 b3 105b42
p6 = b6 (21b2 b5 +35b3 b4 )+(210b22 b4 +280b2 b23 )1260b32 b3 +945b52
p7 = b7 (28b2 b6 +56b3 b5 +35b24 )+(378b22 b5 +1260b2 b3 b4 +280b33 )
(3150b32 b4 +6300b22 b23 )+17325b42 b3 10395b62
p8 = [b8 (36b2 b7 +84b3 b6 +126b4 b5 )
+(630b22 b6 +2520b2 b3 b5 +1575b2 b24 +2100b23 b4 )
(6930b32 b5 +34650b22 b3 b4 +15400b2 b33 )
+(51975b42 b4 +138600b32 b23 )270270b52 b3 +135135b72
(B.25)
Comment: these results are extended to 12th order series in papers by (a) W E Bleieck and
(b) W G Bickley and J C P Miller. See bibliography.
B.9
B.8
In evaluating the inverse of the complex series that arises in the derivation of the transverse
Mercator projection on the ellipsoid (TME) we have the following coefficients
b2 = i s
b3 =
c2 W3
b4 = isc2 W4
b5 =
c4 W5
b6 = i sc4 W6
b7 = c6 W7
b8 = i sc6 W8
(B.26)
where i = 1, s = sin , c = cos , t = tan and the W functions are of the form
W 3 = t2
W4 = 4 2 + t2
W5 = 4 3 (1 6t2 ) + 2 (1 + 8t2 ) 2t2 + t4
W6 = 8 4 (11 24t2 ) 28 3 (1 6t2 ) + 2 (1 32t2 ) 2t2 + t4
W7 = 61 479t2 + 179t4 t6
W8 = 1385 3111t2 + 543t4 t6 ,
(B.27)
where is defined in equation (5.53). Substituting for the b-coefficients in (B.25) gives
p2 = ict [ 1 ]
p3 = c2 W3 +3t2
p4 = ic3 t W4 10W3 15t2
p5 = c4 W5 +(15t2 W4 10W32 )105t2 W3 105t4
p6 = ic5 t W6 (21W5 +35W3 W4 )(210t2 W4 280W32 ) +1260t2 W3 +945t4
h
p7 = c6 W7 +(28t2 W6 56W3 W5 +35t2 W42 ) (378t2 W5 +1260t2 W3 W4 280W33 )
(3150t4 W4 6300t2 W32 )+17325t4 W3 +10395t6
h
p8 = ic7 t W8 (36W7 +84W3 W6 +126W4 W5 )
(630t2 W6 2520W3 W5 +1575t2 W42 2100W32 W4 )
+(6930t2 W5 +34650t2 W3 W4 15400W33 )
+(51975t4 W4 138600t2 W32 ) 270270t4 W3 135135t6
(B.28)
Now substitute for the W . For p2 , . . . p6 we use the expressions given in (B.27). For p7 , p8
we use the spherical approximation (5.58) for all the terms on the right hand sides. That is
we set = 1 in W3 , . . . W6 on the right hand sides using the approximations
W3 W3 = 1 t2 ,
W4 W4 = 5 t2 ,
W5 W5 = 5 18t2 + t4 ,
W6 W6 = 61 58t2 + t4 .
(B.29)
B.10
p6 = ic5 t 8 4 (11 24t2 ) 28 3 (1 6t2 ) + 2 (1 32t2 ) 2t2 + t4
21 4 3 (1 6t2 ) + 2 (1 + 8t2 ) 2t2 + t4
35( t2 )(4 2 + t2 ) 210t2 (4 2 + t2 )
+ 280( t2 )2 + 1260t2 ( t2 ) +945t4
p7 = c6 61 479t2 + 179t4 t6 + 28t2 (61 58t2 + t4 ) 56(1 t2 )(5 18t2 + t4 )
+ 35t2 (5 t2 )2 378t2 (5 18t2 + t4 ) 1260t2 (1 t2 )(5 t2 )
+ 280(1 t2 )3 3150t4 (5 t2 ) + 6300t2 (1 t2 )2 + 17325t4 (1 t2 )
+10395t6
p8 = ic7 t 1385 3111t2 + 543t4 t6 36(61 479t2 + 179t4 t6 )
84(1 t2 )(61 58t2 + t4 ) 126(5 t2 )(5 18t2 + t4 )
630t2 (61 58t2 + t4 ) + 2520(1 t2 )(5 18t2 + t4 )
1575t2 (5 t2 )2 + 2100(1 t2 )2 (5 t2 ) + 6930t2 (5 18t2 + t4 )
+ 34650t2 (1 t2 )(5 t2 ) 15400(1 t2 )3 + 51975t4 (5 t2 )
138600t2 (1 t2 )2 270270t4 (1 t2 ) 135135t6
Note that we have changed p7 , p8 to p7 , p8 to show that these coefficients have been evaluated in the spherical approximation. Finally, simplifying these expressions gives
p2 = ict [ 1 ]
p3 = c2 + 2t2
p4 = ic3 t 4 2 9 6t2
p5 = c4 4 3 (1 6t2 ) 2 (9 68t2 ) 72t2 24t4
p6 = ic5 t 8 4 (11 24t2 ) 84 3 (3 8t2 ) + 225 2 (1 4t2 ) + 600t2 + 120t4
p7 = c6 61 + 662t2 + 1320t4 + 720t6
p8 = ic7 t 1385 7266t2 10920t4 5040t6
(B.30)
B.11
B.9
This derivation of the Lagrange expansion is included since it is to be found only in older
textbookssee bibliography. This account is based on a simplified version of that in Whittakers Modern Analysis (1902!!) where there is a more general statement of the theorem.
The derivation requires an excursion into complex analysis. In particular we require three
results which follow from Cauchys theorem. Since these results can be found in most texts
on complex analysis we quote them without proof.
Definition: a function f (z) is analytic in some domain D if it is single valued and
differentiable within D, except possibly at a finite number of points, the singularities of
f (z). If no point of D is a singularity then we say that f (z) is regular.
Cauchys integral formula: let f (z) be an analytic function, regular within a closed
contour C and continuous within and on C, and let a be a point within C. If in
addition f (z) has derivatives of all orders, then the n-th derivative at a is
I
f (z)
n!
dz.
(B.31)
f (n) (a) =
2i C (z a)n+1
The following result is usually found as a corollary to the proof of the principle of
the argument. If f (z) and g(z) are regular within and on a closed contour C and
f (z) has a simple zero at z = a then
I
g(z)f 0 (z)
1
dz.
(B.32)
g(a) =
2i C f (z)
Rouches theorem: if f (z) and g(z) are two functions regular within and on a closed
contour C, on which f (z) does not vanish and also |g(z)| < |f (z)|, then f (z) and
f (z) + g(z) have the same number of zeroes within C.
Let p(z) be regular within and on a closed contour C and let there be a single simple
zero at the point z = w inside C. Consider the equation
p(z) = t,
(B.33)
(B.34)
By Rouches theorem (with f p and g t) we see that p(z) and p(z)t have the same
number of zeroes inside C, namely one. The zero of p(z) is of course z = w: let the zero
of p(z) t be z = a. Therefore setting f (z) = p(z) t and g(z) = z in equation (B.32),
noting that t is a constant, we find the solution z = a of (B.33) is
I
1
zp0 (z)
z=a=
dz.
(B.35)
2i C p(z) t
B.12
I
C
"
n #
X
zp0 (z)
t
1+
dz.
p(z)
p(z)
(B.36)
Since |t| < |p(z)| on C the series is convergent and we can integrate term by term to find
a=
where
An t n ,
(B.37)
1
A0 =
2i
I
C
zp0 (z)
dz,
p(z)
1
An =
2i
I
C
zp0 (z)
dz
[p(z)]n+1
(n 1).
(B.38)
Now since p(z) has a simple zero at z = w the first integral may be integrated by setting
g(z) = z in equation (B.32).
A0 = w.
(B.39)
For the second integral we integrate by parts. The integral of the total derivative is zero
because the change in a single valued function around a closed curve is zero. Therefore
I
1 1
1
An =
dz, (n 1).
(B.40)
2i n C [p(z)]n
Now set
zw
.
r(z)
(B.41)
[r(z)]n
dz.
(z w)n
(B.42)
p(z) = (z w)q(z) =
so that
1 1
An =
2i n
I
C
p(z) has one zero inside C, at z = w, so q(z) will have no zeroes within C and r(z) will
have no poles within C. Using the Cauchy integral formula (B.31) An becomes
1
An = Dz(n1) [r(z)]n
n!
z=w
1 (n1)
= Dw
[r(w)]n
(n 1).
(B.43)
n!
Therefore the solution of
zw
=t
r(z)
(B.44)
is given by
z=a=w+
n
X
t
1
n!
(n1)
Dw
[r(w)]n .
(B.45)
B.13
Finally we set
f (z) = t r(z),
(B.46)
w = z + f (z),
(B.47)
X
(1)n
1
n!
(n1)
Dw
[f (w)]n .
(B.48)
This is the form of the expansion given in Section B.3. The domain of validity is discussed
in the textbooks. In the current applications we start from convergent series for w(z) and
find that the above series for z(w) is also convergent.
B.14
Appendix
Plane Trigonometry
A brief reminder of some identities from plane trigonometry which are required at various
points in the main text.
basic definition
(C.1)
(C.2)
(C.3)
y y
(C.4)
(C.5)
y y
(C.6)
(C.5) / (C.3)
(C.6) / (C.4)
(C.5) + (C.6)
xy
(C.3) + (C.4)
(C.4) (C.3)
x y x, y in (C.9)
x y x, y in (C.10)
x y x, y in (C.11)
x y x, y in (C.12)
1
[sin(x + y) + sin(x y)]
2
1
cos x sin y = [sin(x + y) sin(x y)]
2
1
cos x cos y = [cos(x + y) + cos(x y)]
2
1
sin x sin y = [cos(x y) cos(x + y)]
2
sin x cos y =
x+y
xy
cos
2
2
x+y
xy
sin x sin y = 2 cos
sin
2
2
x+y
xy
cos x + cos y = 2 cos
cos
2
2
x+y
xy
cos x cos y = 2 sin
sin
2
2
sin x + sin y =
2 sin
(C.7)
(C.8)
(C.9)
(C.10)
(C.11)
(C.12)
(C.13)
(C.14)
(C.15)
(C.16)
C.2
1 = cos2 x + sin2 x
y = x in (C.4)
sec2 x = 1 + tan2 x
2
(C.17)
(C.18)
(C.19)
(C.20)
use (C.17)
= 1 2 sin x
(C.21)
use (C.17)
= 2 cos2 x 1
(C.22)
y = x in (C.5)
from (C.21)
from (C.22)
use (C.21)
use (C.9)
from (C.26)
use (C.9)
from (C.24)
use (C.25)
1
sin x [1 cos 2x]
2
1
1
= sin x [sin 3x sin x]
2
4
3
1
= sin x sin 3x
4
4
(C.23)
(C.24)
(C.25)
sin3 x =
1
3
sin x cos x sin 3x cos x
4
4
3
1
= sin 2x [sin 4x + sin 2x]
8
8
1
= [2 sin 2x sin 4x]
8
(C.26)
sin3 x cos x =
sin4 x =
=
=
use (C.25)
cos3 x =
use (C.11)
=
=
1
[1 cos 2x]2
4
1
1 1
1 2 cos 2x + + cos 4x
4
2 2
1
[3 4 cos 2x + cos 4x]
8
1
cos x [1 + cos 2x]
2
1
1
cos x + [cos 3x + cos x]
2
4
3
1
cos x + cos 3x
4
4
(C.27)
(C.28)
(C.29)
C.3
Sk sin kx
(C.30)
Ck cos kx
(C.31)
NOTATION
sin2 x =
Hence
sin3 x =
sin4 x =
sin5 x =
sin6 x =
sin7 x =
sin8 x =
1
[1 C2 ]
2
1
[3S S3 ]
4
1
[3 4C2 + C4 ]
8
1
[10S 5S3 + S5 ]
16
1
[10 15C2 + 6C4 C6 ]
32
1
[35S 21S3 + 7S5 S7 ]
64
1
[35 56C2 + 28C4 8C6 + C8 ]
128
1
[S2 ]
2
1
sin3 x cos x = [2S2 S4 ]
8
1
sin5 x cos x =
[5S2 4S4 + S6 ]
32
1
sin7 x cos x =
[14S2 14S4 + 6S6 S8 ]
128
sinx cos x =
(C.32)
(C.33)
(C.34)
(C.35)
(C.36)
(C.37)
(C.38)
(C.39)
(C.40)
(C.41)
(C.42)
Hyperbolic functions
The basic definitions are
cosh x =
ex + ex
,
2
ex ex
.
2
(C.43)
eix eix
,
2i
(C.44)
sinh x=
eix + eix
,
2
sin x=
tan ix = i tanh x,
tanh ix = i tan x.
(C.45)
(C.46)
These identities can be used to derive all the hyperbolic formulae from the trigonometric
identities simply by replacing x and y by ix and iy. This effectively changes all cosine
C.4
terms to cosh. Each sine term becomes i sinh and where there is a single sinh in each term
of the identity an overall factor of i will cancel. Terms which have a product of two sines
will become a product of two i sinh terms giving an overall sign change. Likewise for the
tangent terms. We list only the identities corresponding to (C.3C.8) and (C.17C.23).
cosh(x y) = cosh x cosh y sinh x sinh y
(C.47)
(C.48)
1 = cosh2 x sinh2 x
(C.50)
(C.49)
sech x = 1 tanh x
(C.51)
cosech2 x = coth2 x 1
(C.52)
(C.53)
(C.54)
Appendix
Spherical trigonometry
D.1
Introduction
A great circle on a sphere is defined by the intersection of any plane through the centre of
the sphere with the surface of the sphere. Any two points on the sphere must lie on some
great circle and the shorter part of that great circle is also the shortest distance between
the points. In general three great circles define a spherical triangle (Figure D.1) and this
appendix develops the trigonometry of such triangles. There are many (old) text books and
we recommend Todhunters book on Spherical Trigonometrysee Bibliography.
Consider the three great circles defining the triangle ABC: they meet again in the points A1 , B1 and
C1 defining the triangle A1 B1 C1 . In fact they define
eight triangles since each pair of geodesics bounds
four triangles but ABC and A1 B1 C1 are counted
three times. (Think of slicing an apple into eight
pieces with three diametral cuts). Note that we do not
consider the improper triangles such as that formed
by the interior arcs BA , BC together with the exterior arc AC1 A1 C. Such improper triangles have one
angle greater than . Their solution presents no difficulty but we refer to Todhunters book for details.
Figure D.1
We now restrict attention to triangles in which the angles are less than or equal to .
The case of all equal to is degenerate for the triangle must then be three points on one
great circle with the sum of the angles equal to 3 and the sum of the sides equal to 2 (on
the unit sphere). The rigorous proof of this last statement is to be found in Euclid: Book 11,
Proposition 21see Bibliography. We shall see that it has as a corollary that the sum of the
angles of a spherical triangle is greater than . This lower bound is approached by small
triangles (sides much less than the radius) that are almost planar.
D.2
D.2
Figure D.2
Geometric proof
In Figure D.2 AD and AE are the tangents to the sides of the spherical triangle at A. As
long as the angles and are strictly less than /2 the tangent to the side AB meets the
radius OB extended to D and the tangent to the side AC meets the radius OC extended to
E. Since any tangent to the sphere is normal to the radius at the point of contact we have
that the triangles OAD and OAE are right angled.
We apply the planar cosine rule to the triangles ODE and ADE:
DE 2 = OD2 + OE 2 2OD.OE cos ,
DE 2 = AD2 + AE 2 2AD.AE cos A.
Subtracting these equations and using Pythagoras theorem to set OD2 AD2 = OA2 and
OE 2 AE 2 = OA2 we obtain
0 = 2OA2 + 2AD.AE cos A 2OD.OE cos
Dividing each term by the product OD.OE and using OA/OD = cos etc. gives
cos = cos cos + sin sin cos A.
It is conventional to express these identities in terms of the actual sides so that we should
set = a/R etc. If we assume that the lengths have been scaled to a unit sphere then the
above, alongwith the two relations obtained by cyclic permutations, becomes
cos a = cos b cos c + sin b sin c cos A,
cos b = cos c cos a + sin c sin a cos B,
cos c = cos a cos b + sin a sin b cos C.
(D.1)
D.3
For a small triangle with a, b, c 1 on the unit sphere, the spherical cosine rules
reduce to the planar cosine rule if we neglect cubic terms. For example the first becomes
a2
c2
b2
1
1
+ bc cos A,
= 1
2
2
2
which simplifies to
a2 = b2 + c2 2bc cos A.
The above proof assumes that the angles and are less than ninety degrees for the constructions as drawn. This restriction may be removed; it is discussed in detail in Todhunters
book (pages 16 to 19). The following alternative proof does not rely on these assumptions.
Cosine rule: vector proof
For any given spherical triangle we can introduce Cartesian axes with the z-axis along OA and the xz-plane defined by the plane OAB. Take the radius of the sphere
as unity and define vectors B and C along the radii OB
and OC respectively. The angle between the planes AOM
and AON is given by 6 M ON = A, so the components of
these unit vectors are
B = (sin c, 0 , cos c),
C = (sin b cos A, sin b sin A, cos b)
(D.2)
Figure D.3
Now the angle between the unit vectors is simply a, the angle subtended at the centre
by the arc BC. Therefore
BC = cos a = sin b sin c cos A + 0 + cos b cos c,
(D.3)
in agreement with our previous result for the cosine rule. This is the simplest proof of the
cosine rule: it needs no restrictions on the angles.
D.3
D.4
Therefore
sin A
= (a, b, c),
sin a
(D.4)
where
sin A
sin B
sin C
=
=
= (a, b, c)
sin a
sin b
sin c
(D.6)
ON 2 = OM 2 + M N 2 .
(D.7)
Figure D.4
D.5
D.4
In general, if we know three elements of a triangle then we might expect to find the other
three elements by direct application of the spherical sine and cosine rules. This is NOT
possible: to complete the solution in many cases we shall need further rules developed in
the ensuing sections.
The six distinct ways in which three elements may be given are shown in Figure D.5
along with a seventh case involving four given elements. In each figure the given elements
are shown below and the given angles are marked with a small arc and the given sides
are marked with a cross bar; each figure has variations given by cyclic permutations. The
solution of such spherical triangles is harder than in the planar case because we do not know
the sum of the angles: given two angles we do not know the third.
Figure D.5
D.6
This is an appropriate point to mention that any determination of an angle or side from
its sine will generally lead to ambiguities since sin x = sin( x). However the angles
and sides on the unit sphere are in the interval (0, ) so their determination from cosines,
secants, tangents or cotangents will be unambiguous. Likewise the sine, cosine or tangent
of any half-angle (or side) is positive and its inverse is also unambiguous. Many of the
formulae that we will derive were established to avoid the sine ambiguity.
D.5
Figure D.6a below shows the three great circles which intersect to form the spherical triangle
ABC. In addition we show the normals to the plane of each great circle; each intersects
the sphere in two points each of which is a pole when a specific great circle is identified
as an equator. As shown some of the poles (small solid circles) are visible and some(open
circles) are on the hidden face. Three of these six poles may be used to define the polar
triangle. The convention is that A and its pole A0 lie on the same side of the diametral plane
containing BC; likewise for the others. We shall now prove the following statements.
The sides of the polar triangle A0 B 0 C 0 are the supplements of the angles of the original triangle ABC. (We assume a unit sphere on which the lengths of the sides are
equal to the radian measure of the angles they subtend at the centre).
The angles of the polar triangle A0 B 0 C 0 are the supplements of the sides of the original triangle ABC.
Figure D.6
Figure D.6b shows the triangle ABC, the pole A0 of A and the corresponding equator
formed by extending the side BC. Note three properties:
1. Any great circle through the pole A0 to its equator BC is a quadrant arc of length /2
(on the unit sphere), i.e. A0 K = A0 L = /2.
2. Any great circle through A0 intersects its equator BC at a right angle, as at K and L.
3. The angle (in radians) between two such quadrant arcs is equal to the length of the
segment cut on the equator by the arcs, i.e. = 6 KA0 L = KL.
D.7
Figure D.7, which is neither an elevation nor a perspective view, shows the schematic relation between the triangle ABC and its polar triangle A0 B 0 C 0 . The sides of
ABC are extended along their great circles to meet the
sides of A0 B 0 C 0 at the points shown. From the the three
properties discussed in the previous paragraph we can
deduce the following results.
Figure D.7
The great circles A0 B 0 and A0 C 0 through the pole A0 intersect the equator corresponding to A0 , that is BC extended, at points D and E. The intersections are right
angles and the distance DE is equal to the angle A0 expressed in radians. Therefore
A0 = DE.
Now consider the intersections of the great circle B 0 C 0 with the great circles defined
by AB and AC. Since the angles at F and G are right angles we deduce that A must
be the pole to the equator B 0 C 0 . Similarly B, C must be the poles of the equators
C 0 A0 and A0 B 0 repectively. We conclude that the polar triangle of the polar triangle
A0 B 0 C 0 must be the original triangle ABC. Consequently (1) since C is the pole
of A0 B 0 we must have CD = /2; (2) since B is the pole of C 0 A0 we must have
BE = /2; (3) since A is the pole of B 0 C 0 we must have F G = A.
A0 = DE = DC + BE BC =
Similar results follow for the other angles and sides of the polar triangle so that:
A0 = a
a0 = A
B0 = b
b0 = B
C 0 = c,
c0 = C.
(D.8)
An important corollary follows from the existence of the polar triangle. We have already
stated that Euclid proves that the sum of the sides of a spherical triangle on the unit sphere
satisfies = a+b+c < 2. Applying this to the polar triangle gives 3ABC < 2
so that , the sum of the angles, is greater than . Since we conventionally take the angles
to be less than then we must have < < 3. (The restriction to angles and sides less
than may be lifted; the so-called improper triangles so formed are discussed in Todhunter.
We have no need to consider them here.)
D.8
(D.9)
Now substitute for angles and sides usung equation (D.8) noting that cos( ) = cos
and sin( ) = sin :
cos A + cos B cos C = sin B sin C cos a,
cos B + cos C cos A = sin C sin A cos b,
cos C + cos A cos B = sin A sin B cos c.
(D.10)
Now these equations, obtained by applying a known rule to the polar triangle, are obviously
new relations between the elements of the original triangle; they are called the supplemental
cosine rules. This is an example of a powerful method of generating a new formula from
any that we have already found.
The supplemental cosine rules clearly provide a way of solving a spherical triangle when
all three angles are given. This is Case 6 in Figure D.5.
Note that a new rule does not always arise. For example, applying the sine rule to
A0 B 0 C 0 gives
sin A0
sin B 0
sin C 0
=
=
.
sin a0
sin b0
sin c0
On substituting (D.8) we have the usual rules simply inverted:
sin b
sin c
sin a
=
=
.
sin A
sin B
sin C
Alternative derivation of the supplemental cosine rules
It is possible to derive the supplemental cosine rules directly without appealing to the polar
triangle. For example, in the first formula of (D.10) substitute for the terms on the left-hand
side using the normal cosine rules:
(cos a cos b cos c) sin2 a + (cos b cos c cos a)(cos c cos a cos b)
sin2 a sin b sin c
2
2
cos a[1 cos a cos b cos2 c + 2 cos a cos b cos c]
=
sin2 a sin b sin c
2
= cos a sin b sin c
D.9
D.6
The six elements of a triangle may be written in an anti-clockwise order as (aCbAcB). The
cotangent, or four-part, formulae relate two sides and two angles forming four consecutive
elements around the triangle, for example (aCbA) or BaCb. The six distinct formulae that
we shall prove are
cos b cos C = cot a sin b cot A sin C,
cos b cos A = cot c sin b cot C sin A,
cos c cos A = cot b sin c cot B sin A,
cos c cos B = cot a sin c cot A sin B,
cos a cos B = cot c sin a cot C sin B,
cos a cos C = cot b sin a cot B sin C,
(a)
(b)
(c)
(d)
(e)
(f )
(aCbA)
(CbAc)
(bAcB)
(AcBa)
(cBaC)
(BaCb),
(D.11)
where the subset of elements involved is shown to the right of every equation. In the first
equation, for the set aCbA, we term a and A the outer elements and C and b the inner
elements. With this notation the general form of the equations is
cos(inner side). cos(inner angle) =
(D.12)
Note that the inner elements of each set formula occur twice and cannot be deduced from
the other elements; only the outer elements of each set may be derived in terms of the
other three. For example in the first equation involving the set aCbA we can only determine
the outer side a in terms of CbA or the outer angle A in terms of aCb. Note also that the
outer angle or side is determined from its cotangent so that there is no ambiguity.
To prove the first formula start from the cosine rule (D.1a) and on the right-hand side
substitute for cos c from (D.1c) and for sin c from (D.6):
cos a = cos b cos c + sin b sin c cos A
= cos b (cos a cos b + sin a sin b cos C) + sin b sin C sin a cot A
2
cos a sin b = cos b sin a sin b cos C + sin b sin C sin a cot A.
The result follows on dividing by sin a sin b. Similar techniques with the other two cosine
rules give D.11c,e. Equations D.11b,d,f follow by applying D.11e,a,c to the polar triangle.
Solution of spherical triangles II
The four-part formulae may be used to give solutions to two of the cases discussed
in Section D.4. In Case 2 in Figure D.5, where we are given (bAc), we can use equations D.11b,c to find the angles C, B from their cotangents: we can then find a from D.11a
without any sine ambiguity. We can now solve Case 4, where we are given (BaC), by using
equations D.11e,f to give the sides c, b and we can then find A from D.11a. We are still left
with the problem solving Case 7 (since Cases 3, 5 can also be reduced to Case 7).
D.10
D.7
1/2
a
cos S cos(SA)
sin =
sin B sin C
2
1/2
a
cos(SB) cos(SC)
cos =
sin B sin C
2
1/2
a
cos S cos(SA)
tan =
cos(SB) cos(SC)
2
(D.13)
To prove the first formula use cos A = 1 2 sin2 (A/2) and the cosine rule (D.1).
sin2
A
1 cos A
=
2
2
1 cos a cos b cos c
cos(bc) cos a
=
2
2 sin b sin c
2 sin b sin c
1
a+bc
ab+c
=
sin
sin
.
sin b sin c
2
2
D.11
cos 21 (ab)
sin 21 (AB)
cos 12 c
cos 21 C
cos 12 (a+b)
cos 12 (AB)
cos 21 c
sin 12 C
sin 12 (ab)
sin 12 (a+b)
sin 12 c
(D.14)
cos 12 (A+B)
sin 21 C
sin 21 c
These are proved by expanding the numerator on the left hand side and using the half angle
formulae. For example, using equations C.5, C.13 and C.23
1
A
B
A
B
sin (A+B) = sin cos + cos sin
2
2
2
2
2
1/2
1/2
2
sin s sin (sb) sin(sc)
sin s sin2 (sa) sin(sc)
=
+
sin a sin b sin2 c
sin a sin b sin2 c
1/2
sin(sb) + sin(sa) sin s sin(sc)
=
sin c
sin a sin b
1
1
sin 2 c cos 2 (ab)
1
=
cos C,
1
1
2
sin 2 c cos 2 c
and hence the required result.
Napiers analogies
Published by Napier in 1614. His methods were purely geometric but we obtain them by
dividing the Delambre formulae.
tan 12 (A+B) =
tan
1
2 (AB)
cos 12 (ab)
cos
1
2 (a+b)
sin 12 (ab)
sin
1
2 (a+b)
cot 12 C
cot
1
2C
tan 21 (a+b) =
tan
1
2 (ab)
cos 12 (AB)
cos
1
2 (A+B)
sin 12 (AB)
sin
1
2 (A+B)
tan 12 c
tan 21 c
(D.15)
D.12
D.8
Right-angled triangles
There are many problems in which one of the angles, say C, is equal to /2. In this case
there are only 5 elements and in general two will suffice to solve the triangle. We shall
show that the solution of such a triangle can be presented as a set of 10 equations involving
3 elements so that every element can be expressed in terms of any pair of the other elements.
The required 10 equations involving C are found from the third cosine rule (D.1), two
sine rules (D.7), four cotangent formulae (D.11) and all three of the supplemental cosine
rules (D.10). Setting C = /2 we obtain (from the equations indicated)
(D.1c)
(D.11b)
(D.7b)
(D.11e)
(D.7c)
(D.10a)
(D.11a)
(D.10b)
(D.11f)
(D.10c)
(D.16)
As an example suppose we are given a and c (and C = /2). Then we can find c, A, B
from the first, fourth and fifth equations.
Figure D.8
(D.17)
For example if we take 21 c as the middle part the first rule gives sin(/2 c) =
tan(/2 A) tan(/2 B) which gives the last of the equations in (D.16); if we apply the
second rule we get sin(/2 c) = cos a cos b which is the first of the equations in (D.16).
D.13
D.9
Quadrantal triangles
The triangle ABC is quadrantal if at least one side subtends an angle of /2 at the centre of
the sphere. Without loss of generality take c = /2. Therefore the angle C 0 = c of the
polar triangle is equal to /2. Now apply Napiers rules to the polar triangle with C 0 = /2
and
a0 = A, b0 = B, A0 = a, B 0 = b.
The circular parts of the polar triangle
a0 ,
b0 ,
A0 ,
c0 ,
B0,
2
2
2
must be replaced by
A,
B,
a ,
C ,
b ,
2
2
2
Noting that sin(x/2) = cos x, cos(x/2) =
sin x and tan(x/2) = cot x we have the following
equations:
Figure D.9
sin A =
sin a sin C,
sin B =
sin b sin C,
cos a =
sin b cos A,
tan A =
tan a sin B,
cos b =
sin a cos B,
tan B =
tan b sin A,
Example
As an example of a quadrantal triangle we consider
a problem arising in the discussion of geodesics on a
sphere in Chapter 11. With the following identifications
a = s,
b = ,
c= ,
2
2
A = ,
B = 0 ,
C = . (D.19)
the equations (D.18) become
Figure D.10
(D.18)
(D.20)
The practical problems are (a) given 0 and s find , and ; (b) given , find 0 and s.
For the first we use the fourth, ninth, and then the first equation. For the second we use the
fifth and eighth equations.
D.14
Figure D.11
Appendix
(z b) 0
(z b)2 00
(z b)3 000
f (b) +
f (b) +
f (b) + .
1!
2!
3!
(E.1)
or, alternatively,
z 0
z2
z3
f (b) + f 00 (b) + f 000 (b) + .
1!
2!
3!
When b = 0 we obtain Maclaurins series.
z2
z3
z
f (z) = f (0) + f 0 (0) + f 00 (0) + f 000 (0) + .
1!
2!
3!
f (b + z) = f (b) +
E.2
(E.2)
(E.3)
z2
z3
z4
sin b
cos b +
sin b +
2!
3!
4!
(E.4)
z2
z3
z4
cos b +
sin b +
cos b +
2!
3!
4!
(E.5)
8
+z = 1 + 2z + 2z 2 + z 3 +
4
3
z3
(1+3 tan2 b) sec2 b+
3
(E.6)
(E.7)
b
1
z2
+
+
(E.8)
1/2
2
2 (1 b2 )3/2
(1 b )
1
z 2 2b
z3
2
8b2
arctan(b+z) = arctan b + z
+
+
1+b2
2! (1+b2 )2
3! (1+b2 )2 (1+b2 )3
z4
24b
48b3
+
+
(E.9)
4! (1+b2 )3 (1+b2 )4
arcsin(b+z) = arcsin b + z
/over
E.2
E.3
Logarithms
1
1
1
z z2 + z3 z4 +
1 < z 1
2
3
4
1
1
1
1 z < 1
ln(1 z) = z z 2 z 3 z 4 +
2
3
4
1+z
2
2
2
ln
= 2z + z 3 + z 5 + z 7 +
|z| < 1
1z
3
5
7
Trigonometric
1
1
1
sin z = z z 3 + z 5 z 7 +
|z| <
3!
5!
7!
1
1
1
|z| <
cos z = 1 z 2 + z 4 z 6 +
2!
4!
6!
1
2
17 7
tan z = z + z 3 + z 5 +
z +
|z| <
3
15
315
2
1 2
5 4
61 6
sec z = 1 + z + z +
z +
|z| <
2
24
720
2
7 3
31 5
1 1
z +
z + 0 < |z| <
csc z = + z +
z 6
360
15120
1 1
1
2 5
cot z = z z 3
z
0 < |z| <
z 3
45
945
Inverse trig
1
3
5 7
arcsin z = z + z 3 + z 5 +
z +
|z| < 1
6
40
112
1
1
1
1
arctan z = z z 3 + z 5 z 7 + z 9
|z| < 1
3
5
7
9
Hyperbolic
1
1
1
sinh z = z + z 3 + z 5 + z 7 +
|z| <
3!
5!
7!
1
1
1
cosh z = 1 + z 2 + z 4 + z 6 +
|z| <
2!
4!
6!
1
2
17 7
tanh z = z z 3 + z 5
z +
|z| <
3
15
315
2
1 2
5 4
61 6
sech z = 1 z + z
z
|z| <
2
24
720
2
ln(1 + z) =
(E.10)
(E.11)
(E.12)
(E.13)
(E.14)
(E.15)
(E.16)
(E.17)
(E.18)
(E.19)
(E.20)
(E.21)
(E.22)
(E.23)
(E.24)
/over
E.4
E.3
Setting f (z) = (1 + z)n , n an integer, in the Maclaurin series gives the standard binomial
series:
(1 + z)n = 1 + nz +
n(n1) 2
n!
z + +
zr + ,
2!
(nr)! r!
(E.25)
(1 + z)1 = 1 z + z 2 z 3 + z 4 ,
(E.26)
(1 + z)2 = 1 2z + 3z 2 4z 3 + 5z 4 .
(E.27)
(E.28)
1
3
5
35 4
(1 + z)1/2 = 1 z + z 2 z 3 +
z ,
2
8
16
128
(E.29)
3
15
35
315 4
(1 + z)3/2 = 1 z + z 2 z 3 +
z ,
2
8
16
128
(E.30)
E.4
Appendix
Calculus of variations
The simplest problem in the calculus of variations is as follows. Let F (x, y, y 0 ) be a function
of x and some unspecified function y(x) and also its derivative. For every y(x) we construct
the following integral between fixed points A and B at which x = a and x = b:
Z
J[y] =
F (x, y, y 0 )dx.
(F.1)
The problem is to find the particular function y(x) which, for a given function F (x, y, y 0 ),
minimises or maximises J[y]. In general we will not be able to say that we have a maximum
or a minimum solution but the context of any particular problem will usually decide the
matter. The following method only guarantees that J[y] will be extremal. The solution here
is valid for twice continuously differentiable functions.
Figure F.1
We first tighten our notation a little. We assume that an extremal function can be found
and that it is denoted by y(x), the heavy path C in the figure; J[y] then refers to the value of
the integral on the extremal path. We consider the set of all paths AB defined by functions
y where
y(x) = y(x) + (x),
(F.2)
where (x) is an arbitrary function such that (a) = (b) = 0, thus guaranteeing that
the end points of all paths are the same. One of these paths is denoted C1 in the figure.
The set of integrals for one given (x) and varying may be considered as generating a
F.2
F (x, y + , y 0 + 0 )dx.
(F.3)
In this notation the value of the integral on the extremal path is (0) and the condition that
it is an extremum is
d
= 0.
(F.4)
d =0
The Taylor series for the integrand is
F (x, y + , y 0 + 0 ) = F + Fy + Fy0 0 + O(2 )
(F.5)
y0
dx
a
a
a
Since the first term vanishes we have proved that for an extremal
Z b
(x) H(x) dx = 0,
(F.7)
(F.8)
d
Fy 0 F y .
(F.9)
dx
We now show that equation (F.8) implies that H(x) = 0. This result rejoices under the
grand name of the fundamental lemma of the calculus of variations. The proof is by
contradiction: first suppose that H(x) 6= 0, say positive, at some point x0 in (a, b). Then
there must be an interval (x1 , x2 ) surrounding x0 in which H(x) > 0. Since (x) can
be any suitably differentiable function we take = (x2 x)4 (x x1 )4 in [x1 , x2 ] and
Rb
zero elsewhere. Clearly, for such a function we must have a H dx > 0, in contradiction
to (F.8). Therefore our hypothesis that H 6= 0 is not valid. Therefore we must have H = 0,
giving the EulerLagrange equations:
d F
F
E ULER L AGRANGE
= 0.
(F.10)
0
dx y
y
where
H(x) =
In this equation the partial derivatives indicate merely the formal operations of differentiating F (x, y, y 0 ) with respect to y and y 0 as if they were independent variables. On the other
hand the operator d/dx is a regular derivative and the above equation expands to
2F
2F 0 2F 00 F
+
y + 02 y
= 0.
0
xy
yy 0
y
y
(F.11)
This is a second order ordinary differential equation for y(x): it has a solution with two
arbitrary constants which must be fitted at the end points.
F.3
(F.13)
(F.14)
(F.15)
Extensions
There are many variants of the above results:
1. F has two (or more) dependent functions: F x, u(x), u0 (x), v(x), v 0 (x)
2. F has two (or more) independent variables: F x, y, u(x, y), u0 (x, y)
3. Both of above: F x, y, u(x, y), u0 (x, y), v(x, y), v 0 (x, y)
4. F involves higher derivatives: F (x, y, y 0 , y 00 , . . .).
5. The end points are not held fixed.
Only the first of the above concerns us here. The proof is along the same lines as above but
we need to make two independent variations and set
u
(x) = u(x) + (u) (x),
v(x) = v(x) + (v) (x).
Equation (F.8) now becomes of the form
Z bh
i
(u) (x) H (u) (x) + (v) (x) H (v) (x) dx = 0.
a
(F.16)
F.4
Since (u) and (v) are arbitrary independent functions we obtain H (u) = H (v) = 0, i.e.
d F
F
= 0,
(F.17)
0
dx u
u
F
d F
= 0.
0
dx v
v
(F.18)
Sufficiency
The EulerLagrange equations have been shown to be a necessary conditions for the existence of an extremal integral. The proof of sufficiency is non-trivial and is discussed in
advanced texts.
Appendix
Complex numbers
A complex number z is a pair of real numbers, x, y combined with the basic imaginary
number i in the expression z = x + iy. Such complex numbers may be manipulated just
as real numbers with the proviso that i2 = 1. We say that x is the real part of the complex
number, x = Re (z), and y is the imaginary part, y = Im (z). From z = x + iy we form its
complex conjugate z = x iy. Note that zz = (x + iy)(x iy) = x2 + y 2 . The complex
number z = x + iy may be be represented as a point (x, y) in a plane which is called the
complex z-plane. It is also useful to introduce polar coordinates in the plane and write
z = x + iy = r(cos + i sin ).
(G.1)
Figure G.1
The simplest complex function we can consider is a finite polynomial such as:
w(z) = 3 + z + z 2 .
If we substitute z = x + iy in this expression, using i2 = 1, we obtain
(
u(x, y) = 3 + x + x2 y 2 ,
w(z) = u(x, y) + iv(x, y)
where
v(x, y) = y + 2xy.
(G.3)
(G.4)
Here we have written w(z) in terms of two real functions of two variables. All complex functions can be split up in this way.
G.2
(G.5)
where the coefficients will, in general, be complex numbers. The real and imaginary
parts of w(z) will be infinite series.
1
1
1
1
1
z + z2 + z3 + z4 + z5 + .
1!
2!
3!
4!
5!
(G.6)
It can be proved that this series is convergent for all values of z. Note that when z is
purely real, z = x, the series reduces to the usual real definition of exp(x). When z
is purely imaginary, z = i say, we get a very interesting result:
exp i = 1 +
i
1
i
1
i
2 3 + 4 5 + .
1!
2!
3!
4!
5!
(G.7)
Now the real terms in this expansion are simply those in the expansion of cos , whilst
the imaginary terms are those that arise in the expansion of sin . Therefore we can
write the polar coordinate expression of z in (G.1) as
z = r(cos + i sin ) = r exp(i) = rei .
(G.8)
(G.9)
We can also define the sine and cosine functions of a complex number by the series
1 2
z +
2!
1
sin z = z z 3 +
3!
cos z = 1
1 4
z ,
4!
1 5
z .
5!
(G.10)
(G.11)
Once again, when z is purely real, z = x, these series reduce to the usual real definition of cos(x) and sin(x). When z is purely imaginary, z = i say, we get the
following important results using (E.21,E.22):
cos i = cosh ,
(G.12)
sin i = i sinh ,
(G.13)
tan i = i tanh .
(G.14)
For a general complex number we can use compound angle formulae (C.3,C.4) to
determine the real and imaginary parts of the sine and cosine functions:
cos(x + iy) = cos x cos iy sin x sin iy = cos x cosh y i sin x sinh y,
(G.15)
sin(x + iy) = sin x cos iy + cos x sin iy = sin x cosh y + i cos x sinh y.
(G.16)
G.3
In Chapter 4 we require the real and imaginary parts of cot z; these follow from the
quotient of the last two equations. Simply multiply top and bottom by the complex
conjugate of the denominator and use the identities for cos 2x, cosh 2x etc. given in
Appendix C. It will be convenient to use different notation for this example.
cot(A + iB) =
sin 2A i sinh 2B
cosh 2B cos 2A
Similarly, we define the hyperbolic sine and cosine functions of a complex numbers.
When z is purely real, z = x, these series reduce to the usual real definition of cosh x
and sinh x. When z is purely imaginary, z = i say, we get the counterparts of
(G.12G.14)
1 2
z +
2!
1
sinh z = z + z 3 +
3!
cosh i = cos ,
cosh z = 1 +
G.2
(G.17)
1 4
z + ,
4!
1 5
z + ,
5!
(G.18)
(G.19)
(G.20)
sinh i = i sin ,
(G.21)
tanh i = i tan .
(G.22)
f 0 (x) = lim
(G.23)
The small print of the definition is that the limit when x tends to zero from above (x
0+) should be equal to the limit when x tends to zero from below (x 0). In principle
these limits could be different and we would then have to define two different derivatives,
say f+0 (x) and f0 (x). A simple example where the limits differ is the function f (x) = |x|,
for which f+0 = +1 and f0 = 1 at the origin. The only point we wish to make is that
even in one dimension we must be careful about directions when defining derivatives.
G.4
f
f
x +
y.
x
y
(G.25)
If the direction of the displacement is taken in the direction of the unit vector n and the
magnitude of the displacement is s, then we can set x = nx s and y = ny s. We can
then define a directional derivative in two dimensions as
df
f
f
= nx
+ ny .
(G.26)
ds n
x
y
The point we wish to stress is that the derivatives of functions of two variables are essentially
dependent on direction.
w0 (z) = lim
(G.27)
The crucial step is that we demand that this limit should exist independent of the direction
in which z tends to zero. If such a limit exists in all points of some region of the complex
plane then we say that w(z) is an analytic (or regular) function of z (in that region). This
restriction on differentiation is very strong and as a result analytic functions are very special,
with many interesting properties.
G.5
(G.29)
1
(uy + ivy ) = vy iuy .
i
(G.30)
If we now demand that these two derivatives w0 (z) are the same we have the Cauchy
Riemann equations:
ux = vy ,
uy = vx
(G.31)
It can also be shown that if the partial derivatives ux etc. are continuous, then the Cauchy
Riemann conditions are sufficient for the derivative w0 (z) to exist.
(G.32)
the real and imaginary parts and their partial derivatives are
u = x3 3xy 2 ,
v = 3x2 y y 3 .
ux = 3x2 3y 2 ,
vx = 6xy,
uy = 6xy,
vy = 3x2 3y 2 .
These equations show that the CauchyRiemann equations (G.31) are indeed satisfied
and we can use either (G.29) or (G.30) to identify the derivative as
w0 (z) = ux + ivx = vy iuy
= 3x2 3y 2 + i6xy = 3(x + iy)2
= 3z 2 .
(G.33)
G.6
Similarly we can prove that w(z) = z n is analytic with a derivative given by w0 (z) =
nz n1 . (The proof is easier if the CauchyRiemann conditions are written in terms
of polar coordinates and z n is written as rn ein .
Consider the function w(z) = sin z. The real and imaginary parts of the sine function
were determined in equation (G.16) so that w(z) = u + iv where
u = sin x cosh y,
v=
cos x sinh y.
ux = cos x cosh y,
vx = sin x sinh y,
uy = sin x sinh y,
vy =
cos x cosh y.
Once again the CauchyRiemann equations are indeed satisfied and we can identify
the derivative from equation (G.15):
w0 (z) = ux + ivx = vy iuy
= cos x cosh y i sin x sinh y
= cos z.
(G.34)
In similar ways we can show that all the derivatives of standard functions parallel
those that arise for functions of one real variable.
Taylors theorem
We state without proof or qualification that under reasonable conditions an analytic function may be represented by a convergent Taylors series. In the following development we
shall use the theorem in the following form.
w(z) = w(z0 )+
1
1
1
(zz0 )w0 (z0 )+ (zz0 )2 w00 (z0 )+ (zz0 )3 w000 (z0 )+ . (G.35)
1!
2!
3!
It is also true that any convergent power series defines an analytic function. Proofs of these
statements are to be found in the standard texts on complex functions.
/continued overleaf
G.7
G.3
Mathematicians and geographers both use the term map in essentially the same way. A
complex function w(z) may be viewed as simply a pair of real functions which define a
map in the sense that it takes a point (x, y) in the complex z-plane into a point (u, v)
in the complex w-plane by virtue of the two functions u(x, y) and v(x, y). Points go to
points, regions go to regions, curves through a point go to curves through the image point,
(Figure G.2). The important result is that if w(z) is an analytic function then , the angle
of intersection of two curves C1 and C2 at P , is equal to 0 the angle of intersection of
the image curves at the image point. Such maps are said to be be conformal. We proceed
immediately to the proof of this statement.
Figure G.2
Let z0 be a fixed point on the curve C in the z-plane. Let z be a nearby point on C and write
z z0 = rei . Note that is the angle between real axis and the chord; in the limit z z0 ,
this angle will approach the angle between the real axis and the tangent to C at z0 . Let w0
Figure G.3
and w be the corresponding image points and set w w0 = r0 ei . Taylors theorem tells
us that
1
1
w(z) = w(z0 ) + (z z0 )w0 (z0 ) + (z z0 )2 w00 (z0 ) + ,
(G.36)
1!
2!
so that we can write
w w0
1
= w0 (z0 ) + (z z0 )w00 (z0 ) + .
(G.37)
z z0
2!
G.8
The derivative of the function w(z) at z0 is a unique complex number which depends only
on the position z0 and we can write it as A(z0 ) exp(i(z0 )) where A(z0 ) and are both
real.
Therefore in the limit as z z0 equation (G.37) becomes
0
0
r
lim
exp i( ) = w0 (z0 ) = A(z0 ) exp(i(z0 )),
zz0
r
(G.38)
since the remaining terms on the RHS vanish in the limit. We deduce that
0
r
lim
= A(z0 ),
zz0
r
exp i(00 0 ) = exp[i(z0 )],
(G.39)
(G.40)
where 0 and 00 are the angles between the tangents and the real axes. Note that the second of these equations can be derived only when A 6= 0. The value of 0 becomes
indeterminate if A = 0 so we must therefore demand that w0 (z0 ) 6= 0.
The second of the above limits, when it exists, shows that 00 = 0 + (z0 ), that is the
tangent at P is rotated by an angle when it is mapped to the w-plane. This will be true
of all curves through P and consequently the angle of intersection of any two curves will be
preserved under the mapping. This is the definition of a conformal mapping.
For a given measurement accuracy we can always find an infinitesimal region around P
in which the variation of A and is imperceptible. Equations (G.39,G.40) then imply that
the small region is scaled and rigidly rotated, so preserving its shape. This is the property
of orthomorphism.
References
Geodesy
In addition to the Wiki pages we select only two.
The Ordnance Survey of Great Britain (OSGB) have an excellent web site on which
can be found A guide to coordinate systems in Great Britain. This article discusses
modern approaches to reference systems. It can be found at:
http://www.ordnancesurvey.co.uk/oswebsite/gps/information/index.html
This article mentions geodesy only briefly in the introduction. Those who wish to go (much)
deeper can consult the following texts:
R.2
M ALING D H,
(1992), Coordinate Systems and Map Projections,
Pergamon, ISBN : 0-08-037234-1.
The following is a comprehensive summary of just about all projections in use, including
of course all the Mercator projections: normal, transverse and oblique. The introductions
to the projections are very readable and moreover each is complemented with an historical
survey. BUT there are no derivations of any projection formulae.
S NYDER J P, (1987),
Map Projections: a Working Manual,
US Geological Survey, Professional Paper 1395
Published by US Government Printing Office but also available on the web at:
http://pubs.er.usgs.gov/usgspubs/pp/pp1395
For a more advanced treatment covering all projections (but with many gaps in the mathematical development requiring large amounts of work to fill in):
The article referred to in the geodesy section, A guide to coordinate systems in Great
Britain, also includes a statement of the transverse Mercator projection formulae (without derivations) with examples of transforming between grid and geographical coordinates.
There is also an excellent (short) Wikipedia page entitled British National Grid Reference
System which gives a direct link to the above article.
R.3
journal was entitled the Empire Survey Review. Note also that this journal is published in
parts (from four to eight a year) and bound in two year periods so both volume and part
number are specified. Sadly this journal is not readily available in general libraries.
Two important papers are listed below. The first, by Lee, is the first article in the journal
to present a correct derivation of the Transverse Mercator projection formulae. The second
paper, by Redfearn, presents a derivation of the series to high enough order to be applicable
to all practical problems. By one of those quirks of fate it is Redfearns name that has entered the literature. The set of papers by Hotine present a much less elegant derivation which
refuses to countenance the existence of complex variable methods: they are not discussed
here.
H OTINE, M (1946, 1948), Survey Review, Vol 8, Part 62 pp 301311 and Vol 9,
Part 63 pp 2935, Part 64 pp 5270, Part 65 pp 112123, Part 66 pp 157166.
The orthomorphic projection of the spheroid, parts IV.
It should be said at once that these papers are fairly terse when it comes to the derivations
of the projection formulae. The present article errs in the other direction and it is fair to say
that no extra details will be found in the original papers. Note also that this article uses a
different convention for the names of axes; basically x and y are exchanged.
The discussion of geodesic problems in Chapter 11 includes a greatly expanded version of
the following paper:
Mathematics
The appendices include all the mathematics we require and a little more besides. They are
derived from first principles and should hopefully not require further background reading.
In their preparation I found that modern texts were not helpful on the whole because they
were too distant from application. The older books were much more useful. A few texts are
listed here.
Spherical Trigonometry
R.4
Differential Geometry
Lagrange Expansions
The derivation of the Lagrange expansions has essentially disappeared from modern texts.
The proof in Appendix B is a combination of
The twelfth order Lagrange series are published in the Philosophical Magazine.
Index
Index
Africa, size, 2.8
Airy(1830) ellipsoid, 1.3
analytic functions, G.4
azimuth, 1.6, 2.6, 3.9, 3.10, 6.2, 8.2
Bond, Henry, 1.9
calculus of variations, F.1
cartography, 1.4, R.2
CauchyRiemann conditions, 3.12, 4.7, 4.13,
7.6, 8.1, G.5
Clarke(1866) ellipsoid, 1.2, 1.3
complex variables, 4.44.13, 8.4, G.1
conformal latitude, 5.4, 6.8, 6.9
conformal projection, 1.5
conformality
analytic functions, G.7
conditions, see CauchyRiemann
scale isotropy, 3.12
convergence, see grid convergence
curvature, 5.7, A.1A.9
curvature quotient , 5.9
datum, 1.2
GRS80, 1.2
ID1830, 1.3
NAD27, 1.3
NAD83, 1.3
OSGB36, 1.3
WGS72, 1.2
degreeradian conversion, 2.1
distance on the sphere, 2.2
double projections, 6.7
eastings and northings, 9.4
eccentricity, 5.2
ellipsoid
Airy1830, 1.2, 2.2, 5.2
Clarke1866, 1.2
curvature, 5.7, A.3
X.2
Index X. Index
Index X. Index
scale factor
NGGB and UTM, 9.7
Snyder J P, R.2
South America, size, 2.8
spherical limit, 5.9, 7.4
spherical trigonometry, 3.2, D.1D.14, R.3
spheroid, 1.1
Survey Review, R.2
topographic surveying, 1.2
true north, 3.10
true origin, 9.2
Universal Transverse Mercator, see UTM
UTM, 1.11, 3.1, 9.1, R.1
Vincent, T, R.3
Vincenty formulae, 1.11, 5.11
WGS(1972) ellipsoid, 1.2
Wright, Edward, 1.8
X.3
X.4
Index X. Index