Convert lat, long to x, y on Mollweide

Question

I have tried to follow the instructions here but I get wild results compared to this site.

Here is my code.

#include <cmath>

double solveNR(double latitude, double epsilon) {
  if (abs(latitude) == M_PI / 2) {
    return latitude;
  }

  double theta = latitude;
  while (true) {
    double nextTheta = theta - (2 * theta * std::sin(2 * theta) - M_PI * std::sin(latitude)) / (2 + 2 * std::cos(2 * theta));
    if (abs(theta - nextTheta) < epsilon) {
      break;
    }
    theta = nextTheta;
  }
  return theta;
}

void convertToXY(double radius, double latitude, double longitude, double* x, double* y) {

  latitude = latitude * M_PI / 180;
  longitude = longitude * M_PI / 180;
  double longitudeZero = 0 * M_PI / 180;
  double theta = solveNR(latitude, 1);

  *x = radius * 2 * sqrt(2) * (longitude - longitudeZero) * std::cos(theta) / M_PI;
  *y = radius * sqrt(2) * std::sin(theta);
}

For instance,

180 longitude = 21

90 latitude = 8.1209e+06

assuming a radius of 5742340.81

I found this resource which seems to calculate the right answer. But I cannot parse how it is different.

Answer 1

In your solveNR() function why do you use

double nextTheta = theta - (2 * theta * std::sin(2 * theta) - PI * 
std::sin(latitude)) / (2 + 2 * std::cos(2 * theta));

instead

double nextTheta = theta - (2 * theta + std::sin(2 * theta) - PI * 
std::sin(latitude)) / (2 + 2 * std::cos(2 * theta));

Seems like you should use "+" instead "*" (after 2 * theta in the numerator), to accord with wiki-instructions.