Mollweide projection

classes.cpp

@@ -607,8 +607,8 @@ const modelinfo mdinf[int(mdPolynomial)+1] = {
   {X3("central inversion"), mf::azimuthal | mf::conformal},
   {X3("two-point azimuthal"), mf::euc_boring | mf::twopoint},
   {X3("two-point hybrid"), mf::euc_boring | mf::twopoint},
-  {X3(""), 0},
-  {X3(""), 0},
+  {X3("geodesic"), 0},
+  {X3("Mollweide"), mf::euc_boring | mf::quasiband | mf::equiarea },
   {X3(""), 0},
   {X3("polynomial"), mf::conformal}
   };

classes.h

@@ -285,7 +285,9 @@ enum eModel {
   mdRotatedHyperboles, mdSpiral, mdPerspective,
   // 20..24
   mdEquivolume, mdCentralInversion, mdSimulatedPerspective, mdTwoHybrid, mdGeodesic,
-  // 25..
+  // 25
+  mdMollweide,
+  // 26..
   mdGUARD, mdUnchanged, mdHyperboloidFlat, mdPolynomial, mdRug, mdFlatten
   };
 

drawing.cpp

@@ -686,20 +686,34 @@ int mercator_coord;
 int mercator_loop_min = 0, mercator_loop_max = 0;
 ld mercator_period;
 
