coordinate system - An example Tissot ellipse for an Equirectangular Projection?

I am trying to calculate distortion of an equirectangular projection via Tissot indicatrices. I have tried following the directions on this post (among other things) but it is beyond my understanding as an amateur.

So, I am wondering if someone would be so kind as to calculate a single Tissot ellipse for an example equirectangular lat/long (whichever is your favorite and is distorted on an equirectangular projection). I don't understand exactly what the variables are and where they come from, so seeing the equations in action would be very useful.

I am trying to understand these equations so I can plug them into a mapping program I am coding. I asked a bunch of general questions in this thread, but I think a specific example will help me figure the rest out.

Many thanks, as always.

NCashew

Answer

For the record, here is a complete, commented implementation of the Tissot indicatrix (and related) calculations in R, with a worked example. The source of the equations is John Snyder's Map Projections--A Working Manual.

tissot <- function(lambda, phi, prj=function(z) z+0, asDegrees=TRUE, A = 6378137, f.inv=298.257223563, ...) {
  #

  # Compute properties of scale distortion and Tissot's indicatrix at location `x` = c(`lambda`, `phi`)
  # using `prj` as the projection.  `A` is the ellipsoidal semi-major axis (in meters) and `f.inv` is
  # the inverse flattening.  The projection must return a vector (x, y) when given a vector (lambda, phi).
  # (Not vectorized.)  Optional arguments `...` are passed to `prj`.
  #
  # Source: Snyder pp 20-26 (WGS 84 defaults for the ellipsoidal parameters).
  # All input and output angles are in degrees.
  #
  to.degrees <- function(x) x * 180 / pi
  to.radians <- function(x) x * pi / 180

  clamp <- function(x) min(max(x, -1), 1)                             # Avoids invalid args to asin
  norm <- function(x) sqrt(sum(x*x))
  #
  # Precomputation.
  #
  if (f.inv==0) {                                                     # Use f.inv==0 to indicate a sphere
    e2 <- 0 
  } else {
    e2 <- (2 - 1/f.inv) / f.inv                                       # Squared eccentricity
  }

  if (asDegrees) phi.r <- to.radians(phi) else phi.r <- phi
  cos.phi <- cos(phi.r)                                               # Convenience term
  e2.sinphi <- 1 - e2 * sin(phi.r)^2                                  # Convenience term
  e2.sinphi2 <- sqrt(e2.sinphi)                                       # Convenience term
  if (asDegrees) units <- 180 / pi else units <- 1                    # Angle measurement units per radian
  #
  # Lengths (the metric).
  #
  radius.meridian <- A * (1 - e2) / e2.sinphi2^3                      # (4-18)
  length.meridian <- radius.meridian                                  # (4-19)

  radius.normal <- A / e2.sinphi2                                     # (4-20)
  length.normal <- radius.normal * cos.phi                            # (4-21)
  #
  # The projection and its first derivatives, normalized to unit lengths.
  #
  x <- c(lambda, phi)
  d <- numericDeriv(quote(prj(x, ...)), theta="x")
  z <- d[1:2]                                                         # Projected coordinates
  names(z) <- c("x", "y")
  g <- attr(d, "gradient")                                            # First derivative (matrix)

  g <- g %*% diag(units / c(length.normal, length.meridian))          # Unit derivatives
  dimnames(g) <- list(c("x", "y"), c("lambda", "phi"))
  g.det <- det(g)                                                     # Equivalent to (4-15)
  #
  # Computation.
  #
  h <- norm(g[, "phi"])                                               # (4-27)
  k <- norm(g[, "lambda"])                                            # (4-28)
  a.p <- sqrt(max(0, h^2 + k^2 + 2 * g.det))                          # (4-12) (intermediate)
  b.p <- sqrt(max(0, h^2 + k^2 - 2 * g.det))                          # (4-13) (intermediate)

  a <- (a.p + b.p)/2                                                  # (4-12a)
  b <- (a.p - b.p)/2                                                  # (4-13a)
  omega <- 2 * asin(clamp(b.p / a.p))                                 # (4-1a)
  theta.p <- asin(clamp(g.det / (h * k)))                             # (4-14)
  conv <- (atan2(g["y", "phi"], g["x","phi"]) + pi / 2) %% (2 * pi) - pi # Middle of p. 21
  #
  # The indicatrix itself.
  # `svd` essentially redoes the preceding computation of `h`, `k`, and `theta.p`.
  #
  m <- svd(g)

  axes <- zapsmall(diag(m$d) %*% apply(m$v, 1, function(x) x / norm(x)))
  dimnames(axes) <- list(c("major", "minor"), NULL)

  return(list(location=c(lambda, phi), projected=z, 
           meridian_radius=radius.meridian, meridian_length=length.meridian,
           normal_radius=radius.normal, normal_length=length.normal,
           scale.meridian=h, scale.parallel=k, scale.area=g.det, max.scale=a, min.scale=b, 
           to.degrees(zapsmall(c(angle_deformation=omega, convergence=conv, intersection_angle=theta.p))),
           axes=axes, derivatives=g))
}

indicatrix <- function(x, scale=1, ...) {
  # Reprocesses the output of `tissot` into convenient geometrical data.
  o <- x$projected
  base <- ellipse(o, matrix(c(1,0,0,1), 2), scale=scale, ...)             # A reference circle
  outline <- ellipse(o, x$axes, scale=scale, ...)
  axis.major <- rbind(o + scale * x$axes[1, ], o - scale * x$axes[1, ])
  axis.minor <- rbind(o + scale * x$axes[2, ], o - scale * x$axes[2, ])
  d.lambda <- rbind(o + scale * x$derivatives[, "lambda"], o - scale * x$derivatives[, "lambda"])
  d.phi <- rbind(o + scale * x$derivatives[, "phi"], o - scale * x$derivatives[, "phi"])
  return(list(center=x$projected, base=base, outline=outline, 

              axis.major=axis.major, axis.minor=axis.minor,
              d.lambda=d.lambda, d.phi=d.phi))
}
ellipse <- function(center, axes, scale=1, n=36, from=0, to=2*pi) {
  # Vector representation of an ellipse at `center` with axes in the *rows* of `axes`.
  # Returns an `n` by 2 array of points, one per row.
  theta <- seq(from=from, to=to, length.out=n)
  t((scale * t(axes))  %*% rbind(cos(theta), sin(theta)) + center)
}
#

# Example: analyzing a GDAL reprojection.
#
library(rgdal)

prj <- function(z, proj.in, proj.out) {
  z.pt <- SpatialPoints(coords=matrix(z, ncol=2), proj4string=proj.in)
  w.pt <- spTransform(z.pt, CRS=proj.out)
  return(w.pt@coords[1, ])
}
r <- tissot(130, 54, prj,                # Longitude, latitude, and reprojection function

       proj.in=CRS("+init=epsg:4267"),   # NAD 27
       proj.out=CRS("+init=esri:54030")) # World Robinson projection

i <- indicatrix(r, scale=10^4, n=71)
plot(i$outline, type="n", asp=1, xlab="Easting", ylab="Northing")
polygon(i$base, col=rgb(0, 0, 0, .025), border="Gray")
lines(i$d.lambda, lwd=2, col="Gray", lty=2)
lines(i$d.phi, lwd=2, col=rgb(.25, .7, .25), lty=2)
lines(i$axis.major, lwd=2, col=rgb(.25, .25, .7))
lines(i$axis.minor, lwd=2, col=rgb(.7, .25, .25))

lines(i$outline, asp=1, lwd=2)