1st working version cleanset 5

include/dimension.h

@@ -19,7 +19,7 @@
 #define NFMAX 	50000    /* maximum number of families */
 #define NPAMAX 	40     // number of free parameters (see #define in structure.h)
 #define NTMAX 	1024
-#define NPOINT 	100000  /* Number of contour points in cleanlens mode*/
+#define NPOINT 	5000  /* Number of contour points in cleanlens mode*/
 #define NPARMAX 50	
 #define NMCMAX  600	
 #define NSRCFIT 30    // Number of points for source plane fitting

src/do_itos.c

@@ -2,6 +2,9 @@
 #include <constant.h>
 #include <fonction.h>
 
+
+static int  readlist(struct point I[NPOINT], char file[NIMAX]);
+
 /*
 *       nom:        do_itos - lesn-tool
 *       auteur:     Jean-Paul Kneib
@@ -132,3 +135,122 @@ void    do_itos(double **im, struct pixlist *pl, int npl, double dlsds, double z
             erreur[is][js] -= source[is][js] * source[is][js];
 
 }
+
+void    newdo_itos(double **im, double dlsds, double zs,
+                double **source, double **erreur,
+                int **imult)
+{
+    extern  struct g_pixel  imFrame;
+    extern struct point gsource_global[NGGMAX][NGGMAX];  
+    register    int ii, jj, k, is, js;
+    const extern  struct  g_mode  M;
+    const extern  struct g_pixel  ps;
+
+
+    double ech=1.0;
+
+    struct point I[10][NPOINT];
+    int npcont[10];
+
+
+    for (k = 0; k < imFrame.ncont; k++)
+    {
+        npcont[k] = readlist(I[k], imFrame.contfile[k]);
+    }
+
+    if(M.flux!=0)  ech=1.0/((double)ps.ech*ps.ech);
+
+    /*
+     * Initialise full grid
+     */
+    e_unlensgrid(gsource_global, dlsds, zs);
+
+    /*
+     * For each pixel in the source plane
+     */
+    for (is = 0; is < ps.nx; is++)
+    {
+	    for(js = 0; js < ps.ny; js++)
+	    {
+		    source[js][is]=0.0;
+                    imult[js][is] = 0;
+		    struct point spos;
+                    struct point pImage[NIMAX]; // list of image positions for a given source position
+
+		    spos.x=ps.xmin+((double)is+1)*ps.pixelx;
+		    spos.y=ps.ymin+((double)js+1)*ps.pixely;
+
+		    /*
+		     * Compute all image positions from this source location
+		     */
+		    int ni = e_lens_P(spos,pImage,dlsds,zs);
+
+		    /*
+		     * Loop over all images
+		     */
+		    for (k = 0; k < ni; k++)
+                    {
+		         /*
+		         * If the image is within one contour then put flux of closest pixel.
+		         */
+			 int kcont;
+                         for (kcont = 0; kcont < imFrame.ncont; kcont++)
+                         {
+                              if(inconvexe(pImage[k],npcont[kcont],I[kcont])>0)
+		              {
+				   ii=((pImage[k].x-imFrame.xmin))/imFrame.pixelx+0.5;
+				   jj=((pImage[k].y-imFrame.ymin))/imFrame.pixely+0.5;
+				   //do bilinear interpolation instead
+                                   imult[js][is] += 1;
+		                   source[js][is]=(source[js][is]*(imult[js][is]-1)+ech*im[(int)jj][(int)ii])/((double)imult[js][is]);
+				   if((is==31)&&(js==35))
+			           {
+					   printf("%d %d %f\n",ii,jj,im[(int)jj][(int)ii]);
+					   printf("x=%f y=%f\n",pImage[k].x,pImage[k].y);
+					   printf("xmin=%f ymin=%f\n",imFrame.xmin,imFrame.ymin);
+					   printf("pixelx=%f pixely=%f\n",imFrame.pixelx,imFrame.pixely);
+			           }
+			      }
+			 }
+		    }
+	    }
+    }
+            
+}
+
+/****************************************************************
+ * In the cleanlens mode, read the contour file
+ * Return
+ *  - I the array of point filled with the contour
+ *  - k the number of elements in I
+ */
+
+static int  readlist(struct point I[NPOINT], char file[NIMAX])
+{
+    int n, k = 0;
+    int stop = 0;
+    FILE    *IN;
+
+    IN = fopen(file, "r");
+    if ( IN != NULL )
+        while (stop != -1)
+        {
+            stop = fscanf(IN, "%d%lf%lf\n", &n, &I[k].x, &I[k].y);
+
+            if (stop == 3)
+                k++;
+            else if (stop > 0)
+            {
+                fprintf(stderr, "ERROR: %s is not in the right format\n", file);
+                exit(-1);
+            }
+        }
+    else
+    {
+        fprintf( stderr, "ERROR: Opening the contour file %s\n", file);
+        exit(-1);
+    }
+    fclose(IN);
+    return(k);
+}
+

src/imtosou.c

@@ -16,8 +16,6 @@
 *
 */
 
-struct g_mode   M;
-
 void    imtosou(double zimage, char *sname )
 {
 //  register    int i,j,k;
@@ -91,7 +89,7 @@ void    imtosou(double zimage, char *sname )
             imFrame.nx, imFrame.ny, ps.nx, ps.ny);
 
     if(M.iclean==1) do_itos(im, pl, npl, dlsds, zimage, source, erreur, imult);
-    if(M.iclean==5) newdo_itos(im, pl, npl, dlsds, zimage, source, erreur, imult);
+    if(M.iclean==5) newdo_itos(im, dlsds, zimage, source, erreur, imult);
 
     /*
     *  save the results