Commit d53bb21f authored by Johan Richard's avatar Johan Richard
Browse files

Corrected issue in pixel simulation for multiplane lensing

parent 29a77378
......@@ -141,6 +141,8 @@ void e_pixel(int np, char *iname, char *sname, struct galaxie *source)
* Called by o_chi() and e_pixel()
* zz is filled in a field defined by champ section with xmin, ymin, xmax, ymax
*
*
* Notes: not optimal because of multiplane (gradient is computed multiple times per pixel, one for each source).
*/
void o_pixel(double **zz, int nx, int ny, double scale, double xmin,
double xmax, double ymin, double ymax, struct galaxie *source, double **dpl_x, double **dpl_y)
......@@ -150,42 +152,28 @@ void o_pixel(double **zz, int nx, int ny, double scale, double xmin,
const extern struct g_source S;
extern struct pot lens[];
int j, k, l;
struct point pi, ps;
double zs_ref = 100.;
double dlsds_ref = dratio(lens[0].z, zs_ref);
if (dpl_x == NULL)
for (j = 0; j < ny; j++)
{
dpl_x = (double **) alloc_square_double(ny, nx);
dpl_y = (double **) alloc_square_double(ny, nx);
for (j = 0; j < ny; j++)
for (k = 0; k < nx; k++)
{
pi.y = ymin + j * scale;
for (k = 0; k < nx; k++)
{
pi.x = xmin + k * scale;
e_dpl(&pi, dlsds_ref, zs_ref, &ps);
dpl_x[j][k] = pi.x - ps.x;
dpl_y[j][k] = pi.y - ps.y;
}
zz[j][k] = 0.0;
}
}
//#pragma omp parallel for private(ps,pi,j,k,l)
for (j = 0; j < ny; j++)
for (l = 0 ; l < S.ns ; l++ )
{
pi.y = ymin + j * scale;
for (k = 0; k < nx; k++)
for (j = 0; j < ny; j++)
{
pi.x = xmin + k * scale;
zz[j][k] = 0.;
for (l = 0 ; l < S.ns ; l++ )
pi.y = ymin + j * scale;
for (k = 0; k < nx; k++)
{
ps.x = pi.x - dpl_x[j][k] * source[l].dr / dlsds_ref;
ps.y = pi.y - dpl_y[j][k] * source[l].dr / dlsds_ref;
pi.x = xmin + k * scale;
e_dpl(&pi,source[l].dr,source[l].z,&ps);
double f = d_profil(ps.x, ps.y, &source[l]);
/* if ((f>sqrt(O.SKY/O.gain)) || (f>source[l].I0/100.))
zz[j][k]+=d_integrer(source[l],pi.x,pi.y,scale/2.,f,0);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment