diff --git a/graph.cpp b/graph.cpp
index ef56e9e3..85c4439c 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -5026,11 +5026,6 @@ void queuecircleat(cell *c, double rad, int col) {
 
 void drawMarkers() {
 
-  if(pmodel == mdTwoPoint) {
-    queuechr(xpush(+vid.twopoint_param) * C0, 8, 'X', 0xFFFF00);
-    queuechr(xpush(-vid.twopoint_param) * C0, 8, 'X', 0xFFFF00);
-    }
-  
   if(!(cmode & sm::NORMAL)) return;
   
   for(cell *c1: crush_now) 
@@ -5500,18 +5495,24 @@ void drawfullmap() {
     
   ptds.clear();
 
-  if(!stereo::active() && sphere && pmodel == mdTwoPoint) {
+  if(pmodel == mdTwoPoint) {
+    queuechr(xpush(+vid.twopoint_param) * C0, vid.xres / 100, 'X', 0xFF0000);
+    queuechr(xpush(-vid.twopoint_param) * C0, vid.xres / 100, 'X', 0xFF0000);
+    }
+  
+  if(!twopoint_do_flips && !stereo::active() && sphere && pmodel == mdTwoPoint) {
     queuereset(vid.usingGL ? mdDisk : mdUnchanged, PPR_CIRCLE);
 
-    for(int a=0; a<=180; a++) {
+    for(int b=-1; b<=1; b+=2)
+    for(int a=-90; a<=90; a++) {
       using namespace hyperpoint_vec;
-      ld x = cos(a * M_PI / 90);
+      ld x = sin(a * vid.twopoint_param * b / 90);
       ld y = 0;
       ld z = -sqrt(1 - x*x);
       hyperpoint h1;
       applymodel(hpxyz(x,y,z), h1);
-      if(h1[1] < 0) h1[1] = -h1[1];
-      if(a >= 90) h1[1] = -h1[1];
+      
+      h1[1] = abs(h1[1]) * b;
       curvepoint(h1 * vid.radius);
       }
 
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 35d95560..ced9108a 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -66,6 +66,15 @@ ld asin_auto(ld x) {
     }
   }
 
+ld asin_auto_clamp(ld x) {
+  switch(cgclass) {
+    case gcEuclid: return x;
+    case gcHyperbolic: return asinh(x);
+    case gcSphere: return x>1 ? M_PI/2 : x<-1 ? -M_PI/2 : isnan(x) ? 0 : asin(x);
+    default: return x;
+    }
+  }
+
 ld cos_auto(ld x) {
   switch(cgclass) {
     case gcEuclid: return 1;
diff --git a/hypgraph.cpp b/hypgraph.cpp
index 68e50a60..e1f9d973 100644
--- a/hypgraph.cpp
+++ b/hypgraph.cpp
@@ -123,6 +123,8 @@ bool hypot_zlev(bool zlev_used, ld& d, ld zlev, ld& df, ld& zf) {
   return hypot_zlev(zlev_used, d, zlev, df, zf, z);
   }
 
+int twopoint_sphere_flips;
+bool twopoint_do_flips;
 
 void applymodel(hyperpoint H, hyperpoint& ret) {
   
@@ -266,7 +268,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
     else {
       ld x, y, yf, zf;
       y = asin_auto(H[1]);
-      x = asin_auto(H[0] / cos_auto(y));
+      x = asin_auto_clamp(H[0] / cos_auto(y));
       if(sphere) {
         if(H[2] < 0 && x > 0) x = M_PI - x;
         else if(H[2] < 0 && x <= 0) x = -M_PI - x;
@@ -276,6 +278,19 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
         auto p = vid.twopoint_param;
         ld dleft = hypot_auto(x-p, y);
         ld dright = hypot_auto(x+p, y);
+        if(sphere) {
+          int tss = twopoint_sphere_flips;
+          if(tss&1) { tss--; 
+            dleft = 2*M_PI - 2*p - dleft;
+            dright = 2*M_PI - 2*p - dright;
+            swap(dleft, dright);
+            y = -y;
+            }
+          while(tss) { tss -= 2;
+            dleft = 2*M_PI - 4*p + dleft;
+            dright = 2*M_PI - 4*p + dright;
+            }
+          }
         x = (dright*dright-dleft*dleft) / 4 / p;
         y = (y>0?1:-1) * sqrt(dleft * dleft - (x-p)*(x-p) + 1e-9);
         }
diff --git a/polygons.cpp b/polygons.cpp
index b690fbd2..5b9e2fc8 100644
--- a/polygons.cpp
+++ b/polygons.cpp
@@ -469,6 +469,10 @@ unsigned char& part(int& col, int i) {
 
 bool in_twopoint = false;
 
+ld glhypot2(glvertex a, glvertex b) { 
+  return (a[0]-b[0]) * (a[0]-b[0]) + (a[1]-b[1]) * (a[1]-b[1]) + (a[2]-b[2]) * (a[2]-b[2]);
+  }
+
 void drawpolyline(polytodraw& p) {
   auto& pp = p.u.poly;
 
@@ -485,43 +489,79 @@ void drawpolyline(polytodraw& p) {
     }
 
   if(sphere && pmodel == mdTwoPoint && !in_twopoint) {
-    static vector<glvertex> v0, v1;
-    v0.clear(); v1.clear();
-    for(int i=0; i<pp.cnt; i++) {
-      glvertex h = glhr::pointtogl(pp.V * glhr::gltopoint((*pp.tab)[pp.offset+i]));
-      if(h[2] >= 0)
-        v0.push_back(h), v1.push_back(h);
-      else if(h[1] < 0)
-        v0.push_back(h);
-      else if(h[1] > 0)
-        v1.push_back(h);
+    #define MAX_PHASE 4
+    vector<glvertex> phases[MAX_PHASE];
+    extern int twopoint_sphere_flips;
+    extern bool twopoint_do_flips;
+    int pha;
+    if(twopoint_do_flips) {
+      for(int i=0; i<pp.cnt; i++) {
+        hyperpoint h1 = pp.V * glhr::gltopoint((*pp.tab)[pp.offset+i]);
+        for(int j=0; j<MAX_PHASE; j++) {
+          twopoint_sphere_flips = j;
+          hyperpoint h2; applymodel(h1, h2);
+          using namespace hyperpoint_vec;
+          glvertex h = glhr::pointtogl(h2 * vid.radius);
+          if(i == 0)
+            phases[j].push_back(h);
+          else {
+            int best = -1;
+            ld bhypot = 1e60;
+            for(int j0=0; j0<MAX_PHASE; j0++)
+              if(size(phases[j0]) == i) {
+                ld chypot = glhypot2(phases[j0].back(), h);
+                if(chypot < bhypot || best == -1) bhypot = chypot, best = j0;
+                }
+            phases[best].push_back(h);
+            }
+          }
+        }
+      twopoint_sphere_flips = 0;
+      pha = MAX_PHASE-1;
+      }
+    else {
+      pha = 1;
+      if(true) {
+        // a
+        // b
+        // lin(a,b) is of form (x, 0, z)
+        int cpha = 0;
+        for(int i=0; i<pp.cnt; i++) {
+          using namespace hyperpoint_vec;
+
+          hyperpoint h1 = pp.V * glhr::gltopoint((*pp.tab)[pp.offset+i]);
+          hyperpoint mh1; applymodel(h1, mh1);
+          phases[cpha].push_back(glhr::pointtogl(mh1 * vid.radius));
+
+          // check if the i-th edge intersects the boundary of the ellipse
+          // (which corresponds to the segment between the antipodes of foci)
+          // if yes, switch cpha to the opposite
+          hyperpoint h2 = pp.V * glhr::gltopoint((*pp.tab)[pp.offset+(i+1)%pp.cnt]);
+          if(h1[1] * h2[1] > 0) continue;
+          ld c1 = h1[1], c2 = -h2[1];
+          if(c1 < 0) c1 = -c1, c2 = -c2;
+          hyperpoint h = h1 * c1 + h2 * c2;
+          h /= hypot3(h);
+          if(h[2] < 0 && abs(h[0]) < sin(vid.twopoint_param)) cpha = 1-cpha, pha = 2;
+          }
+        if(cpha == 1) pha = 0;
+        }
       }
-    auto tV = pp.V;
-    auto toffset = pp.offset;
-    auto ttab = pp.tab;
-    auto tcnt = pp.cnt;
     vector<glvertex> tv;
     if(pp.tinf) {
       for(int i=0; i<pp.cnt; i++)
         tv.push_back(pp.tinf->tvertices[pp.offset+i]);
       swap(pp.tinf->tvertices, tv);
       }
-    pp.V = Id; pp.offset = 0;
-    in_twopoint = true;
-    if(size(v0) == pp.cnt) {
-      pp.tab = &v0;
-      drawpolyline(p);      
-      }
-    else if(size(v1) == pp.cnt) {
-      pp.tab = &v1;
+    dynamicval<eModel> d1(pmodel, mdUnchanged);
+    dynamicval<transmatrix> d2(pp.V, Id);
+    dynamicval<int> d3(pp.offset, 0);
+    dynamicval<decltype(pp.tab)> d4(pp.tab, pp.tab);
+    for(int j=0; j<pha; j++) {
+      dynamicval<int> d5(pp.cnt, size(phases[j]));
+      pp.tab = &phases[j];
       drawpolyline(p);
       }
-    else {
-      if(size(v0) >= 3) { pp.tab = &v0; pp.cnt = size(v0); drawpolyline(p); }
-      if(size(v1) >= 3) { pp.tab = &v1; pp.cnt = size(v1); drawpolyline(p); }
-      }
-    in_twopoint = false;
-    pp.V = tV; pp.offset = toffset; pp.tab = ttab; pp.cnt = tcnt;
     if(pp.tinf) swap(pp.tinf->tvertices, tv);
     return;
     }