if ( lens[i].type == 5 and fabs(im_position.x - lens[i].position.x) <= dx/2. and fabs(im_position.y - lens[i].position.y) <= dx/2.){
skip_image = 1;
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==0){
//checking whether a closest image has already been found
// Loop over all images in the set that are not yet allocated to a theoretical image
// and allocate the closest one
- for(int i=0; i<nimages_strongLensing[bid] && nimagesfound[image_id][i]==0; i++) // we search for an observed image not already linked (nimagesfound=0)
+ for(int i=0; i<nimages_strongLensing[source_id] && nimagesfound[image_id][i]==0; i++) // we search for an observed image not already linked (nimagesfound=0)
{
if(im_dist[i]<im_dist[second_closest_id])
{
second_closest_id=i;
im_index=i; // im_index value changes only if we found a not linked yet image
tim[image_id][im_index]=temp; // if we found an observed and not already linked image, we allocate the theoretical image temp
}
}
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
}
}
}
}
//////////////////////////////////////computing the local chi square//////////////////////////////////////
double chiimage;
for( int iter = 0; iter < nimages_strongLensing[source_id]*nimages_strongLensing[source_id]; iter++){
- int i=tid/nimages_strongLensing[iter];
- int j=tid % nimages_strongLensing[iter];
+ int i=iter/nimages_strongLensing[source_id];
+ int j=iter % nimages_strongLensing[source_id];
+
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
*chi += chiimage;
- //printf("chi %f %f %f bid %d tid %d X %f Y %f\n", chi,temp_chi,chiimage,bid,tid ,tim[i][j].x,tim[i][j].y);
}
else if(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];
- //
- //printf("chi %f %f %f bid %d tid %d X %f Y %f\n", chi,temp_chi,chiimage,bid,tid ,tim[i][j].x,tim[i][j].y);
}
}
}
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];
}
}
/** @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
*
*
*/
void module_chiClassic_transformImageToSourcePlane(const runmode_param *runmode, const struct point *image_point, double dlsds, const struct Potential *lens, struct point *source_point)
{ // dlsds is the distance between lens and source divided by the distance observer-source
struct point Grad; // gradient
Grad = module_potentialDerivatives_totalGradient(runmode,image_point, lens);