diff --git a/classes.cpp b/classes.cpp
index 29952d89..a6cb83c7 100644
--- a/classes.cpp
+++ b/classes.cpp
@@ -973,8 +973,8 @@ enum eModel : int {
   // 20..24
   mdEquivolume, mdCentralInversion, mdSimulatedPerspective, mdTwoHybrid, mdGeodesic,
   // 25
-  mdMollweide, mdCentralCyl, mdCollignon, mdHorocyclic,
-  // 29..
+  mdMollweide, mdCentralCyl, mdCollignon, mdHorocyclic, mdQuadrant, mdAxial, mdAntiAxial,
+  // 32..
   mdGUARD, mdPixel, mdHyperboloidFlat, mdPolynomial, mdManual
   };
 #endif
@@ -1016,6 +1016,9 @@ EX vector<modelinfo> mdinf = {
   {X3("central cylindrical"), mf::euc_boring | mf::band, DEFAULTS},
   {X3("Collignon"), mf::pseudoband | mf::equiarea, DEFAULTS},
   {X3("horocyclic coordinates"), mf::euc_boring, DEFAULTS},
+  {X3("quadrant coordinates"), mf::euc_boring, DEFAULTS},
+  {X3("axial coordinates"), mf::euc_boring, DEFAULTS},
+  {X3("anti-axial coordinates"), mf::euc_boring, DEFAULTS},
   {X3("guard"), 0, DEFAULTS},
   {X3("polynomial"), mf::conformal, DEFAULTS},
   };
diff --git a/drawing.cpp b/drawing.cpp
index a4a33a2c..2cde2b00 100644
--- a/drawing.cpp
+++ b/drawing.cpp
@@ -262,14 +262,19 @@ void add1(const hyperpoint& H) {
   glcoords.push_back(glhr::pointtogl(H)); 
   }  
 
+int axial_sign() {
+  return ((axial_x ^ axial_y)&1) ? -1:1;
+  }
+
 bool is_behind(const hyperpoint& H) {
+  if(pmodel == mdAxial && sphere) return axial_sign() * H[2] <= BEHIND_LIMIT;
   return pmodel == mdDisk && (hyperbolic ? H[2] >= 0 : true) && (nonisotropic ? false : pconf.alpha + H[2] <= BEHIND_LIMIT);
   }
 
 hyperpoint be_just_on_view(const hyperpoint& H1, const hyperpoint &H2) {
   // H1[2] * t + H2[2] * (1-t) == BEHIND_LIMIT - pconf.alpha
   // H2[2]- BEHIND_LIMIT + pconf.alpha = t * (H2[2] - H1[2])
-  ld t = (H2[2] - BEHIND_LIMIT + pconf.alpha) / (H2[2] - H1[2]);
+  ld t = (axial_sign() * H2[2] - BEHIND_LIMIT + (pmodel == mdAxial ? 0 : pconf.alpha)) / (H2[2] - H1[2]) * axial_sign();
   return H1 * t + H2 * (1-t);
   }
 
@@ -2136,6 +2141,16 @@ EX void draw_main() {
       }
     #endif
     }
