Added grid source precomputing and test functions

include/fonction.h

@@ -417,9 +417,11 @@ double einasto_gamma( double r, double rs, int n, double rhos, double kappa_crit
 struct point   einasto_gamma_eps(double r, double rs, double n, double theta, double kappa_crit, double rhos,double eps);
 float einasto_phi(double r, double rs, int n, double rhos, double kappa_crit);
 float einasto_alpha( double r, double rs, int n, double rhos, double kappa_crit);
-
 void read_table_einasto();
 
+void precompsource();
+void testpixel();
+
 //#ifdef __cplusplus
 //};
 //#endif

src/o_global.c

@@ -31,6 +31,7 @@ double chi_im, chip, chix, chiy, chia, chil, chis, chi_vel,chi_mass;    //I adde
 //double amplifi[NFMAX][NIMAX];
 double **imo, **wo, **soo, **ero;
 double ***cubeo,***ercubeo;
+double **fluxmap;
 double drima;   /* distance ratio between the image to inverse and the main lens*/
 double **sm_p;      /* optimization spline map : initial save */
 
@@ -279,6 +280,7 @@ void o_global_free()
 
         free_square_double(imo, imFrame.ny);
         free_square_double(wo, wFrame.ny);
+        free_square_double(fluxmap, imFrame.ny*imFrame.nx);
 
         // Reset source positions to absolute positions
         if( I.n_mult > 0 )

src/o_pixsource.c

@@ -13,6 +13,114 @@
 /*      place:      Lyon                */
 /****************************************************************/
 
-double ** precompsource(){
-extern struct g_frame F;
+void precompsource(){
+   const extern  struct  pot lens[];
+   extern struct g_mode M; // M.iclean
+   const extern struct g_observ  O;
+   extern struct g_pixel imFrame,ps,PSF;
+   extern double **fluxmap;
+
+   int ni=imFrame.nx*imFrame.ny;
+   int ns=ps.nx*ps.ny;
+
+   fluxmap=(double **)alloc_square_double(ni,ns);
+
+   double dlsds = dratio(lens[0].z, M.zclean);
+   double sizeratio=imFrame.pixelx/ps.pixelx; // pixel size ratio between image and source plane
+   double xi,yi,xs,ys,xs2,ys2;
+
+   double **tmpimage=(double **)alloc_square_double(imFrame.nx,imFrame.ny);
+   int is,js,i,j;
+
+   for(is=0;is<ps.nx;is++)
+   {
+      for(js=0;js<ps.ny;js++)
+      {
+        xs = is * (ps.xmax - ps.xmin) / (ps.nx - 1) + ps.xmin;
+        ys = js * (ps.ymax - ps.ymin) / (ps.ny - 1) + ps.ymin;
+        for(i=0;i<imFrame.nx;i++)
+        {
+          for(j=0;j<imFrame.ny;j++)
+          {
+            xi = i * (imFrame.xmax - imFrame.xmin) / (imFrame.nx - 1) + imFrame.xmin;
+            yi = j * (imFrame.ymax - imFrame.ymin) / (imFrame.ny - 1) + imFrame.ymin;
+            //e_dpl(xi,yi, dlsds, xs2,ys2);
+            if((xs2-xs)*(xs2-xs)+(ys2-ys)*(ys2-ys)<0.25) tmpimage[i][j]=sizeratio;
+            else tmpimage[i][j]=0.0;
+          }
+        }
+        //convolve tmpimage with PSF
+        if (O.setseeing)
+               d_seeing_omp(tmpimage, imFrame.nx, imFrame.ny, imFrame.pixelx);
+
+        //store tmpimage in flux matrix
+        for(i=0;i<imFrame.nx;i++)
+        {
+          for(j=0;j<imFrame.ny;j++)
+          {
+              fluxmap[i*imFrame.ny+j][is*ps.ny+js]=tmpimage[i][j];
+          }   
+        }
+     }
+   }
+   free_square_double(tmpimage,imFrame.ny);
+}
+
+void    testpixel(double *sflux)
+{
+    const extern struct g_observ  O;
+    extern struct g_mode    M;
+    extern struct g_pixel   ps,imFrame;
+    extern double **fluxmap;
+
+    double xmin, ymin, xmax, ymax, dx, dy;
+    double f, scale;
+    int nx, ny, ny_sav;
+    double **zz;
+    int i,j,k;
+
+    xmin = imFrame.xmin;
+    xmax = imFrame.xmax;
+    ymin = imFrame.ymin;
+    ymax = imFrame.ymax;
+
+    dx = xmax - xmin;
+    dy = ymax - ymin;
+    nx = imFrame.nx;    // default: assign np to nx and adjust ny
+    scale = dx / (nx-1);
+    ny = dy / scale + 1;  // warning: trunc(ny) behavior 
+
+    dx = (nx - 1) * scale;
+    dy = (ny - 1) * scale;
+    xmax = xmin + dx;
+    ymax = ymin + dy;
+
+    NPRINTF(stderr, "\timage (%d,%d) s=%.3lf [%.3lf,%.3lf] [%.3lf,%.3lf]\n",
+            nx, ny, scale, xmin, xmax, ymin, ymax);
+
+    zz = (double **) alloc_square_double(ny, nx);
+
+    int ns=ps.nx*ps.ny;
+    for(i=0;i<nx;i++)
+    {
+       for(j=0;j<ny;j++)
+       {
+          zz[i][j]=0.0;
+          for(k=0;k<ns;k++)
+          {
+             zz[i][j]+=fluxmap[k][i*ny+j]*sflux[k];
+          }
+       }
+    }
+
+    if (O.setseeing)
+        d_seeing(zz, nx, ny, scale);
+
+    if ( M.iref > 0 )
+        wrf_fits_abs("testpixel.fits", zz, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
+    else
+        wrf_fits("testpixel.fits", zz, nx, ny, xmin, xmax, ymin, ymax);
+
+    free_square_double(zz, ny_sav);
+
 }