[GRASSLIST:3244] help with theory of map projections

[chris2]

I have been trying to learn the theory behind map projections from the book "Map Projections, A reference Manual" by Bugayevskiy and Snyder.  I am stuck on page 2 in Chapter 1, if you can believe it, where they take an ellipsoid of revolution, and a point on its surface at latitude phi, and mention that the radius of curvature of the meridian through the point is:

M= a*(1-e^2) / [(1-e^2*sin(phi)^2)^(3/2)]

Where e, a, and b are the eccentricity and semimajor and semiminor axes of the ellipse of revolution.

I set out to verify this formula as a way of reviewing my vector caclulus and getting started with the book.

I used the following parameterization of the ellipse:

x = a*cos(t)
y = b*sin(t)

My t, by the way, is equivalent to their phi, I am pretty sure.  That is, t is the angle between a normal line to the ellipse through the given point, and the x-axis.

And, using the notation of my calculus book (Thomas and Finney), went through the following formulas:

R = a*cos(t) i + b*sin(t) j

V = dR/dt
T = dR/ds = V / |V|
k = curvature = | dT/ds | = | dT/dt | / | ds/dt | = |dT/dt| / |V|
radius of curvature = 1/k

I couldn't reduce the expression I got to the one from the book, so I wrote a little R script to graph the radius of curvature from 0 to 2*pi using Bugayevskiy's formula and the one I derived. I have included my R code below to produce the plots.

Can anyone spot what I am doing wrong?  The black graph is Bugayevskiy's and the red one is mine.

Also, can anyone recommend a book with a little more explanation than "Map Projections"?  My ultimate goal is to understand the UTM coordinate system well enough to write code to do projections with it.  B and S don't seem to feel the need to derive their equations. 

I have included my R code below.  If you stuck it in a file "curve.R", you could produce the plots by typing "source("curve.R")" from the R prompt.

Chris Marshall

pow <- function(x,y){
   exp(y*log(x))
}

a <- 3
b <- 0.5
e <- sqrt((a*a-b*b))/a
t <- (1:100)*2*pi/100
x <- a*cos(t)
y <- b*sin(t)

# B and S curvature formulas
k1 <- e*sin(t)
M <- a*(1-e*e)/pow(1-k1*k1,3/2)
N <- a/sqrt(1-k1*k1)

# my brute force approach
st <- sin(t)
s2t <- st*st
ct <- cos(t)
c2t <- ct*ct
abs.v <- sqrt(a*a*s2t+b*b*c2t)
j <- 1/abs.v
dj.dt <- (b*b-a*a)*ct*st/(abs.v*abs.v*abs.v)
dg.dt <- -a*ct*j - a*st*dj.dt
dh.dt <- -b*st*j + b*ct*dj.dt
abs.dT.dt <- sqrt(dg.dt*dg.dt+dh.dt*dh.dt)
k <- abs.dT.dt/abs.v

# graph the curvature formulas
graphics.off()
ylim <- c(0,max(M,N,1/k))
#plot(x,y,type="l")
#X11()
plot(t,M,type="l",ylim=ylim)
#X11()
#plot(t,N,type="l",ylim=ylim)
#X11()
#plot(t,1/k,type="l",ylim=ylim)
lines(t,1/k,type="l",col=2)