diff --git a/geom-exp.cpp b/geom-exp.cpp
index 37e13213..21008d0b 100644
--- a/geom-exp.cpp
+++ b/geom-exp.cpp
@@ -1256,6 +1256,11 @@ int read_geom_args() {
     set_variation(eVariation::bch);
     shift(); reg3::subcube_count = argi();
     }
+  else if(argis("-coxeter")) {
+    PHASEFROM(2);
+    set_variation(eVariation::coxeter);
+    shift(); reg3::coxeter_param = argi();
+    }
   #endif
   #if CAP_FIELD
   else if(argis("-fi")) {
diff --git a/geometry.cpp b/geometry.cpp
index 179069bb..65f49f9e 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -507,6 +507,10 @@ EX bool is_subcube_based(eVariation var) {
   return among(var, eVariation::subcubes, eVariation::dual_subcubes, eVariation::bch);
   }
 
+EX bool is_reg3_variation(eVariation var) {
+  return var == eVariation::coxeter;
+  }
+
 void geometry_information::prepare_basics() {
 
   DEBBI(DF_INIT | DF_POLY | DF_GEOM, ("prepare_basics"));
diff --git a/reg3.cpp b/reg3.cpp
index 535bdfeb..3e7e0a89 100644
--- a/reg3.cpp
+++ b/reg3.cpp
@@ -17,6 +17,12 @@ namespace hr {
 EX namespace reg3 {
 
   EX int subcube_count = 1;
+  
+  EX flagtype coxeter_param = 0;
+
+  const flagtype cox_othercell = 1;
+  const flagtype cox_midedges = 2;
+  const flagtype cox_vertices = 4;
 
   #if HDR
   inline short& altdist(heptagon *h) { return h->emeraldval; }
@@ -301,6 +307,70 @@ EX namespace reg3 {
       }
     }
   
+  EX void generate_coxeter(flagtype f) {
+    auto& ssh = cgi.subshapes;
+    for(auto& fac: cgi.cellshape) {
+      hyperpoint facectr = Hypc;
+      vector<hyperpoint> ring;
+      hyperpoint last = fac.back();
+      ring.push_back(last);
+      for(hyperpoint h: fac) {
+        if(f & cox_midedges)
+          ring.push_back(mid(last, h));
+        ring.push_back(last = h);
+        facectr += h;
+        }
+      facectr = normalize(facectr);
+
+      hyperpoint fc2 = rspintox(facectr) * xpush0(2*hdist0(facectr));
+      
+      if(f & cox_vertices) {
+        for(int i=1; i<isize(ring); i++) {
+          ssh.emplace_back();
+          auto &ss = ssh.back();
+          auto h = (f & cox_othercell) ? facectr : fc2;
+          ss.faces.push_back({C0, h, ring[i-1]});
+          ss.faces.push_back({C0, h, ring[i]});
+          ss.faces.push_back({C0, ring[i-1], ring[i]});
+          ss.faces.push_back({h, ring[i-1], ring[i]});
+          }
+        }
+      else if(f & cox_midedges) {
+        ring.push_back(ring[1]);
+        for(int i=3; i<isize(ring); i+=2) {
+          ssh.emplace_back();
+          auto &ss = ssh.back();
+          auto h = (f & cox_othercell) ? facectr : fc2;
+          ss.faces.push_back({C0, ring[i-2], ring[i-1]});
+          ss.faces.push_back({C0, ring[i-1], ring[i-0]});
+          ss.faces.push_back({C0, h, ring[i-2]});
+          ss.faces.push_back({C0, h, ring[i-0]});
+          if(f & cox_othercell) {
+            ss.faces.push_back({facectr, ring[i-2], ring[i-1], ring[i-0]});
+            }
+          else {
+            ss.faces.push_back({fc2, ring[i-1], ring[i-0]});
+            ss.faces.push_back({fc2, ring[i-2], ring[i-1]});
+            }
+          }
+        }
+      else {
+        ssh.emplace_back();
+        auto &ss = ssh.back();
+        for(int i=1; i<isize(ring); i++)
+          ss.faces.push_back({C0, ring[i-1], ring[i]});
+        if(f & cox_othercell) {
+          ring.pop_back();
+          ss.faces.push_back(ring);
+          }
+        else {
+          for(int i=1; i<isize(ring); i++)
+            ss.faces.push_back({fc2, ring[i-1], ring[i]});
+          }
+        }
+      }
+    }
+  
   EX void generate_special_subcubes(bool bch) {
     const int sub = subcube_count;
     if(1) {
@@ -510,6 +580,10 @@ EX namespace reg3 {
         generate_special_subcubes(true);
         break;
       
+      case eVariation::coxeter:
+        generate_coxeter(coxeter_param);
+        break;
+    
       case eVariation::pure: {
         cgi.subshapes.emplace_back();
         cgi.subshapes[0].faces = cgi.cellshape;
diff --git a/system.cpp b/system.cpp
index 7cb2a031..26a39fd8 100644
--- a/system.cpp
+++ b/system.cpp
@@ -1382,6 +1382,11 @@ EX void set_variation(eVariation target) {
       variation = target;
       return;
       }
+    if(is_reg3_variation(target)) {
+      if(!reg3::in()) geometry = hyperbolic ? gSpace435 : gCell8;
+      variation = target;
+      return;
+      }
     if(bt::in() || sol || kite::in() || WDIM == 3) if(!prod) geometry = gNormal;
     auto& cd = ginf[gCrystal];
     if(target == eVariation::bitruncated && cryst && cd.sides == 8 && cd.vertex == 4) {