#-------------------------------------------------------------------------------
# INSTRUCTIONS
# To patch from mangle1.4.1 to mangle1.4.3:
#
# cd <mangle_directory>
# patch -p0 < path_to/mangle1.4.1_to_1.4.3.patch
#
#-------------------------------------------------------------------------------
--- CHANGES.orig	2004-03-31 13:17:43.000000000 -0700
+++ CHANGES	2006-09-14 10:01:31.000000000 -0600
@@ -1,6 +1,104 @@
 --------------------------------------------------------------------------------
+mangle1.4.3
+14 Sep 2006
+
+Thanks to Colin Hill, working with Max Tegmark and Molly Swanson at MIT,
+for finding a polygon that defeated mangle.
+
+1. The problem was numerical round-off,
+and the immediate fix was to replace line 71 in split_poly.c
+	if (area != area_tot) {			/* boundary intersects poly1 */
+with
+	if (area > area_tot) {			/* boundary intersects poly1 */
+It should never happen that area < area_tot, since removing a boundary
+of a polygon should never decrease its area, but thanks to numerical round-off
+the area did decrease.  Once that happened, mangle span its wheels.
+
+2. The killer polygon that was responsible for the numerical failure of item 1
+is discussed below.  To improve mangle's defence against such polygons,
+I changed the strategy for modifying the tolerance angle to multiple
+intersections.  Originally, if mangle detected an inconsistency in the
+topology of the distribution of vertices around multiple intersections,
+then mangle would double the tolerance angle and try again.
+The new strategy is to try tightening as well as loosening the tolerance angle.
+Mangle now see-saws between and tighter and looser tolerance angles
+successive factors of 2 away from the original input tolerance angle.
+
+3. There was also a bug on line 182 of balkanize.c.
+    for (i = npoly; i < npolys; i++) {
+should have been
+    for (i = 0; i < npolys; i++) {
+This bug became evident as a result of a compiler change.
+When tripped, the bug produces a segmentation violation, so if you never
+experienced this bug before, it should never have caused an error.
+
+This is Colin Hill's killer polygon:
+polygon 6 ( 5 caps, 1.0 weight, 0.000003265281337 str):
+ -0.2972995896945569 -0.5791945044073898 -0.7590432662449000 0.5786075961488696
+ -0.1092635942528740  0.6422417176858356 -0.7586745303718128 0.5279437897883537
+ -0.3927631897959760 -0.5199166454922474 -0.7585669110114012 -0.4794583692304901
+ -0.0457867636138749  0.6498588533234551 -0.7586745303718132 -0.5279437897887224
+ -0.3467216633109536 -0.5515406106135746 -0.7586745303718130 0.5279437897887302
+It is a long (4.3 degrees) thin (36 arcsec) rectangle, split by a diagonal
+which is almost, but not quite, tangent at each end to the long direction.
+Each near tangent end is both almost multiply intersecting, and almost kissing.
+
+If you balkanize the killer polygon with the default tolerance angle
+for multiple intersections of 10^-5 arsec, then mangle1.4.1 will discard the
+polygon as having zero area:
+
+% balkanize killer.poly -
+---------------- balkanize ----------------
+snap angles: axis 2s latitude 2s edge 2s
+multiple intersections closer than 1e-05s will be treated as coincident
+...
+1 polygons read from j6
+warning from balkanize: following polygons have zero area & are being discarded:
+ 0
+...
+ 0 polygons written to output
+
+Well, the polygon is thin, 36 arcsec wide, so maybe you don't mind losing
+the polygon.  But with snap tolerance angles of 2 arcsec, you'd think mangle
+should keep it.  It turns out that the multiple intersections at each end
+of the polygon have separations that vary from 10^-8 arcsec in the thin
+direction to 80 arcsec in the long direction, and this distribution of
+separations conspires so that mangle finds no consistent topology
+(satisfying the 64-bit check number) for tolerance angles to multiple
+intersections anywhere between 10^-8 arcsec and 40 arcsec.  But 40 arcsec
+exceeds the 36 arcsec width of the polygon, so mangle concludes that the
+area is zero.
+
+If the tolerance angle for multiple intersections is tightened to
+5 x 10^-9 arcsec:
+
+% balkanize -m5e-9s killer.poly -
+
+then mangle does find the correct solution.  That is, since the polygon
+is already a valid polygon, balkanize spits back the killer polygon
+unchanged.
+
+Thus the solution to dealing with this polygon is not to loosen the
+tolerance, but to tighten it.  The new mangle strategy does this.
+
+Max Tegmark suggests switching the entire mangle code from double (64-bit)
+to quadruple (128-bit) precision.  This may be a sensible thing to do
+at some point.  The need for high precision was always evident.
+Inevitably, people are messing with masks that have arcsecond precision.
+A 1 x 1 arcsec rectangle has an area of 2 x 10^-11 str, which is 2 x 10^-12
+of the whole sky.  This is a factor of 10^4 larger than the effective
+precision 2 x 10^-16 of 64-bit floats.  The factor of 10^4 is ok,
+but with ever more complicated masks appearing, it may no longer be
+providing a sufficient safety margin.
+
+--------------------------------------------------------------------------------
+5 Sep 2005
+
+Added '\r' to characters considered blank, so polygon files are dos-compatible.
+
+--------------------------------------------------------------------------------
 mangle1.4.1
-30 Sep 2004
+30 Mar 2004
 
 Many thanks to Michael Blanton, who is working on the SDSS mask,
 for bringing these two problems to light.
--- src/balkanize.c.orig	2004-03-30 12:34:19.000000000 -0700
+++ src/balkanize.c	2006-09-10 09:57:19.000000000 -0600
@@ -8,6 +8,11 @@
 #include "manglefn.h"
 #include "defaults.h"
 
+/* define CARRY_ON_REGARDLESS if you want balkanize() to continue even when the number of polygons hits NPOLYSMAX;
+   if CARRY_ON_REGARDLESS is defined, then balkanize() will create a possibly incomplete polygon file of polygons */
+#undef	CARRY_ON_REGARDLESS
+//#define	CARRY_ON_REGARDLESS
+
 /* getopt options */
 const char *optstr = "dqa:b:t:y:m:s:e:v:p:i:o:";
 
@@ -174,7 +179,7 @@
     msg("balkanizing %d polygons ...\n", np);
 
     /* nullify all output polygons */
-    for (i = npoly; i < npolys; i++) {
+    for (i = 0; i < npolys; i++) {
 	polys[i] = 0x0;
     }
 
@@ -244,7 +249,12 @@
 		if (n > npolys) {
 		    fprintf(stderr, "balkanize: total number of polygons exceeded maximum %d\n", npoly + npolys);
 		    fprintf(stderr, "if you need more space, enlarge NPOLYSMAX in defines.h, and recompile\n");
+		    n = npolys;
+#ifdef	CARRY_ON_REGARDLESS
+		    break;
+#else
 		    return(-1);
+#endif
 		}
 	    }
 	    /* copy down non-null polygons */
@@ -262,6 +272,8 @@
 	    n = m + dm;
 	    if (dm == 0) break;
 	}
+	/* too many polygons */
+	if (n >= npolys) break;
 	ip++;
     }
     msg("added %d polygons to make %d\n", dnp, np);
@@ -296,7 +308,12 @@
 	if (n > npolys) {
 	    fprintf(stderr, "balkanize: total number of polygons exceeded maximum %d\n", npoly + npolys);
 	    fprintf(stderr, "if you need more space, enlarge NPOLYSMAX in defines.h, and recompile\n");
+	    n = npolys;
+#ifdef	CARRY_ON_REGARDLESS
+	    break;
+#else
 	    return(-1);
+#endif
 	}
 	ip++;
     }
--- src/garea.s.f.orig	2003-09-08 14:23:16.000000000 -0600
+++ src/garea.s.f	2006-09-14 09:28:59.000000000 -0600
@@ -26,7 +26,7 @@
 C     logical warn
       logical whole
       real*8 bik,cmi,cmik,cmk,d,darea,dph,ikchk,ikran,
-     *  ph,phm,php,psi,si,xi(3),yi(3)
+     *  ph,phm,php,psi,si,tolin,xi(3),yi(3)
 c *
 c * Area of surface of sphere of unit radius bounded by
 c *    1 - r.rp(i) < cm(i)  (if cm(i).ge.0)
@@ -57,8 +57,10 @@
 c        possible multiple intersection when dph < dphmin
       data dphmin /1.d-8/
 c
+c        input tolerance to multiple intersections
+      tolin=tol
 C     print *,'--------------------'
-c        come here with enlarged tolerance to multiple intersections
+c        come here with modified tolerance to multiple intersections
   100 continue
 C     write (*,'(a3,a20,4a24)')
 C    *  ' ','x','y','z',
@@ -129,6 +131,8 @@
         elseif (ni.gt.0) then
 c        find ordering of intersection angles around i circle
           call findbot(phi,2*np,iord,ni)
+C         write (*,'("phi")')
+C         write (*,'(i4,2g24.16)') (j,(phi(k,j),k=1,2),j=1,np)
 c........contribution to area from each segment of i circle
           jpl=0
 c        come here to do another segment
@@ -150,7 +154,7 @@
 C             warn=.true.
 C             print *,
 C    *          '*** warning from garea: near multiple intersection at',
-C    *          i,': edge',km,kp,' dph=',dph
+C    *          i,': edge',km,kp,' dph =',dph
             endif
 c........segment satisfies conditions
             nbd=nbd+1
@@ -229,26 +233,18 @@
 c       write (*,'(i3,4g24.16)')
 c    *    (j,(rp(i,j),i=1,3),
 c    *    cm(j),j=1,np)
-c        retry with enlarged tolerance
-        if (tol.le.0.d0) then
-          tol=1.d-15
-        else
-          tol=tol*2.d0
-        endif
+c        retry with modified tolerance
+        call gtol(tol,tolin)
         goto 100
       elseif (tol.gt.0.d0) then
 C       print *,'... from garea: success at tol =',tol
       endif
 c--------add/subtract 2*pi's to area
       retry=garpi(area,iarea,rp,cm,np,whole,nbd0m,nbd0p,nbd,nmult)
-c        retry with enlarged tolerance
+c        retry with modified tolerance
       if (retry.eq.1) then
 C       warn=.true.
-        if (tol.le.0.d0) then
-          tol=1.d-15
-        else
-          tol=tol*2.d0
-        endif
+        call gtol(tol,tolin)
         goto 100
       endif
 c--------done
--- src/gspher.s.f.orig	2003-09-08 14:24:48.000000000 -0600
+++ src/gspher.s.f	2006-09-14 09:12:31.000000000 -0600
@@ -28,7 +28,7 @@
       logical whole
       real*8 bik,ci,cmi,cmik,cmk,cti,ctk,ctpsi,
      *  d,darea,dbound(2),dph,dvert(2),ikchk,ikran,
-     *  ph,phi,phii,phm,php,psi,psip,ri,rii,si,sk,sqrt4pi,t,
+     *  ph,phi,phii,phm,php,psi,psip,ri,rii,si,sk,sqrt4pi,t,tolin,
      *  xi(3),yi(3)
 c *
 c * Spherical transform of region W of sphere of unit radius bounded by
@@ -173,8 +173,10 @@
 c        possible multiple intersection when dph < dphmin
       data dphmin /1.d-8/
 c
+c        input tolerance to multiple intersections
+      tolin=tol
 C     print *,'--------------------'
-c        come here with enlarged tolerance
+c        come here with modified tolerance
   100 continue
 c        initialise error flag to no error
       ldegen=.false.
@@ -534,11 +536,7 @@
 c       print *,'*** from gspher: at tol =',tol,
 c    *    ', ikchk=',ikchk,' should be 0'
 c        retry with enlarged tolerance
-        if (tol.le.0.d0) then
-          tol=1.d-15
-        else
-          tol=tol*2.d0
-        endif
+        call gtol(tol,tolin)
         goto 100
       endif
 c--------add/subtract 2*pi's to area
@@ -550,14 +548,10 @@
         i=nint((w(1,1)*sqrt4pi+area)/TWOPI)
       endif
       w(1,1)=w(1,1)-i*TWOPI/sqrt4pi
-c        retry with enlarged tolerance
+c        retry with modified tolerance
       if (retry.eq.1) then
 C       warn=.true.
-        if (tol.le.0.d0) then
-          tol=1.d-15
-        else
-          tol=tol*2.d0
-        endif
+        call gtol(tol,tolin)
         goto 100
       endif
 c--------done
--- src/gsubs.s.f.orig	2003-09-08 14:25:51.000000000 -0600
+++ src/gsubs.s.f	2006-09-14 09:26:42.000000000 -0600
@@ -1261,3 +1261,29 @@
       return
       end
 c
+c-----------------------------------------------------------------------
+      subroutine gtol(tol,tolin)
+      real*8 tol,tolin
+c *
+c * Modify tolerance tol to multiple intersections.
+c * The tolerance tol is changed by successive factors of 2
+c * from the original input tolerance tolin,
+c * see-sawing between smaller and larger values.
+c *
+      if (tolin.le.0.d0) then
+        if (tol.le.0.d0) then
+          tol=1.d-15
+        else
+          tol=tol*2.d0
+        endif
+c        see-saw tolerance between smaller and larger values
+      else
+        if (tol.ge.tolin) then
+          tol=tolin*tolin/tol/2.d0
+        else
+          tol=tolin*tolin/tol
+        endif
+      endif
+      return
+      end
+c
--- src/gvert.s.f.orig	2003-09-08 14:26:01.000000000 -0600
+++ src/gvert.s.f	2006-09-14 09:12:35.000000000 -0600
@@ -26,7 +26,7 @@
       integer i,ii,ik,iseg,iv,ive,j,jm,jml,jmu,jp,jpl,jpu,km,kp,
      *  mve,ni,scmi
 C     logical warn
-      real*8 amve,cmi,dph,ikchk,ph,phm,php,si,xi(3),xv,yi(3),yv,zv
+      real*8 amve,cmi,dph,ikchk,ph,phm,php,si,tolin,xi(3),xv,yi(3),yv,zv
 c *
 c * Vertices, plus points on edges, of the polygon defined by
 c *    1 - r.rp(i) < cm(i)  (if cm(i).ge.0)
@@ -120,8 +120,10 @@
 c     data dphmin /1.d-8/
       data big /1.d6/
 c
+c        input tolerance to multiple intersections
+      tolin=tol
 C     print *,'--------------------'
-c        come here with enlarged tolerance
+c        come here with modified tolerance
   100 continue
 c        initialise error flag to no error
       ldegen=.false.
@@ -308,12 +310,8 @@
 C       warn=.true.
 C       print *,'*** from gvert: at tol =',tol,
 C    *    ', ikchk=',ikchk,' should be 0'
-c        retry with enlarged tolerance
-        if (tol.le.0.d0) then
-          tol=1.d-15
-        else
-          tol=tol*2.d0
-        endif
+c        retry with modified tolerance
+        call gtol(tol,tolin)
         goto 100
       elseif (tol.gt.0.d0) then
 C       print *,'... from gvert: success at tol =',tol
--- src/gvlim.s.f.orig	2003-09-08 14:26:15.000000000 -0600
+++ src/gvlim.s.f	2006-09-14 09:12:40.000000000 -0600
@@ -28,7 +28,7 @@
       integer i,ii,ik,iphi,iseg,iv,j,jm,jml,jmu,jp,jpl,jpu,km,kp,ni,scmi
 C     logical warn
       real*8 cmi,cmik,dph,ikchk,ph,phin,phif,phm,phmax,phmin,php,
-     *  si,sik,xi(3),xv,yi(3),yv,zv
+     *  si,sik,tolin,xi(3),xv,yi(3),yv,zv
 c *
 c * Lifted mostly from gvert and gcmlim.
 c *
@@ -97,7 +97,9 @@
 c     data dphmin /1.d-8/
       data big /1.d6/
 c
-c        come here with enlarged tolerance
+c        input tolerance to multiple intersections
+      tolin=tol
+c        come here with modified tolerance
   100 continue
 c        initialise error flag to no error
       ldegen=.false.
@@ -321,12 +323,7 @@
       if (ikchk.ne.0.d0) then
 c       print *,'*** from gvlim: ikchk=',ikchk,' should be 0'
 C       warn=.true.
-c        retry with enlarged tolerance
-        if (tol.le.0.d0) then
-          tol=1.d-15
-        else
-          tol=tol*2.d0
-        endif
+        call gtol(tol,tolin)
         goto 100
       endif
 c--------order vertices right-handedly about polygon
--- src/rdangle.c.orig	2003-06-03 10:29:33.000000000 -0600
+++ src/rdangle.c	2005-09-05 09:13:01.000000000 -0600
@@ -20,7 +20,7 @@
 */
 int rdangle(char *word, char **next, char unit, double *angle)
 {
-    const char *blank = " \t\n";
+    const char *blank = " \t\n\r";
     const char *number = "+-.0123456789eE";
 
     int i, ird, sgn;
--- src/rdmask.c.orig	2004-03-31 13:17:43.000000000 -0700
+++ src/rdmask.c	2005-09-05 09:12:51.000000000 -0600
@@ -114,7 +114,7 @@
 */
 int rdmask(char *name, format *fmt, int npolys, polygon *polys[/*npolys*/])
 {
-    char input[] = "input";
+    char *input = "input";
     char *line_rest, *word;
     int ird;
     int npoly = 0;
@@ -205,7 +205,7 @@
 */
 char *get_keyword(char *str, char **str_rest, format *fmt)
 {
-    const char *blank = " \t\n";
+    const char *blank = " \t\n\r";
     char c;
     int id, ird, n;
 
@@ -267,7 +267,7 @@
     const char *end_fmt = "%d";
     const char *unit_fmt = " %c";
     char *word;
-    char blank[] = " \t\n";
+    char *blank = " \t\n\r";
     int ird, iscan, nholes;
     size_t word_len;
 
@@ -483,8 +483,8 @@
 int get_n(format *fmt, char *lin)
 {
 #define WARNMAX		3
-    char blank[] = " \t\n";
-    char plural[] = "s";
+    char *blank = " \t\n\r";
+    char *plural = "s";
     static int warn = 0;
     char *ch;
     int n;
@@ -778,7 +778,7 @@
     static vertices *vert = 0x0;
     static int *ev = 0x0;
 
-    const char *blank = " \t\n";
+    const char *blank = " \t\n\r";
     char unit;
     char *next, *word;
     int anti, i, iang, iedge, iev, ird, iv, ive, ivert, nedge, nv, nve, nvert, reverse;
--- src/split_poly.c.orig	2003-06-15 15:35:31.000000000 -0600
+++ src/split_poly.c	2006-09-12 08:52:50.000000000 -0600
@@ -68,8 +68,7 @@
 	poly->cm[np1 + ip] = 2.;		/* suppress boundary to be tested */
 	ier = garea(poly, &tol, verb, &area);	/* area of intersection sans boundary */
 	poly->cm[np1 + ip] = cm;		/* restore tested boundary */
-	if (area != area_tot) {			/* boundary intersects poly1 */
-
+	if (area > area_tot) {			/* boundary intersects poly1 */
 	    /* poly2 splits poly1, but do not actually split */
 	    if (!poly3) return(2);
 
@@ -112,6 +111,10 @@
 
 	    /* poly1 successfully split into poly1 and poly3 */
 	    return(2);
+	} else if (area < area_tot) {
+	    /* area should be >= area_tot because suppressing a boundary of poly should always increase its area;
+	       but this can happen because of roundoff */
+	    //fprintf(stderr, "split_poly: area %.16g of polygon with boundary %d suppressed should be >= area %.16g of polygon\n", area, ip, area_tot);
 	}
     }
 
