diff --git a/tests/test_ggg.py b/tests/test_ggg.py
index 176552fc..f382f0da 100644
--- a/tests/test_ggg.py
+++ b/tests/test_ggg.py
@@ -4766,13 +4766,11 @@ def test_map3_logmultipole():
 
     # Skip the lowest r values, to make sure the integral has s and t values sufficiently
     # smaller than R for the integral to come out right.
-    # Interestingly, with the multpole method, we don't need to scale up the true_map
-    # by the (1+5r0/L) factor to get the middle values right like we did in test_map3_logsas.
-    # But we need to cut out more low R values to get to the ones where rtol=0.1 works.
-    # And the B-mode values are noisier, so I had to increase the tolerance for those tests.
-    r = ggg.rnom1d[14:]
+    r = ggg.rnom1d[7:]
     print('r = ',r)
     true_map3 = 2816./243. *np.pi * gamma0**3 * r0**12 * r**6 / (L**2 * (r**2+r0**2)**8)
+    true_map3 *= (1+5*r0/L)  # cf. comment in test_map3_logsas
+    print('scale true_map3 by ',1+5*r0/L)
     print('true map3 = ',true_map3)
 
     ggg_map3 = ggg.calculateMap3(R=r)
@@ -4785,9 +4783,7 @@ def test_map3_logmultipole():
     print('max diff = ',max(abs(map3 - true_map3)))
     np.testing.assert_allclose(map3, true_map3, rtol=0.1, atol=2.e-9)
     for mx in (mapmapmx, mapmxmap, mxmapmap, mxmxmap, mxmapmx, mapmxmx, mx3):
-        #print('mx = ',mx)
-        #print('max = ',max(abs(mx)))
-        np.testing.assert_allclose(mx, 0., atol=4.e-9)
+        np.testing.assert_allclose(mx, 0., atol=2.e-9)
 
     # Target the range where we expect good results.
     # This is the same as we had for the regular LogSAS tests.  Both methods work well here.
diff --git a/treecorr/gggcorrelation.py b/treecorr/gggcorrelation.py
index c64897d5..d6eec079 100644
--- a/treecorr/gggcorrelation.py
+++ b/treecorr/gggcorrelation.py
@@ -816,17 +816,19 @@ def toSAS(self, **kwargs):
         # p3 is on the x axis:
         # s = p3 - p1
         # t = p2 - p1
+        # u = angle bisector of s, t
         # q1 = (s+t)/3.  (this is the centroid)
         # q2 = q1-t
         # q3 = q1-s
         s = sas.meand2
         t = sas.meand3 * np.exp(1j * sas.meanphi * sas._phi_units)
+        u = (1 + t/np.abs(t))/2
         q1 = (s+t)/3.
         q2 = q1-t
         q3 = q1-s
 
         # Currently the projection is as follows:
-        # g1 is projected along q1
+        # g1 is projected along u
         # g2 is projected along t
         # g3 is projected along s
         #
@@ -836,22 +838,25 @@ def toSAS(self, **kwargs):
         # g3 projected along q3
         #
         # The phases to multiply by are exp(2iphi_current) * exp(-2iphi_target). I.e.
-        # g1phase = 1
+        # g1phase = (u conj(q1))**2 / |u conj(q1)|**2
         # g2phase = (t conj(q2))**2 / |t conj(q2)|**2
         # g3phase = (s conj(q3))**2 / |s conj(q3)|**2
+        g1phase = (u * np.conj(q1))**2
         g2phase = (t * np.conj(q2))**2
         g3phase = (s * np.conj(q3))**2
+        g1phase /= np.abs(g1phase)
         g2phase /= np.abs(g2phase)
         g3phase /= np.abs(g3phase)
 
         # Now just multiply each gam by the appropriate combination of phases.
-        # Note: gam0,1 have the same phase correction and gam2,3 are just conjuagates.
-        gam0phase = g2phase * g3phase
-        gam2phase = np.conj(g2phase) * g3phase
+        gam0phase = g1phase * g2phase * g3phase
+        gam1phase = np.conj(g1phase) * g2phase * g3phase
+        gam2phase = g1phase * np.conj(g2phase) * g3phase
+        gam3phase = g1phase * g2phase * np.conj(g3phase)
         gam0 *= gam0phase
-        gam1 *= gam0phase
+        gam1 *= gam1phase
         gam2 *= gam2phase
-        gam3 *= np.conj(gam2phase)
+        gam3 *= gam3phase
 
         sas.gam0r = np.real(gam0)
         sas.gam0i = np.imag(gam0)