+  else if(pmodel == mdAxial && sphere) {
+    for(auto& ptd: ptds) if(ptd->prio == PPR::OUTCIRCLE)
+      ptd->draw();
+    for(axial_x=-4; axial_x<=4; axial_x++) 
+    for(axial_y=-4; axial_y<=4; axial_y++) 
+    for(auto& ptd: ptds) if(ptd->prio != PPR::OUTCIRCLE) {
+      ptd->draw();
+      }
+    glflush();
+    }
   else {
     DEBB(DF_GRAPH, ("draw_main1"));
     if(ray::in_use && !ray::comparison_mode) {
diff --git a/hypgraph.cpp b/hypgraph.cpp
index 108a1ba9..2a7925c5 100644
--- a/hypgraph.cpp
+++ b/hypgraph.cpp
@@ -343,6 +343,8 @@ hyperpoint compute_hybrid(hyperpoint H, int rootid) {
 
 EX ld signed_sqrt(ld x) { return x > 0 ? sqrt(x) : -sqrt(-x); }
 
+EX int axial_x, axial_y;
+
 EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
 
   hyperpoint H = H_orig.h;
@@ -486,6 +488,73 @@ EX void applymodel(shiftpoint H_orig, hyperpoint& ret) {
       break;
       }
     
+    case mdAxial: {
+      models::apply_orientation_yz(H[1], H[2]);
+      models::apply_orientation(H[0], H[1]);
+      
+      ld& mt = pconf.model_transition;
+      
+      ld z = H[LDIM];
+      if(mt != 1) z += (1-mt) * pconf.alpha;
+
+      ret[0] = H[0] / z;
+      ret[1] = H[1] / z;
+      if(GDIM == 3) ret[2] = H[2] / z;
+      else ret[2] = 0;
+      ret[3] = 1;
+      
+      if(mt) for(int i=0; i<LDIM; i++)  {
+        if(mt < 1) 
+          ret[i] *= mt;
+        ret[i] = atan_auto(ret[i]);
+        if(mt < 1) 
+          ret[i] /= mt;
+        }
+      
+      if(sphere) ret[0] += axial_x * M_PI, ret[1] += axial_y * M_PI;
+
+      models::apply_orientation(ret[1], ret[0]);
+      models::apply_orientation_yz(ret[2], ret[1]);
+      break;
+      }
+    
+    case mdAntiAxial: {
+      models::apply_orientation_yz(H[1], H[2]);
+      models::apply_orientation(H[0], H[1]);
+      
+      ret[0] = asin_auto(H[0]);
+      ret[1] = asin_auto(H[1]);
+
+      ret[2] = 0; ret[3] = 1;
+
+      models::apply_orientation(ret[1], ret[0]);
+      models::apply_orientation_yz(ret[2], ret[1]);
+      break;
+      }
+    
+    case mdQuadrant: {
+      H = space_to_perspective(H);
+      models::apply_orientation_yz(H[1], H[2]);
+      models::apply_orientation(H[0], H[1]);
+      
+      tie(H[0], H[1]) = make_pair((H[0] + H[1]) / sqrt(2), (H[1] - H[0]) / sqrt(2));
+
+      H[1] += 1;
+      double rad = sqhypot_d(GDIM, H);
+      H /= -rad;
+      H[1] += .5;
+      
+      H *= 2;
+      
+      ld x = exp(-H[0]/2);
+      ret[0] = -H[1] * x - 1;
+      ret[1] = H[1] / x + 1;
+
+      models::apply_orientation(ret[1], ret[0]);
+      models::apply_orientation_yz(ret[2], ret[1]);
+      break;
+      }
+    
     case mdHorocyclic: {
 
       find_zlev(H);
diff --git a/models.cpp b/models.cpp
index 3ae281cd..50097d35 100644
--- a/models.cpp
+++ b/models.cpp
@@ -205,7 +205,7 @@ EX namespace models {
   
   EX bool has_orientation(eModel m) {
     return
-      among(m, mdHalfplane, mdPolynomial, mdPolygonal, mdTwoPoint, mdJoukowsky, mdJoukowskyInverted, mdSpiral, mdSimulatedPerspective, mdTwoHybrid, mdHorocyclic) || mdBandAny();
+      among(m, mdHalfplane, mdPolynomial, mdPolygonal, mdTwoPoint, mdJoukowsky, mdJoukowskyInverted, mdSpiral, mdSimulatedPerspective, mdTwoHybrid, mdHorocyclic, mdAxial, mdAntiAxial, mdQuadrant) || mdBandAny();
     }
   
   EX bool is_perspective(eModel m) {
@@ -218,7 +218,7 @@ EX namespace models {
     }
   
   EX bool has_transition(eModel m) {
-    return among(m, mdJoukowsky, mdJoukowskyInverted, mdBand) && GDIM == 2;
+    return among(m, mdJoukowsky, mdJoukowskyInverted, mdBand, mdAxial) && GDIM == 2;
     }
   
   EX bool product_model(eModel m) {