Updated zonemult for MarchingSquares

src/criticnew.c

@@ -573,10 +573,10 @@ static void cl_connectdots(int verbose)
     extern  struct  biline  radial[], tangent[];
     extern  int nrline, ntline;
 //#ifdef DEBUG
-        fprintf( stderr, "I connect the dots here\n" );
+        fprintf( stderr, "Connecting the marchingSquares segments....\n" );
 //#endif
 
-    float currentdist,distmin;
+    float currentdist,distmin,disttest;
     struct biline zpos,currentpos;
     int pos,minpos,iline;
     
@@ -590,6 +590,7 @@ static void cl_connectdots(int verbose)
 	 {
 		 iline++;
 		 currentpos=tangent[i];
+                 tangent[i].i=iline;
 		 continue;
 	 }
 	 pos=i;
@@ -605,6 +606,15 @@ static void cl_connectdots(int verbose)
 	    }
 	    pos++;
          }
+
+	 if(distmin>pow(CL.limitHigh,2.0))
+         {
+              iline++;
+	      currentpos=tangent[i];
+              tangent[i].i=iline;
+	      continue;
+	 }
+
 	 //swap position i and minpos
 	 zpos=tangent[i];
 	 tangent[i]=tangent[minpos];
@@ -623,6 +633,7 @@ static void cl_connectdots(int verbose)
 	 {
 		 iline++;
 		 currentpos=radial[i];
+                 radial[i].i=iline;
 		 continue;
 	 }
 	 pos=i;
@@ -630,7 +641,7 @@ static void cl_connectdots(int verbose)
 	 minpos=i;
 	 while((radial[pos].i!=1)&&(pos<nrline))
          {
-	    currentdist=pow(currentpos.I.x-radial[pos].I.x,2)+pow(currentpos.I.y-radial[pos].I.y,2);
+	    currentdist=pow(currentpos.I.x-radial[pos].I.x,2.0)+pow(currentpos.I.y-radial[pos].I.y,2.0);
 	    if((distmin<0)||(currentdist<distmin))
 	    {  
 		    minpos=pos;
@@ -638,6 +649,17 @@ static void cl_connectdots(int verbose)
 	    }
 	    pos++;
          }
+
+	 //disttest=pow(currentpos.I.x-radial[startpos].I.x,2)+pow(currentpos.I.y-radial[startpos].I.y,2);
+	 //if((disttest<distmin)&&(disttest>0))
+	 if(distmin>pow(CL.limitHigh,2.0))
+         {
+              iline++;
+	      currentpos=radial[i];
+              radial[i].i=iline;
+	      continue;
+	 }
+	 
 	 //swap position i and minpos
 	 zpos=radial[i];
 	 radial[i]=radial[minpos];

src/zonemult.c

@@ -24,30 +24,22 @@ void    zonemult()
     extern int nrline, ntline;//,flagr,flagt;
 
     register int    i, j, k;
+    register int inttest;
     int iline, nline[NLMAX], ncount, flag;
     double  dlsds;
-    struct  point   ps, pi, line[NIMAX][NPOINT];
+    //struct  point   ps, pi, line[NLMAX][NPOINT];
+    struct  point   ps, pi, **line;
     //int   size[4];
     double  xmin, ymin, xmax, ymax;
     int **nimage;
 
-    /* Definition de la fenetre de calcul */
 
-    if (CL.dmax != 0.)
-    {
-        xmin = -CL.dmax;
-        xmax = CL.dmax;
-        ymin = -CL.dmax;
-        ymax = CL.dmax;
-    }
-    else
-    {
-        xmin = F.xmin;
-        xmax = F.xmax;
-        ymin = F.ymin;
-        ymax = F.ymax;
-    }
+    /* Definition de la fenetre de calcul */
 
+    xmin = F.xmin;
+    xmax = F.xmax;
+    ymin = F.ymin;
+    ymax = F.ymax;
 
     /*    verification du nombre de plan en z pour le calcul des zones */
 
@@ -57,17 +49,23 @@ void    zonemult()
         return;
     }
 
-    if ( strcmp(CL.algorithm, "SNAKE") )
-    {
-        NPRINTF(stderr, "WARNING: image zone not computed. Critical line algorithm must be SNAKE\n");
-        return;
-    }
+    //if ( strcmp(CL.algorithm, "SNAKE") )
+    //{
+    //    NPRINTF(stderr, "WARNING: image zone not computed. Critical line algorithm must be SNAKE\n");
+    //    return;
+    //}
 
     /* calcul des zones images pour zs=CL.cz[0] */
     NPRINTF(stderr, "COMP: multiple images area in the Image plane for sources at z=%.3lf\n", CL.cz[0]);
 
+
     dlsds = dratio(lens[0].z, CL.cz[0]);
 
+    line = (struct point **)malloc(NLMAX*sizeof(struct point *));
+    for (int i=0;i<NLMAX;i++)
+	    line[i]= (struct point *)malloc(NPOINT*sizeof(struct point));
+
+
     iline = 0;
     if (ntline > 2)
     {
@@ -87,6 +85,7 @@ void    zonemult()
             }
         }
 
+
         if (iline != 2)
             nline[iline++] = ncount;
     }
@@ -111,16 +110,23 @@ void    zonemult()
         nline[iline++] = ncount;
     }
 
+
 #ifdef DEBUG
     NPRINTF(stderr, "DBG: Before alloc in zone_mult\n");
     NPRINTF(stderr, "DBG: CL.npzone %d\n", CL.npzone);
 #endif
 
     nimage = (int **) alloc_square_int(CL.npzone, CL.npzone);
+    if(nimage==NULL)
+    {
+            fprintf(stderr, "FATAL ERROR: allocation of zonemult\n");
+	    exit(-1);
+    }
 #ifdef DEBUG
     NPRINTF(stderr, "DBG: After alloc in zone_mult\n");
 #endif
 
+
     for (j = 0; j < CL.npzone; j++)
     {
         pi.y = j * (ymax - ymin) / (CL.npzone - 1) + ymin;
@@ -133,11 +139,12 @@ void    zonemult()
             nimage[j][i] = 1;
 
             for (k = 0; k < iline; k++)
-                nimage[j][i] += 2 * inconvexe(ps, nline[k], line[k]);
+                  nimage[j][i] += 2 * inconvexe(ps, nline[k], line[k]);
 
         };
     };
 
+
     if (M.iref >0) wri_fits_abs(CL.zonefile, nimage, CL.npzone, CL.npzone, F.xmin, F.xmax, F.ymin, F.ymax,M.ref_ra,M.ref_dec);
     else wri_fits(CL.zonefile, nimage, CL.npzone, CL.npzone, F.xmin, F.xmax, F.ymin, F.ymax);