printf("WARNING: You are using SIE potentials. An image to close to one of the potential centers has been classified as numerical error and removed \n");
}
}
if(!skip_image)
{
//checking whether a closest image has already been found
nimagesfound[image_id][im_index]++;// increasing the total number of images found (If we find more than 1 theoretical image linked to 1 real image, these theoretical
// images are included in this number)
}
}
thread_found_image=0;// for next iteration
}
//@@p2_time += myseconds();
}
}
//inner_time += myseconds();
}
//////////////////////////////////////computing the local chi square//////////////////////////////////////
printf("WARNING: You are using SIE potentials. An image to close to one of the potential centers has been classified as numerical error and removed \n");
tim[image_id][im_index]=temp;// if we found an observed and not already linked image, we allocate the theoretical image temp
}
}
#pragma omp atomic
nimagesfound[image_id][im_index]++;// increasing the total number of images found (If we find more than 1 theoretical image linked to 1 real image, these theoretical
// images are included in this number)
}
}
thread_found_image=0;// for next iteration
}
//p2_time += myseconds();
}
}
//#pragma omp barrier
loop_time+=myseconds();
}
//____________________________ end of image loop
//
//____________________________ computing the local chi square
//
doublechiimage;
//
chi2_time-=myseconds();
int_nimages=nimages_strongLensing[source_id];
//
for(intiter=0;iter<_nimages*_nimages;iter++)
{
inti=iter/nimages_strongLensing[source_id];
intj=iter%nimages_strongLensing[source_id];
//printf("nimagesfound %d %d = %d\n", i, j, nimagesfound[i][j]);
if(i!=j)
{
// In the current method, we get the source in the source plane by ray tracing image in nimagesfound[i][i]. If we ray trace back,
// we arrive again at the same position and thus the chi2 from it is 0. Thus we do not calculate the chi2 (-> if i!=j)
if(nimagesfound[i][j]>0)
{
chiimage=pow(images[index+j].center.x-tim[i][j].x,2)+pow(images[index+j].center.y-tim[i][j].y,2);// compute the chi2
//printf("%d %d = %.15f\n", i, j, chiimage);
*chi+=chiimage;
}
elseif(nimagesfound[i][j]==0)
{
// If we do not find a correpsonding image, we add a big value to the chi2 to disfavor the model
*chi+=100.*nimages_strongLensing[source_id];
}
}
}
//
//____________________________ end of computing the local chi square
//
//printf("%d: chi = %.15f\n", source_id, *chi);
chi2_time+=myseconds();
/*
for (int i=0; i < nimages_strongLensing[source_id]; ++i){
for (int j=0; j < nimages_strongLensing[source_id]; ++j){
printf(" %d",nimagesfound[i][j]);
}
printf("\n");
}*/
//Incrementing Index: Images already treated by previous source_id
index+=nimages_strongLensing[source_id];
}
time+=myseconds();
intnthreads=1;
#pragma omp parallel
nthreads=omp_get_num_threads();
printf(" chi2 time = %f s. using %d threads\n",time,nthreads);
printf(" - trans time = %f s.\n",trans_time);
printf(" - loop time = %f s.\n",loop_time/nthreads);
printf(" - trans2 time = %f s.\n",trans2_time/nthreads);
printf(" - p1 time = %f s.\n",p1_time/nthreads);
printf(" - p2 time = %f s.\n",p2_time/nthreads);
printf(" - chi2 update time = %f s.\n",chi2_time);
printf(" images found: %d out of %d\n",images_found,images_total);
free(grid_gradient_x);
free(grid_gradient_y);
}
/** @brief Tranform a point from image to source plane. Result stored in sourcepoint argument
*
* Tranform a point from image to source plane using lensequation
*
* @param image_point image position
* @param dlsds dls/ds
* @param nhalos number of halos
* @param potential_param gravitational potential information
* @param source_point address where source information will be stored