+ld period_at(ld y) {
+  
+  ld m = current_display->radius;
+  y /= (m * vid.stretch);
+  
+  switch(pmodel) {
+    case mdBand:
+      return m * 4;
+    case mdSinusoidal:
+      return m * 2 * cos(y * M_PI);
+    case mdMollweide:
+      return (abs(y) > .5-1e-6) ? m * 2 : m * 2 * sqrt(1 - y*y*4);
+    default:
+      return m * 2;
+    }     
+  }
+
 void fixMercator(bool tinf) {
 
-  if(pmodel == mdBand)
-    mercator_period = 4 * current_display->radius;
-  else
-    mercator_period = 2 * current_display->radius;
+  mercator_period = period_at(0);
   
   if(!models::model_straight)
     for(auto& g: glcoords)
       models::apply_orientation(g[0], g[1]);
     
-  if(pmodel == mdSinusoidal)
+  if(among(pmodel, mdSinusoidal, mdMollweide))
     for(int i = 0; i<isize(glcoords); i++) 
-      glcoords[i][mercator_coord] /= cos(glcoords[i][1] / current_display->radius / vid.stretch * M_PI);
+      glcoords[i][mercator_coord] *= mercator_period / period_at(glcoords[i][1-mercator_coord]);
 
   ld hperiod = mercator_period / 2;
     
@@ -720,10 +734,10 @@ void fixMercator(bool tinf) {
   if(pmodel == mdBandEquiarea)
     dmin = -vid.stretch * current_display->radius / M_PI, dmax = vid.stretch * current_display->radius / M_PI;
 
-  for(int i = 0; i<isize(glcoords); i++) {
-    while(glcoords[0][mercator_coord] < hperiod) glcoords[0][mercator_coord] += mercator_period;
-    while(glcoords[0][mercator_coord] > hperiod) glcoords[0][mercator_coord] -= mercator_period;
-    }
+ // for(int i = 0; i<isize(glcoords); i++) {
+  while(glcoords[0][mercator_coord] < hperiod) glcoords[0][mercator_coord] += mercator_period;
+  while(glcoords[0][mercator_coord] > hperiod) glcoords[0][mercator_coord] -= mercator_period;
+ //   }
     
   ld first = glcoords[0][mercator_coord];
   ld next = first;
@@ -754,9 +768,9 @@ void fixMercator(bool tinf) {
       mercator_loop_min--, mincoord -= mercator_period;
     while(maxcoord < cmax)
       mercator_loop_max++, maxcoord += mercator_period;
-    if(pmodel == mdSinusoidal)
+    if(among(pmodel, mdSinusoidal, mdMollweide))
       for(int i = 0; i<isize(glcoords); i++) 
-        glcoords[i][mercator_coord] *= cos(glcoords[i][1] / current_display->radius / vid.stretch * M_PI);    
+        glcoords[i][mercator_coord] *= period_at(glcoords[i][1-mercator_coord]) / mercator_period;
     if(!models::model_straight)
       for(auto& g: glcoords)
         models::apply_orientation(g[1], g[0]);
@@ -784,9 +798,10 @@ void fixMercator(bool tinf) {
         }
       minto += mercator_period;
       }
-    if(pmodel == mdSinusoidal)
+    if(among(pmodel, mdSinusoidal, mdMollweide))
       for(int i = 0; i<isize(glcoords); i++) 
-        glcoords[i][mercator_coord] *= cos(glcoords[i][1] / current_display->radius / vid.stretch * M_PI);
+        glcoords[i][mercator_coord] *= period_at(glcoords[i][1-mercator_coord]) / mercator_period;
+
     glcoords.push_back(glcoords.back());
     glcoords.push_back(glcoords[0]);
     for(int u=1; u<=2; u++) {
@@ -1092,10 +1107,10 @@ void dqi_poly::draw() {
   
     if(l || lastl) { 
       for(int i=0; i<isize(glcoords); i++) {
-        if(pmodel == mdSinusoidal) {
+        if(among(pmodel, mdSinusoidal, mdMollweide)) {
           ld y = glcoords[i][1], x = glcoords[i][0];
           models::apply_orientation(x, y);
-          mercator_period = 2 * current_display->radius * cos(y / current_display->radius / vid.stretch * M_PI);
+          mercator_period = period_at(y);
           }
         glcoords[i][mercator_coord] += models::ocos * mercator_period * (l - lastl);
         glcoords[i][1-mercator_coord] += models::osin * mercator_period * (l - lastl);

hypgraph.cpp

@@ -656,6 +656,27 @@ EX void applymodel(hyperpoint H, hyperpoint& ret) {
       makeband(H, ret, make_twopoint);
       break;
     
+    case mdMollweide: 
+      makeband(H, ret, [] (ld& x, ld& y) { 
+        ld theta = 
+          hyperbolic ? min(y / 2 + 0.572365, y * 0.78509) :
+          euclid ? y :
+          y > 0 ? max(y * 0.012/0.015, M_PI/2 - (M_PI/2-y) * 0.066262/0.015708) :
+          min(y * 0.012/0.015, -M_PI/2 + (M_PI/2+y) * 0.066262/0.015708);
+          
+        if(sphere && abs(theta) >= M_PI/2 - 1e-6) ;
+        else {
+          for(int it=0; it<4; it++) {
+            auto a = (sin_auto(2*theta) +2*theta - M_PI * sin_auto(y));
+            auto b = (2 + 2 * cos_auto(2*theta));
+            theta = theta - a / b;
+            // theta = theta - (sinh(2*theta) +2*K*theta - M_PI * sinh(y)) / (2 * K + 2 * cosh(2*theta));
+          } }
+        y = M_PI * sin_auto(theta) / 2;
+        x = x * cos_auto(theta); 
+        });
+      break;
+    
     case mdBandEquiarea: 
       makeband(H, ret, [] (ld& x, ld& y) { y = sin_auto(y); });
       break;
@@ -1629,7 +1650,7 @@ EX void draw_boundary(int w) {
   dynamicval<ld> dw(vid.linewidth, vid.linewidth * (svg::in ? svg::divby : 1));
   #endif
 
-  if(elliptic && !among(pmodel, mdBand, mdBandEquidistant, mdBandEquiarea, mdSinusoidal))
+  if(elliptic && !among(pmodel, mdBand, mdBandEquidistant, mdBandEquiarea, mdSinusoidal, mdMollweide))
     circle_around_center(M_PI/2, periodcolor, 0, PPR::CIRCLE);
   
   switch(pmodel) {
@@ -1658,7 +1679,7 @@ EX void draw_boundary(int w) {
       return;
       }
     
-    case mdBand: case mdBandEquidistant: case mdBandEquiarea: case mdSinusoidal: {
+    case mdBand: case mdBandEquidistant: case mdBandEquiarea: case mdSinusoidal: case mdMollweide: {
       if(DIM == 3) return;
       if(pmodel == mdBand && models::model_transition != 1) return;
       bool bndband = ((pmodel == mdBand) ? hyperbolic : sphere);

models.cpp

@@ -108,7 +108,7 @@ EX namespace polygonal {
 
 #if HDR
 inline bool mdAzimuthalEqui() { return among(pmodel, mdEquidistant, mdEquiarea, mdEquivolume); }
-inline bool mdBandAny() { return among(pmodel, mdBand, mdBandEquidistant, mdBandEquiarea, mdSinusoidal); }
+inline bool mdBandAny() { return among(pmodel, mdBand, mdBandEquidistant, mdBandEquiarea, mdSinusoidal, mdMollweide); }
 #endif
 
 EX namespace models {