#define MAX_DEVIATION_AT_FIRST_STEP 6.0/3600.0 // 5.0/3600.0 //1.8/3600.0
//#define MAX_DEVIATION_AT_SECOND_STEP 0.5*MAX_DEVIATION_AT_FIRST_STEP //5.0/6.0*MAX_DEVIATION_AT_FIRST_STEP // 3.0/3600.0 //0.9/3600.0
#define REFERENCE_LOCAL_SOLUTION_RADIUS_DEG 1.0

#define MIN_APASS_MAG 1.0
#define MAX_APASS_MAG 18.0
#define MIN_APASS_MAG_ERR 0.0
#define MAX_APASS_MAG_ERR 1.0

#define DEFAULT_APPROXIMATE_FIELD_OF_VIEW_ARCMIN 40.0

//#define VIZIER_SITE "vizier.u-strasbg.fr"
//#define VIZIER_SITE "vizier.cfa.harvard.edu"
#define VIZIER_SITE "`lib/choose_vizier_mirror.sh`"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <libgen.h> // for basename()
#include <sys/types.h> // for getpid()
#include <unistd.h> // also for getpid() and unlink() ... 
#include <math.h>

#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>

#include "fit_plane_lin.h"
#include "fitsio.h"
#include "fitsfile_read_check.h"
#include "vast_limits.h"
#include "ident.h"

#include "wpolyfit.h"

struct str_catalog_search_parameters{
 double search_radius_deg;
 double search_radius_second_step_deg;
 double brightest_mag;
 double faintest_mag;
};

struct detected_star{
 int n_current_frame; // star number in the input SExtractor catalog
 double x_pix; // star position
 double y_pix;
 double ra_deg_measured;
 double dec_deg_measured;
 double flux;
 double flux_err;
 double mag;
 double mag_err;
 int flag;
 /////////////////////////
 double ra_deg_measured_orig;
 double dec_deg_measured_orig;
 //double aperture;
 double distance_from_image_edge;
 /////////////////////////
 int matched_with_catalog;
 char ucac4id[32];
 double d_ra;
 double d_dec;
 double computed_d_ra;
 double computed_d_dec;
 double corrected_ra;
 double corrected_dec;
 double corrected_mag_ra;
 double corrected_mag_dec;
 double local_correction_ra;
 double local_correction_dec;
 double corrected_ra_local;
 double corrected_dec_local;
 double catalog_ra;
 double catalog_dec;
 double catalog_mag;
 double catalog_mag_err;
 int good_star;
 /////////////////////////
 double APASS_B;
 double APASS_B_err;
 double APASS_V;
 double APASS_V_err;
 double APASS_r;
 double APASS_r_err;
 double APASS_i;
 double APASS_i_err;
 /////////////////////////
 double Rc_computed_from_APASS_ri;
 double Rc_computed_from_APASS_ri_err;
 /////////////////////////
 double estimated_local_correction_accuracy;
};

/*
void remove_outliers_from_an_array(double * a, int * N_good){
 int i,j;
 double median;
 double sigma,MAD,M;
 double * copy_a=malloc((*N_good)*sizeof(double));
 // Compute median
 gsl_sort(a,1,(*N_good));
 median=gsl_stats_median_from_sorted_data(a,1,(*N_good));
 //MAD=gsl_stats_absdev_m(a,1,(*N_good),median);
 sigma=gsl_stats_sd_m(a,1,(*N_good),median);
 for(i=0,j=0;i<(*N_good);i++){
  // taken from http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h.htm
  //M=0.6745*(a[i]-median)/MAD;
  M=(a[i]-median)/sigma;
  if( fabs(M)<3.5 ){
   copy_a[j]=a[i];
   j++;
  }
 }
 for(i=0;i<j;i++)
  a[i]=copy_a[j];
 (*N_good)=j;
 free(copy_a);
 return;
}
*/

void remove_outliers_from_a_pair_of_arrays(double * a, double * b, int * N_good){
 int i,j;
 double median1,MAD1,M1;
 double median2,MAD2,M2;
 double * copy_a=malloc((*N_good)*sizeof(double));
 // Compute median
 gsl_sort(a,1,(*N_good));
 median1=gsl_stats_median_from_sorted_data(a,1,(*N_good));
 MAD1=gsl_stats_absdev_m(a,1,(*N_good),median1);
 gsl_sort(b,1,(*N_good));
 median2=gsl_stats_median_from_sorted_data(b,1,(*N_good));
 MAD2=gsl_stats_absdev_m(b,1,(*N_good),median2);
 for(i=0,j=0;i<(*N_good);i++){
  // taken from http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h.htm
  M1=0.6745*(a[i]-median1)/MAD1;
  M2=0.6745*(b[i]-median2)/MAD2;
  //fprintf(stderr,"%lf %lf  %lf %lf  %lf %lf\n",a[i]*3600,b[i]*3600,median1*3600,median2*3600,M1,M2);
  if( fabs(M1)<3.5 && fabs(M2)<3.5 ){
   copy_a[j]=a[i];
   j++;
  }
 }
 for(i=0;i<j;i++)
  a[i]=copy_a[j];
 (*N_good)=j;
 free(copy_a);
 return;
}


void set_catalog_search_parameters(double approximate_field_of_view_arcmin, struct str_catalog_search_parameters * catalog_search_parameters){
 catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP*approximate_field_of_view_arcmin/60.0;
 catalog_search_parameters->brightest_mag=1.0;
 catalog_search_parameters->faintest_mag=9.0;
 if( approximate_field_of_view_arcmin<500.0 ){
  catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP*approximate_field_of_view_arcmin/60.0;
  catalog_search_parameters->brightest_mag=2.0;
  catalog_search_parameters->faintest_mag=12.0; 
 }
 if( approximate_field_of_view_arcmin<400.0 ){
  catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP*approximate_field_of_view_arcmin/60.0;
  catalog_search_parameters->brightest_mag=5.0;
  catalog_search_parameters->faintest_mag=13.0; 
 }
 if( approximate_field_of_view_arcmin<240.0 ){
  catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP*approximate_field_of_view_arcmin/60.0;
  catalog_search_parameters->brightest_mag=6.0;
  catalog_search_parameters->faintest_mag=14.0; 
 }
 if( approximate_field_of_view_arcmin<120.0 ){
  catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP;
  catalog_search_parameters->brightest_mag=9.0;
  catalog_search_parameters->faintest_mag=15.0; 
 }
 if( approximate_field_of_view_arcmin<60.0 ){
  catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP;
  catalog_search_parameters->brightest_mag=10.0;
  catalog_search_parameters->faintest_mag=17.0; 
 }
 if( approximate_field_of_view_arcmin<30.0 ){
  catalog_search_parameters->search_radius_deg=MAX_DEVIATION_AT_FIRST_STEP;
  catalog_search_parameters->brightest_mag=12.0;
  catalog_search_parameters->faintest_mag=20.0; 
 }
 // HST
 if( approximate_field_of_view_arcmin<5.0 ){
  catalog_search_parameters->search_radius_deg=1.5/3600.0;
  catalog_search_parameters->brightest_mag=12.0;
  catalog_search_parameters->faintest_mag=22.0; 
 }
 catalog_search_parameters->search_radius_second_step_deg=0.6*catalog_search_parameters->search_radius_deg;
 return;
}

int blind_plate_solve_with_astrometry_net(char * fits_image_filename,double approximate_field_of_view_arcmin){
 char cmdstr[2*FILENAME_LENGTH];
 sprintf(cmdstr,"util/wcs_image_calibration.sh %s %.1lf",fits_image_filename,approximate_field_of_view_arcmin);
 fprintf(stderr,"Trying to blindly solve the plate...\n%s\n",cmdstr);
 if( 0!=system(cmdstr) ){fprintf(stderr,"ERROR solving the plate!\n");return 1;}
 return 0;
}

void guess_wcs_catalog_filename( char * wcs_catalog_filename, char * fits_image_filename){
 char test_string[FILENAME_LENGTH];
 strcpy(test_string,basename(fits_image_filename));
 if( test_string[0]=='w' && test_string[1]=='c' && test_string[2]=='s' && test_string[3]=='_' )
  wcs_catalog_filename[0]='\0';
 else
  strncpy(wcs_catalog_filename,"wcs_",FILENAME_LENGTH);wcs_catalog_filename[FILENAME_LENGTH-1]='\0';
 strncat(wcs_catalog_filename,test_string,FILENAME_LENGTH-strlen(wcs_catalog_filename));wcs_catalog_filename[FILENAME_LENGTH-1]='\0';
 strncat(wcs_catalog_filename,".cat",FILENAME_LENGTH-strlen(wcs_catalog_filename));wcs_catalog_filename[FILENAME_LENGTH-1]='\0';
// fprintf(stderr,"test_string=#%s#\n",test_string);
// fprintf(stderr,"fits_image_filename=#%s#\n",fits_image_filename);
// fprintf(stderr,"wcs_catalog_filename=#%s#\n",wcs_catalog_filename);
 return;
}

void write_wcs_catalog( char * fits_image_filename, struct detected_star * stars, int number_of_stars_in_wcs_catalog){
 FILE *f;
 int i;
 char wcs_catalog_filename[FILENAME_LENGTH];
 guess_wcs_catalog_filename( wcs_catalog_filename, fits_image_filename);
 strcat(wcs_catalog_filename,".ucac4");
 f=fopen(wcs_catalog_filename,"w");
 for(i=0;i<number_of_stars_in_wcs_catalog;i++){
  //if( stars[i].matched_with_catalog==1 && stars[i].good_star==1 )
   fprintf(f,"%6d %11.7lf %+10.7lf %9.4lf %9.4lf  %9.1lf %9.1lf %+7.4lf %6.4lf %3d    %6.3lf %5.3lf  %6.3lf %5.3lf  %6.3lf %5.3lf  %6.3lf %5.3lf  %6.3lf %5.3lf  %6.3lf %5.3lf\n",stars[i].n_current_frame,stars[i].corrected_ra_local,stars[i].corrected_dec_local,stars[i].x_pix,stars[i].y_pix,stars[i].flux,stars[i].flux_err,stars[i].mag,stars[i].mag_err,stars[i].flag,  stars[i].catalog_mag,stars[i].catalog_mag_err,  stars[i].APASS_B,stars[i].APASS_B_err,  stars[i].APASS_V,stars[i].APASS_V_err, stars[i].APASS_r,stars[i].APASS_r_err, stars[i].APASS_i,stars[i].APASS_i_err, stars[i].Rc_computed_from_APASS_ri,stars[i].Rc_computed_from_APASS_ri_err);
 }
 fclose(f);
 return;
}

void write_matched_stars_to_ds9_region( char * fits_image_filename, struct detected_star * stars, int number_of_stars_in_wcs_catalog){
 FILE *f;
 int i;
 char wcs_catalog_filename[FILENAME_LENGTH];
 guess_wcs_catalog_filename( wcs_catalog_filename, fits_image_filename);
 strcat(wcs_catalog_filename,".ds9.reg");
 f=fopen(wcs_catalog_filename,"w");
 fprintf(f,"# Region file format: DS9 version 4.0\n");
 fprintf(f,"# Filename:\n");
 fprintf(f,"global color=green font=\"sans 10 normal\" select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\n");
 fprintf(f,"image\n");    
 for(i=0;i<number_of_stars_in_wcs_catalog;i++)
  if(stars[i].matched_with_catalog==1)
   fprintf(f,"# text(%lf,%lf) text={%5.2lf}\ncircle(%lf,%lf,4.0)\n",stars[i].x_pix,stars[i].y_pix,(stars[i].catalog_dec-stars[i].corrected_dec_local)*3600,stars[i].x_pix,stars[i].y_pix);
 fclose(f);
 return;
}

int read_wcs_catalog( char * fits_image_filename, struct detected_star * stars, int * number_of_stars_in_wcs_catalog){
 FILE *f;
 int i;
 char wcs_catalog_filename[FILENAME_LENGTH];
 double double_garbage;
 int int_garbage;
 char char_garbage[4096];
 double X_im_size;
 double Y_im_size;
 gettime( fits_image_filename, &double_garbage, &int_garbage, 0, &X_im_size, &Y_im_size, char_garbage, char_garbage, 0, 0); // This is just an overkill way to get X_im_size Y_im_size
 guess_wcs_catalog_filename( wcs_catalog_filename, fits_image_filename);
 fprintf(stderr,"WCS catalog name: %s \n",wcs_catalog_filename);
 f=fopen(wcs_catalog_filename,"r");
 if( f==NULL )return 1;
 i=0;
 while(0<fscanf(f,"%d %lf %lf %lf %lf  %lf %lf %lf %lf %d",&stars[i].n_current_frame,&stars[i].ra_deg_measured,&stars[i].dec_deg_measured,&stars[i].x_pix,&stars[i].y_pix,&stars[i].flux,&stars[i].flux_err,&stars[i].mag,&stars[i].mag_err,&stars[i].flag)){
  ///
  if( stars[i].flux==0 )continue;   
  if( stars[i].flux_err==999999 )continue;
  if( stars[i].mag==99.0000 )continue;
  if( stars[i].mag_err==99.0000 )continue;
  if( stars[i].flux<MIN_SNR*stars[i].flux_err )continue;
  /// 
  if( stars[i].flag==0 )
   stars[i].good_star=1;
  else
   stars[i].good_star=0;
  //
  stars[i].ra_deg_measured_orig=stars[i].ra_deg_measured;
  stars[i].dec_deg_measured_orig=stars[i].dec_deg_measured;
  // set default values of the derived parameters
  stars[i].matched_with_catalog=0; // no catalog match at first
  stars[i].distance_from_image_edge=MIN(stars[i].x_pix,stars[i].y_pix);
  stars[i].distance_from_image_edge=MIN(stars[i].distance_from_image_edge,X_im_size-stars[i].x_pix);
  stars[i].distance_from_image_edge=MIN(stars[i].distance_from_image_edge,Y_im_size-stars[i].y_pix);
  if( stars[i].distance_from_image_edge<FRAME_EDGE_INDENT_PIXELS )stars[i].good_star=0;
  i++;
  if( i>=MAX_NUMBER_OF_STARS ){fprintf(stderr,"ERROR: too many stars in the SExtractor catalog file %s\n",wcs_catalog_filename);fclose(f);return 1;}
 }
 (*number_of_stars_in_wcs_catalog)=i;
 fclose(f);
 fprintf(stderr,"Got %d stars from the SExtractor catalog %s \n",i,wcs_catalog_filename);
 return 0;
}

int parse_string_with_APASS_magnitudes(char * string_with_APASS_magnitudes, double * APASS_B, double * APASS_B_err, double * APASS_V, double * APASS_V_err, double * APASS_r, double * APASS_r_err, double * APASS_i, double * APASS_i_err){
 char str[32];
 // check if the string is not too long
 if( 1023<strlen(string_with_APASS_magnitudes) )return 1;
 if( 9>strlen(string_with_APASS_magnitudes) )return 1;
// fprintf(stderr,"#%s#\n",string_with_APASS_magnitudes);
 str[0]=string_with_APASS_magnitudes[0];
 str[1]=string_with_APASS_magnitudes[1];
 str[2]=string_with_APASS_magnitudes[2];
 str[3]=string_with_APASS_magnitudes[3];
 str[4]=string_with_APASS_magnitudes[4];
 str[5]=string_with_APASS_magnitudes[5];
 str[6]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[5]==' ' )return 1;
 (*APASS_B)=atof(str);
// fprintf(stderr,"#%6.3lf#\n",(*APASS_B));
 if( (*APASS_B)>MAX_APASS_MAG || (*APASS_B)<MIN_APASS_MAG )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[7];
 str[1]=string_with_APASS_magnitudes[8];
 str[2]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[1]==' ' )return 1;
 (*APASS_B_err)=atof(str)/100.0;
// fprintf(stderr,"#%6.3lf#\n",(*APASS_B_err));
 if( (*APASS_B_err)>MAX_APASS_MAG_ERR || (*APASS_B_err)<MIN_APASS_MAG_ERR )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[10];
 str[1]=string_with_APASS_magnitudes[11];
 str[2]=string_with_APASS_magnitudes[12];
 str[3]=string_with_APASS_magnitudes[13];
 str[4]=string_with_APASS_magnitudes[14];
 str[5]=string_with_APASS_magnitudes[15];
 str[6]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[5]==' ' )return 1;
 (*APASS_V)=atof(str);
// fprintf(stderr,"#%6.3lf#\n",(*APASS_V));
 if( (*APASS_V)>MAX_APASS_MAG || (*APASS_V)<MIN_APASS_MAG )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[17];
 str[1]=string_with_APASS_magnitudes[18];
 str[2]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[1]==' ' )return 1;
 (*APASS_V_err)=atof(str)/100.0;
// fprintf(stderr,"#%6.3lf#\n",(*APASS_V_err));
 if( (*APASS_V_err)>MAX_APASS_MAG_ERR || (*APASS_V_err)<MIN_APASS_MAG_ERR )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[20];
 str[1]=string_with_APASS_magnitudes[21];
 str[2]=string_with_APASS_magnitudes[22];
 str[3]=string_with_APASS_magnitudes[23];
 str[4]=string_with_APASS_magnitudes[24];
 str[5]=string_with_APASS_magnitudes[25];
 str[6]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[5]==' ' )return 1;
 (*APASS_r)=atof(str);
// fprintf(stderr,"#%6.3lf#\n",(*APASS_r));
 if( (*APASS_r)>MAX_APASS_MAG || (*APASS_r)<MIN_APASS_MAG )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[27];
 str[1]=string_with_APASS_magnitudes[28];
 str[2]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[1]==' ' )return 1;
 (*APASS_r_err)=atof(str)/100.0;
// fprintf(stderr,"#%6.3lf#\n",(*APASS_r_err));
 if( (*APASS_r_err)>MAX_APASS_MAG_ERR || (*APASS_r_err)<MIN_APASS_MAG_ERR )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[30];
 str[1]=string_with_APASS_magnitudes[31];
 str[2]=string_with_APASS_magnitudes[32];
 str[3]=string_with_APASS_magnitudes[33];
 str[4]=string_with_APASS_magnitudes[34];
 str[5]=string_with_APASS_magnitudes[35];
 str[6]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[5]==' ' )return 1;
 (*APASS_i)=atof(str);
// fprintf(stderr,"#%6.3lf#\n",(*APASS_i));
 if( (*APASS_i)>MAX_APASS_MAG || (*APASS_i)<MIN_APASS_MAG )return 1; // simple check
 str[0]=string_with_APASS_magnitudes[37];
 str[1]=string_with_APASS_magnitudes[38];
 str[2]='\0';
// fprintf(stderr,"#%s#\n",str);
 if( str[1]==' ' )return 1;
 (*APASS_i_err)=atof(str)/100.0;
// fprintf(stderr,"#%6.3lf#\n",(*APASS_i_err));
 if( (*APASS_i_err)>MAX_APASS_MAG_ERR || (*APASS_i_err)<MIN_APASS_MAG_ERR )return 1; // simple check
 return 0;
}


int read_UCAC4_from_vizquery( struct detected_star * stars, int N, char *vizquery_output_filename, struct str_catalog_search_parameters * catalog_search_parameters){
 FILE *f;
 char string[1024];
 char string_with_APASS_magnitudes[1024];
 int i;
 double measured_ra,measured_dec,distance,catalog_ra,catalog_dec,catalog_mag,catalog_mag_err;
 char ucac4id[32];
 double cos_delta;
 int N_stars_matched_with_catalog=0;
 
 double APASS_B;
 double APASS_B_err;
 double APASS_V;
 double APASS_V_err;
 double APASS_r;
 double APASS_r_err;
 double APASS_i;
 double APASS_i_err;
 
 f=fopen(vizquery_output_filename,"r");
 while(NULL!=fgets(string, 1024, f)){
  if( string[0]=='#' )continue;
  if( string[0]=='\n' )continue;
  if( string[0]=='-' )continue;
  if( string[0]=='_' )continue;
  if( string[0]==' ' )continue;
  if( string[0]!=' ' && string[0]!='0' && string[0]!='1' && string[0]!='2' && string[0]!='3' && string[0]!='4' && string[0]!='5' && string[0]!='6' && string[0]!='7' && string[0]!='8' && string[0]!='9' )continue;
  
  string_with_APASS_magnitudes[0]='\0'; // reset the string with APASS magnitudes
  // We expect that all UCAC4-specific data including UCAC magnitude error are available  
  //if( 8>sscanf(string,"%lf %lf %lf %s %lf %lf %lf %lf %39[^\t\n]",&distance,&measured_ra,&measured_dec,ucac4id,&catalog_ra,&catalog_dec,&catalog_mag,&catalog_mag_err,string_with_APASS_magnitudes) )continue;
  // Evil CDS people changed the column order in vizquery output AGAIN!!!!
  if( 8>sscanf(string,"%lf %lf %lf %s %lf %lf %lf %lf %39[^\t\n]",&measured_ra,&measured_dec,&distance,ucac4id,&catalog_ra,&catalog_dec,&catalog_mag,&catalog_mag_err,string_with_APASS_magnitudes) )continue;
  if( catalog_mag_err>MAX_APASS_MAG_ERR )continue; // make sure APASS_B is not read instead of catalog_mag_err
  if( 0!=parse_string_with_APASS_magnitudes(string_with_APASS_magnitudes, &APASS_B, &APASS_B_err, &APASS_V, &APASS_V_err, &APASS_r, &APASS_r_err, &APASS_i, &APASS_i_err) )continue;
  
  cos_delta=cos(measured_dec*M_PI/180.0);
  for(i=0;i<N;i++){
   if( stars[i].matched_with_catalog==1 )continue;
   if( fabs(stars[i].dec_deg_measured-measured_dec)<catalog_search_parameters->search_radius_deg )
    if( fabs(stars[i].ra_deg_measured-measured_ra)*cos_delta<catalog_search_parameters->search_radius_deg ){
     if( distance>catalog_search_parameters->search_radius_deg*3600 )continue;
     stars[i].matched_with_catalog=1;
     stars[i].d_ra=catalog_ra-measured_ra;
     stars[i].d_dec=catalog_dec-measured_dec;
     stars[i].catalog_ra=catalog_ra;
     stars[i].catalog_dec=catalog_dec;
     stars[i].catalog_mag=catalog_mag;
     stars[i].catalog_mag_err=catalog_mag_err;
     strncpy(stars[i].ucac4id,ucac4id,32);stars[i].ucac4id[32-1]='\0';
     
     stars[i].APASS_B=APASS_B;
     stars[i].APASS_B_err=APASS_B_err;
     stars[i].APASS_V=APASS_V;
     stars[i].APASS_V_err=APASS_V_err;
     stars[i].APASS_r=APASS_r;
     stars[i].APASS_r_err=APASS_r_err;
     stars[i].APASS_i=APASS_i;
     stars[i].APASS_i_err=APASS_i_err;
     
     // Jester et al. (2005)
     // All stars with Rc-Ic < 1.15    
     // V-R    =    1.09*(r-i) + 0.22        0.03
     stars[i].Rc_computed_from_APASS_ri=APASS_V - 1.09*(APASS_r-APASS_i) - 0.22;
     stars[i].Rc_computed_from_APASS_ri_err=sqrt( 1.09*1.09*(APASS_r_err*APASS_r_err+APASS_i_err*APASS_i_err) + 0.03*0.03 );
          
     N_stars_matched_with_catalog++;
     //fprintf(stderr,"DEBUG MATCHED: stars[i].x_pix= %8.3lf\n",stars[i].x_pix);
    }
   //if( fabs(stars[i].dec_deg_measured-measured_dec)<MAX_DEVIATION_AT_FIRST_STEP )
  }//for(i=0;i<N;i++)
 }
 fclose(f);
 fprintf(stderr,"Matched %d stars with UCAC4.\n",N_stars_matched_with_catalog);
 if( N_stars_matched_with_catalog<10 ){
  fprintf(stderr,"ERROR: too few stars matched!\n");
  return 1;
 }
 return 0;
}

int usno_hack_read_USNOB1_from_vizquery( struct detected_star * stars, int N, char *vizquery_output_filename, struct str_catalog_search_parameters * catalog_search_parameters){
 FILE *f;
 char string[1024];
 //char string_with_APASS_magnitudes[1024];
 int i;
 double measured_ra,measured_dec,distance,catalog_ra,catalog_dec;//,catalog_mag,catalog_mag_err;
 char ucac4id[32];
 double cos_delta;
 int N_stars_matched_with_catalog=0;
 
 double APASS_B;
 double APASS_B_err=0.1;
 double APASS_V;
 double APASS_V_err=0.1;
 double APASS_r;
 double APASS_r_err=0.1;
 double APASS_i;
 double APASS_i_err=0.1;

 double USNO_B1mag;
 double USNO_R1mag;
 double USNO_B2mag;
 double USNO_R2mag;

 for(i=0;i<N;i++)stars[i].matched_with_catalog=0; // reset any possible previous UCAC match
 
 f=fopen(vizquery_output_filename,"r");
 while(NULL!=fgets(string, 1024, f)){
  if( string[0]=='#' )continue;
  if( string[0]=='\n' )continue;
  if( string[0]=='-' )continue;
  if( string[0]=='_' )continue;
  if( string[0]==' ' )continue;
  if( string[0]!=' ' && string[0]!='0' && string[0]!='1' && string[0]!='2' && string[0]!='3' && string[0]!='4' && string[0]!='5' && string[0]!='6' && string[0]!='7' && string[0]!='8' && string[0]!='9' )continue;
  
  //string_with_APASS_magnitudes[0]='\0'; // reset the string with APASS magnitudes
  // We expect that all UCAC4-specific data including UCAC magnitude error are available  
  //if( 8>sscanf(string,"%lf %lf %lf %s %lf %lf %lf %lf %39[^\t\n]",&distance,&measured_ra,&measured_dec,ucac4id,&catalog_ra,&catalog_dec,&catalog_mag,&catalog_mag_err,string_with_APASS_magnitudes) )continue;
  // Evil CDS people changed the column order in vizquery output AGAIN!!!!
  if( 10>sscanf(string,"%lf %lf %lf %s %lf %lf %lf %lf %lf %lf",&measured_ra,&measured_dec,&distance,ucac4id,&catalog_ra,&catalog_dec,&USNO_B1mag,&USNO_R1mag,&USNO_B2mag,&USNO_R2mag) )continue;
  
  cos_delta=cos(measured_dec*M_PI/180.0);
  for(i=0;i<N;i++){
   if( stars[i].matched_with_catalog==1 )continue;
   if( fabs(stars[i].dec_deg_measured-measured_dec)<catalog_search_parameters->search_radius_deg ){
    if( fabs(stars[i].ra_deg_measured-measured_ra)*cos_delta<catalog_search_parameters->search_radius_deg ){
     if( distance>catalog_search_parameters->search_radius_deg*3600 )continue;
     stars[i].matched_with_catalog=1;
     stars[i].d_ra=catalog_ra-measured_ra;
     stars[i].d_dec=catalog_dec-measured_dec;
     stars[i].catalog_ra=catalog_ra;
     stars[i].catalog_dec=catalog_dec;
     stars[i].catalog_mag=USNO_B2mag;
     stars[i].catalog_mag_err=0.1;
     strncpy(stars[i].ucac4id,ucac4id,32);stars[i].ucac4id[32-1]='\0';
     
     APASS_B=(USNO_B1mag+USNO_B2mag)/2.0;
     APASS_V= 0.444*USNO_B1mag + 0.556*USNO_R1mag; // John Greaves, see http://www.aerith.net/astro/color_conversion/JG/USNO-B1.0.html
     stars[i].Rc_computed_from_APASS_ri=APASS_i=APASS_r=(USNO_R1mag+USNO_R2mag)/2.0;
     stars[i].Rc_computed_from_APASS_ri_err=APASS_i_err=APASS_r_err=0.1;

     stars[i].APASS_B=APASS_B;
     stars[i].APASS_B_err=APASS_B_err;
     stars[i].APASS_V=APASS_V;
     stars[i].APASS_V_err=APASS_V_err;
     stars[i].APASS_r=APASS_r;
     stars[i].APASS_r_err=APASS_r_err;
     stars[i].APASS_i=APASS_i;
     stars[i].APASS_i_err=APASS_i_err;
     
     N_stars_matched_with_catalog++;
    }
   }
   //else{
   // fprintf(stderr,"DEBUG MAX_DEVIATION_AT_FIRST_STEP: stars[i].x_pix= %8.3lf\n",stars[i].x_pix);
   //}//if( fabs(stars[i].dec_deg_measured-measured_dec)<MAX_DEVIATION_AT_FIRST_STEP )
  }//for(i=0;i<N;i++)
 }
 fclose(f);
 fprintf(stderr,"Matched %d stars with USNO-B1.0.\n",N_stars_matched_with_catalog);
 if( N_stars_matched_with_catalog<10 ){
  fprintf(stderr,"ERROR: too few stars matched!\n");
  return 1;
 }
 return 0;
}

void usno_hack_create_vizquery_input_from_ucac4_output( char *vizquery_output_filename, char *vizquery_input_filename){
 FILE *f;
 FILE *f2;
 char string[1024];
 char string_with_APASS_magnitudes[1024];
 //int i;
 double measured_ra,measured_dec,distance,catalog_ra,catalog_dec,catalog_mag,catalog_mag_err;
 char ucac4id[32];
 //double cos_delta;
 //int N_stars_matched_with_catalog=0;
 f2=fopen(vizquery_input_filename,"w");
 f=fopen(vizquery_output_filename,"r");
 while(NULL!=fgets(string, 1024, f)){
  if( string[0]=='#' )continue;
  if( string[0]=='\n' )continue;
  if( string[0]=='-' )continue;
  if( string[0]=='_' )continue;
  if( string[0]==' ' )continue;
  if( string[0]!=' ' && string[0]!='0' && string[0]!='1' && string[0]!='2' && string[0]!='3' && string[0]!='4' && string[0]!='5' && string[0]!='6' && string[0]!='7' && string[0]!='8' && string[0]!='9' )continue;
  
  string_with_APASS_magnitudes[0]='\0'; // reset the string with APASS magnitudes
  // We expect that all UCAC4-specific data including UCAC magnitude error are available  
  //if( 8>sscanf(string,"%lf %lf %lf %s %lf %lf %lf %lf %39[^\t\n]",&distance,&measured_ra,&measured_dec,ucac4id,&catalog_ra,&catalog_dec,&catalog_mag,&catalog_mag_err,string_with_APASS_magnitudes) )continue;
  // Evil CDS people changed the column order in vizquery output AGAIN!!!!
  if( 8>sscanf(string,"%lf %lf %lf %s %lf %lf %lf %lf %39[^\t\n]",&measured_ra,&measured_dec,&distance,ucac4id,&catalog_ra,&catalog_dec,&catalog_mag,&catalog_mag_err,string_with_APASS_magnitudes) )continue;
  if( catalog_mag_err>MAX_APASS_MAG_ERR )continue; // make sure APASS_B is not read instead of catalog_mag_err
  //fprintf(f2,"%11.7lf %11.7lf\n",catalog_ra,catalog_dec);
  fprintf(f2,"%11.7lf %11.7lf\n",measured_ra,measured_dec);
 }
 fclose(f);
 fclose(f2);
 return;
}


int search_UCAC4_with_vizquery( struct detected_star * stars, int N, struct str_catalog_search_parameters * catalog_search_parameters){
 char command[1024];
 FILE *vizquery_input;
 int i;
 int pid=getpid();
 char vizquery_input_filename[FILENAME_LENGTH];
 char vizquery_output_filename[FILENAME_LENGTH];
 int vizquery_run_success;
 sprintf(vizquery_input_filename,"vizquery_%d.input",pid);
 sprintf(vizquery_output_filename,"vizquery_%d.output",pid);
 vizquery_input=fopen(vizquery_input_filename,"w");
 for(i=0;i<MIN(N,6000);i++){
  if( stars[i].good_star==1 ){
   fprintf(vizquery_input,"%lf %lf\n",stars[i].ra_deg_measured,stars[i].dec_deg_measured);
   //fprintf(stderr,"DEBUG  %lf %lf\n",stars[i].ra_deg_measured,stars[i].dec_deg_measured);
  }
 }
 fclose(vizquery_input);
 
 
 // yes, sorting in magnitude works
 sprintf(command, "`lib/find_timeout_command.sh` 180 lib/vizquery -site=%s -mime=text -source=UCAC4 -out.max=1 -out.add=_1 -out.add=_r -out.form=mini -out=UCAC4,RAJ2000,DEJ2000,a.mag,e_a.mag,Bmag,e_Bmag,Vmag,e_Vmag,rmag,e_rmag,imag,e_imag a.mag=%.1lf..%.1lf -sort=a.mag -c.rs=%.1lf -list=%s > %s",VIZIER_SITE,catalog_search_parameters->brightest_mag,catalog_search_parameters->faintest_mag,catalog_search_parameters->search_radius_deg*3600,vizquery_input_filename,vizquery_output_filename);
 
 fprintf(stderr, "%s\n", command);
 vizquery_run_success=system(command);
 if(vizquery_run_success!=0){
  fprintf(stderr,"WARNING: some problem running lib/vizquery script. Is this an internet connection problem? Retrying...\n");
  fprintf(stderr, "%s\n", command);
  vizquery_run_success=system(command);
  if(vizquery_run_success!=0){
   fprintf(stderr,"WARNING: some problem running lib/vizquery script. Is this an internet connection problem? Retrying...\n");
   fprintf(stderr, "%s\n", command);
   vizquery_run_success=system(command);
   if(vizquery_run_success!=0){fprintf(stderr,"ERROR: problem running lib/vizquery script :(\n");exit(1);}
  }
 }

 if(0!=read_UCAC4_from_vizquery( stars, N, vizquery_output_filename, catalog_search_parameters)){
  fprintf(stderr,"Problem getting data from VizieR. :(\n");
  fprintf(stderr,"Maybe this sky area is not covered by APASS yet?\nTrying the USNO-B1.0 hack...\n\nWARNING: using USNO-B1.0 instead of UCAC4!!!!\n\n");
  usno_hack_create_vizquery_input_from_ucac4_output( vizquery_output_filename, vizquery_input_filename);
  //sprintf(command, "lib/vizquery -site=vizier.u-strasbg.fr -mime=text -source=USNO-B1.0 -out.max=1 -out.add=_1 -out.add=_r -out.form=mini -out=USNO-B1.0,RAJ2000,DEJ2000,B1mag,R1mag,B2mag,R2mag -sort=B2mag B2mag=%.1lf..%.1lf B1mag=%.1lf..%.1lf R2mag=%.1lf..%.1lf -c.rs=%.1lf -list=%s > %s",MIN_APASS_MAG,MAX_APASS_MAG,MIN_APASS_MAG,MAX_APASS_MAG,catalog_search_parameters->brightest_mag,catalog_search_parameters->faintest_mag, catalog_search_parameters->search_radius_deg*3600,vizquery_input_filename,vizquery_output_filename);
  sprintf(command, "`lib/find_timeout_command.sh` 180 lib/vizquery -site=%s -mime=text -source=USNO-B1.0 -out.max=1 -out.add=_1 -out.add=_r -out.form=mini -out=USNO-B1.0,RAJ2000,DEJ2000,B1mag,R1mag,B2mag,R2mag -sort=B2mag B2mag=%.1lf..%.1lf B1mag=%.1lf..%.1lf R2mag=%.1lf..%.1lf -c.rs=%.1lf -list=%s > %s",VIZIER_SITE,MIN_APASS_MAG,MAX_APASS_MAG,MIN_APASS_MAG,MAX_APASS_MAG,catalog_search_parameters->brightest_mag,catalog_search_parameters->faintest_mag, catalog_search_parameters->search_radius_deg*3600,vizquery_input_filename,vizquery_output_filename);
  fprintf(stderr, "%s\n", command);
  vizquery_run_success=system(command);
  if(vizquery_run_success!=0){
   fprintf(stderr,"WARNING: some problem running lib/vizquery script. Is this an internet connection problem? Retrying...\n");
   fprintf(stderr, "%s\n", command);
   vizquery_run_success=system(command);
   if(vizquery_run_success!=0){
    fprintf(stderr,"WARNING: some problem running lib/vizquery script. Is this an internet connection problem? Retrying...\n");
    fprintf(stderr, "%s\n", command);
    vizquery_run_success=system(command);
    if(vizquery_run_success!=0){fprintf(stderr,"ERROR: problem running lib/vizquery script :(\n");exit(1);}
   }
  }
  if(0!=usno_hack_read_USNOB1_from_vizquery( stars, N, vizquery_output_filename, catalog_search_parameters)){fprintf(stderr,"Problem getting data from VizieR. :(\n");return 1;}
 }
 // delete temporary files only on success
 if( vizquery_run_success==0 ){
  if(0!=unlink(vizquery_input_filename))fprintf(stderr,"WARNING! Cannot delete temporary file %s\n",vizquery_input_filename);
  if(0!=unlink(vizquery_output_filename))fprintf(stderr,"WARNING! Cannot delete temporary file %s\n",vizquery_output_filename);
 }
 
 return 0;
}

static inline int compare_star_on_mag_solve(const void *a1, const void *a2) {
 struct detected_star *s1, *s2;
 s1 = (struct detected_star *)a1;
 s2 = (struct detected_star *)a2;
 if( s1->mag<s2->mag )
  return 0;
 else
  return 1;
}

static inline double compute_distance_on_sphere(double RA1, double DEC1, double target_ra, double target_dec){
 double distance;
 double deg2rad_conversion=M_PI/180.0;
 double RA1_rad=RA1*deg2rad_conversion;
 double DEC1_rad=DEC1*deg2rad_conversion;
 double target_ra_rad=target_ra*deg2rad_conversion;
 double target_dec_rad=target_dec*deg2rad_conversion;
 distance=deg2rad_conversion*acos(cos(DEC1_rad)*cos(target_dec_rad)*cos( MAX(RA1_rad,target_ra_rad)-MIN(RA1_rad,target_ra_rad) )+sin(DEC1_rad)*sin(target_dec_rad));
 return distance;
}

void correct_measured_positions( struct detected_star * stars, int N, double search_radius, int process_only_stars_matched_with_catalog, struct str_catalog_search_parameters * catalog_search_parameters){
 int i,j,N_good;

 double A1,B1,C1,A2,B2,C2;

 double *x;
 double *y;
 double *z1;
 double *z2;

 double distance;
 
 x=malloc(N*sizeof(double));
 y=malloc(N*sizeof(double));
 z1=malloc(N*sizeof(double));
 z2=malloc(N*sizeof(double));

 // *** Find global linear solution  ***
 for(i=0,N_good=0;i<N;i++){
  if( stars[i].good_star==1 && stars[i].matched_with_catalog==1 ){
//   x[N_good]=stars[i].ra_deg_measured;
//   y[N_good]=stars[i].dec_deg_measured;
   x[N_good]=stars[i].x_pix;
   y[N_good]=stars[i].y_pix;
   z1[N_good]=stars[i].d_ra;
   z2[N_good]=stars[i].d_dec;
   N_good++;
  }
 }
 
// for(i=0;i<N_good;i++)
//  fprintf(stderr,"OGAOGA %lf %lf %lf\n",x[i],y[i],z1[i]);
 
 fit_plane_lin( x, y, z1, (unsigned int)N_good, &A1, &B1, &C1);  
 fit_plane_lin( x, y, z2, (unsigned int)N_good, &A2, &B2, &C2);

 // apply the linear correction
 /*
 #ifdef VAST_ENABLE_OPENMP                            
  #ifdef _OPENMP                                      
   #pragma omp parallel for private(i)
  #endif
 #endif
 */
 for(i=0;i<N;i++){
  //stars[i].computed_d_ra=A1*stars[i].ra_deg_measured+B1*stars[i].dec_deg_measured+C1;
  stars[i].computed_d_ra=A1*stars[i].x_pix+B1*stars[i].y_pix+C1;
  //stars[i].computed_d_dec=A2*stars[i].ra_deg_measured+B2*stars[i].dec_deg_measured+C2;
  stars[i].computed_d_dec=A2*stars[i].x_pix+B2*stars[i].y_pix+C2;
  stars[i].corrected_ra=stars[i].ra_deg_measured+stars[i].computed_d_ra;
  stars[i].corrected_dec=stars[i].dec_deg_measured+stars[i].computed_d_dec;
  // clean-up outliers after the linear fit
  if( stars[i].matched_with_catalog==1 )
   if( compute_distance_on_sphere(stars[i].catalog_ra,stars[i].catalog_dec,stars[i].corrected_ra,stars[i].corrected_dec)>catalog_search_parameters->search_radius_second_step_deg )
    stars[i].matched_with_catalog=0;
  // fprintf(stderr,"corrected-catalog=%lf measured-catalog=%lf correction=%lf\n",(stars[i].corrected_dec-stars[i].catalog_dec)*3600,(stars[i].dec_deg_measured-stars[i].catalog_dec)*3600,stars[i].computed_d_dec*3600);
 }

 // *** Find astrometric correction as a function of magnitude  ***
 double poly_coeff[10];

 for(N_good=0,i=0;i<N;i++){
  if( stars[i].good_star==1 && stars[i].matched_with_catalog==1 ){
   x[N_good]=stars[i].mag;
   y[N_good]=0.1; // fake error, same for all stars since we want an unweighted fit
   z1[N_good]=stars[i].catalog_ra-stars[i].corrected_ra;
   z1[N_good]=stars[i].catalog_dec-stars[i].corrected_dec;
   N_good++;
  }
 }
 wlinearfit(x, z1, y, N_good, poly_coeff);
 C1=poly_coeff[0];
 A1=poly_coeff[1];
 wlinearfit(x, z2, y, N_good, poly_coeff);
 C2=poly_coeff[0];
 A2=poly_coeff[1];
 // apply the correction
 /*
 #ifdef VAST_ENABLE_OPENMP                            
  #ifdef _OPENMP                                      
   #pragma omp parallel for private(i)
  #endif
 #endif
 */
 for(i=0;i<N;i++){
  stars[i].corrected_mag_ra=stars[i].corrected_ra+(A1*stars[i].mag+C1);
  stars[i].corrected_mag_dec=stars[i].corrected_dec+(A2*stars[i].mag+C2);
 }

 // *** Find local corrections  ***
 double local_correction_ra;
 double local_correction_dec;
 //double target_mag;
 double target_ra;
 double target_dec;
 double target_x_pix;
 double target_y_pix;
 double current_accuracy, current_accuracy_ra, current_accuracy_dec;
 double best_accuracy;
 double current_search_radius;
 double best_search_radius=REFERENCE_LOCAL_SOLUTION_RADIUS_DEG;
 double best_local_correction_ra;
 double best_local_correction_dec;
 double cos_delta;

 // And finally make a correction using the few nearby stars
 // if this is not a global solution
 //if( target_ra!=0.0 && target_dec!=0.0 ){
/*
 #ifdef VAST_ENABLE_OPENMP                            
  #ifdef _OPENMP                                      
   #pragma omp parallel for private(j, target_x_pix, target_ra, target_dec, cos_delta, best_accuracy, best_search_radius, current_search_radius, i, distance, z1, z2, N_good, current_accuracy_ra, current_accuracy_dec, current_accuracy)
  #endif
 #endif
*/
 for(j=0;j<N;j++){
 
  if( process_only_stars_matched_with_catalog==1 && stars[j].matched_with_catalog!=1 )continue;
  // set target star
  //target_mag=stars[j].mag;
  target_x_pix=stars[j].x_pix;
  target_y_pix=stars[j].y_pix;
  target_ra=stars[j].ra_deg_measured;
  target_dec=stars[j].dec_deg_measured;
  cos_delta=cos(stars[j].dec_deg_measured*M_PI/180.0);

  // try various corrections
  best_accuracy=99.9;
  best_search_radius=99.9; // just so we can distinguish if the value comes from a previous star or not
  //for(current_search_radius=search_radius;current_search_radius>0.05*search_radius;current_search_radius=current_search_radius-0.1*search_radius){
  for(current_search_radius=search_radius;current_search_radius>0.05*search_radius;current_search_radius=current_search_radius-0.1*current_search_radius){
   // determine the best search radius for local correction
   //for(i=0,N_good=0;i<N;i++){
   N_good=0;
   for(i=N;i--;){
    if( stars[i].matched_with_catalog!=1 )continue;
    if( stars[i].good_star!=1 )continue;
    if( fabs(target_x_pix-stars[i].x_pix)>500 )continue; // a miserable attempt to optimize
    if( fabs(target_y_pix-stars[i].y_pix)>500 )continue; // a miserable attempt to optimize
    if( fabs(target_dec-stars[i].dec_deg_measured)>current_search_radius )continue; // a miserable attempt to optimize
    //
    //if( target_mag<stars[i].mag-1.5 )continue;
    ///
    distance=compute_distance_on_sphere(stars[i].ra_deg_measured, stars[i].dec_deg_measured, target_ra, target_dec);
    if( distance<current_search_radius ){
     //if( compute_distance_on_sphere(stars[i].catalog_ra,stars[i].catalog_dec,stars[i].corrected_ra,stars[i].corrected_dec)<catalog_search_parameters->search_radius_second_step_deg ){
      z1[N_good]=stars[i].catalog_ra-stars[i].corrected_mag_ra;
      z2[N_good]=stars[i].catalog_dec-stars[i].corrected_mag_dec;
      N_good++;
      if( N_good>501 )break; // too many stars
     //}
    }
   }
   if( N_good>500 )continue; // too many stars
   ///
   //remove_outliers_from_an_array(z1, &N_good);
   remove_outliers_from_a_pair_of_arrays( z1, z2, &N_good);
   //remove_outliers_from_an_array(z2, &N_good);
   ///
   if( N_good<5 )break; // too few stars
   current_accuracy_ra=gsl_stats_sd(z1,1,N_good);
   current_accuracy_dec=gsl_stats_sd(z2,1,N_good);
   current_accuracy=sqrt(current_accuracy_ra*cos_delta*current_accuracy_ra*cos_delta+current_accuracy_dec*current_accuracy_dec);
   if( current_accuracy<best_accuracy ){
    best_accuracy=current_accuracy;
    best_search_radius=current_search_radius;
   }
  }


  // determine the local correction using the best search radius
  for(i=0,N_good=0;i<N;i++){
   if( stars[i].matched_with_catalog!=1 )continue;
   if( stars[i].good_star!=1 )continue;
   if( fabs(target_x_pix-stars[i].x_pix)>500 )continue; // a miserable attempt to optimize
   if( fabs(target_y_pix-stars[i].y_pix)>500 )continue; // a miserable attempt to optimize
   if( fabs(target_dec-stars[i].dec_deg_measured)>best_search_radius )continue; // a miserable attempt to optimize
   //
   //if( target_mag<stars[i].mag-1.5 )continue;
   //
   distance=compute_distance_on_sphere(stars[i].ra_deg_measured, stars[i].dec_deg_measured, target_ra, target_dec);
   if( distance<best_search_radius ){
    //if( distance==0.0 )continue;//
    z1[N_good]=stars[i].catalog_ra-stars[i].corrected_mag_ra;
    z2[N_good]=stars[i].catalog_dec-stars[i].corrected_mag_dec;
    //if( z1[N_good]==0.0 || z2[N_good]==0.0 )continue;//
    N_good++;
   }
  }

  ///
  //remove_outliers_from_an_array(z1, &N_good);
  remove_outliers_from_a_pair_of_arrays( z1, z2, &N_good);
  //remove_outliers_from_an_array(z2, &N_good);
  ///

  if( N_good>=3 ){
   gsl_sort(z1,1,N_good);
   gsl_sort(z2,1,N_good);
   best_local_correction_ra=gsl_stats_median_from_sorted_data(z1,1,N_good);
   best_local_correction_dec=gsl_stats_median_from_sorted_data(z2,1,N_good);
  }
  else{
   //fprintf(stderr,"DEBUG: best_accuracy=%lf best_search_radius=%lf (arcmin)\n",best_accuracy,best_search_radius*60);
   best_local_correction_ra=0.0;
   best_local_correction_dec=0.0;
   best_accuracy=0.0;
  }

  //fprintf(stderr,"DEBUG: best_accuracy=%lf best_search_radius=%lf stars[j].matched_with_catalog=%d best_local_correction_ra=%lf best_local_correction_dec=%lf\n",best_accuracy*3600,best_search_radius,stars[j].matched_with_catalog,best_local_correction_ra*3600,best_local_correction_dec*3600);

  local_correction_ra=best_local_correction_ra;
  local_correction_dec=best_local_correction_dec;

  // apply the local correction
  stars[j].local_correction_ra=local_correction_ra;
  stars[j].local_correction_dec=local_correction_dec;
  stars[j].corrected_ra_local=stars[j].corrected_mag_ra+stars[j].local_correction_ra;
  stars[j].corrected_dec_local=stars[j].corrected_mag_dec+stars[j].local_correction_dec;
  stars[j].estimated_local_correction_accuracy=best_accuracy;
  
  // clean-up outliers
  if( stars[j].matched_with_catalog==1 ){
   distance=compute_distance_on_sphere(stars[j].corrected_ra_local, stars[j].corrected_dec_local, stars[j].catalog_ra, stars[j].catalog_dec);
   if( distance>catalog_search_parameters->search_radius_second_step_deg )
    stars[j].matched_with_catalog=0;
  }
  //

 }//for(j=0;j<N;j++){
 //}
 
 // Estimate accuracy
 for(i=0,j=0;j<N;j++){
  if( stars[j].estimated_local_correction_accuracy!=0.0 ){
   z1[i]=stars[j].estimated_local_correction_accuracy;
   i++;
  }
 }
 gsl_sort(z1,1,i);
 fprintf(stderr,"Estimated accuracy of the plate solution: %.2lf\"\n",3600*gsl_stats_median_from_sorted_data(z1,1,i));
 

 free(x);
 free(y);
 free(z1);
 free(z2);
 
 return;
}


int main(int argc, char **argv){
 int number_of_stars_in_wcs_catalog;
 struct detected_star *stars;
 char fits_image_filename[FILENAME_LENGTH];
 int i;
 stars=malloc(MAX_NUMBER_OF_STARS*sizeof(struct detected_star));  
 
 double approximate_field_of_view_arcmin;
 
 int stars_matched_at_previous_iteration,stars_matched_at_this_iteration;
 
 struct str_catalog_search_parameters catalog_search_parameters;

 int solution_iteration;
 
 FILE * pipe_for_try_to_guess_image_fov;
 char command_string[2*FILENAME_LENGTH];


 if( argc<2 || argc>3 ){
  fprintf(stderr,"Usage: %s myfitsimage.fits [APPROXIMATE_FIELD_OF_VIEW_ARCMIN]\n\nExamples: %s myfitsimage.fits\n          %s myfitsimage.fits 40\n",argv[0],argv[0],argv[0]);
  return 1;
 }
 

 strncpy(fits_image_filename,argv[1],FILENAME_LENGTH);fits_image_filename[FILENAME_LENGTH-1]='\0';
 // Test if whatever is provided as argv[1] is actually a readable FITS file
 if( 0!=fitsfile_read_check(fits_image_filename) ){fprintf(stderr,"ERROR reading FITS file %s\n",fits_image_filename);return 1;}

 // set the field of view, if needed
 if( argc==3 ){
  approximate_field_of_view_arcmin=atof(argv[2]);
  if( approximate_field_of_view_arcmin<2.0 || approximate_field_of_view_arcmin>600.0 ){
   fprintf(stderr,"Warining! The supplied approximate field of view (%s = %lf) arcmin seems wrong. Resorting to the default value of %lf.\n",argv[2],approximate_field_of_view_arcmin,DEFAULT_APPROXIMATE_FIELD_OF_VIEW_ARCMIN);
   approximate_field_of_view_arcmin=DEFAULT_APPROXIMATE_FIELD_OF_VIEW_ARCMIN;
  }
 }
 else{
  //approximate_field_of_view_arcmin=DEFAULT_APPROXIMATE_FIELD_OF_VIEW_ARCMIN;
  sprintf(command_string,"lib/try_to_guess_image_fov %s",fits_image_filename);
  pipe_for_try_to_guess_image_fov=popen(command_string, "r");
  if( NULL==pipe_for_try_to_guess_image_fov ){
   fprintf(stderr,"WARNING: failed to run command: %s\n",command_string);
   approximate_field_of_view_arcmin=DEFAULT_APPROXIMATE_FIELD_OF_VIEW_ARCMIN;
  }
  else{
   fscanf(pipe_for_try_to_guess_image_fov,"%lf",&approximate_field_of_view_arcmin);
   pclose(pipe_for_try_to_guess_image_fov);
  }
 }

 set_catalog_search_parameters(approximate_field_of_view_arcmin, &catalog_search_parameters);


 // **** Blind plate solution with Astrometry.net ****
 if( 0!=blind_plate_solve_with_astrometry_net(fits_image_filename,approximate_field_of_view_arcmin) ){fprintf(stderr,"ERROR: cannot perform blind plate solution.\n");return 1;}
 
 // **** Read the star catalog ****
 if( 0!=read_wcs_catalog( fits_image_filename, stars, &number_of_stars_in_wcs_catalog) ){fprintf(stderr,"ERROR: reading SExtractor catalog file...\n");return 1;}
 qsort(stars, number_of_stars_in_wcs_catalog, sizeof(struct detected_star), compare_star_on_mag_solve);
 
 // Make sure all stars have a flag that they are not matched with a catalog yet
 //for(i=0;i<number_of_stars_in_wcs_catalog;i++)stars[i].matched_with_catalog=0;
 
 // **** Querry UCAC4 ****
 fprintf(stderr,"\nITERATION 01\n");
 if( 0!=search_UCAC4_with_vizquery( stars, number_of_stars_in_wcs_catalog, &catalog_search_parameters) ){fprintf(stderr,"ERROR running vizquery\n");return 1;}

 // Check if there is any hope
 // compute the number of matched stars
 for(stars_matched_at_this_iteration=0,i=0;i<number_of_stars_in_wcs_catalog;i++){
  if(stars[i].matched_with_catalog==1)
   stars_matched_at_this_iteration++;
 }
 if( stars_matched_at_this_iteration<50 ){
  fprintf(stderr,"\nRetrying with a larger catalog search radius.\n");
  catalog_search_parameters.search_radius_deg=2.0*catalog_search_parameters.search_radius_deg;
  catalog_search_parameters.search_radius_second_step_deg=2.0*catalog_search_parameters.search_radius_second_step_deg;
  if( 0!=search_UCAC4_with_vizquery( stars, number_of_stars_in_wcs_catalog, &catalog_search_parameters) ){fprintf(stderr,"ERROR running vizquery\n");return 1;}
 }

 
 // Correct the measured positions
 fprintf(stderr,"Correcting the measured astrometric positions...\n");
 correct_measured_positions( stars, number_of_stars_in_wcs_catalog, REFERENCE_LOCAL_SOLUTION_RADIUS_DEG, 1, &catalog_search_parameters);
 
 fprintf(stderr,"Removing outliers...\n");
 // Remove outliers (likley wrong identifications)
 for(i=0;i<number_of_stars_in_wcs_catalog;i++)
  if(stars[i].matched_with_catalog==1)
   if( compute_distance_on_sphere(stars[i].catalog_ra,stars[i].catalog_dec,stars[i].corrected_ra_local,stars[i].corrected_dec_local)>catalog_search_parameters.search_radius_second_step_deg )stars[i].matched_with_catalog=0;
 
 // re-compute correct the measured positions
 fprintf(stderr,"Correcting the measured astrometric positions...\n");
 correct_measured_positions( stars, number_of_stars_in_wcs_catalog, REFERENCE_LOCAL_SOLUTION_RADIUS_DEG, 0, &catalog_search_parameters);
 
 // compute the number of matched stars
 for(stars_matched_at_this_iteration=0,i=0;i<number_of_stars_in_wcs_catalog;i++){
  if(stars[i].matched_with_catalog==1)
   stars_matched_at_this_iteration++;
 }
 stars_matched_at_previous_iteration=stars_matched_at_this_iteration;

 // Apply the local corrections and iterate to find a better plate solution
 for(solution_iteration=2;solution_iteration<=10;solution_iteration++){
 
  fprintf(stderr,"\nITERATION %02d\n",solution_iteration);
  for(i=0;i<number_of_stars_in_wcs_catalog;i++){
   stars[i].ra_deg_measured=stars[i].corrected_ra_local;
   stars[i].dec_deg_measured=stars[i].corrected_dec_local;
   // Make sure all stars have a flag that they are not matched with a catalog yet
   stars[i].matched_with_catalog=0;
  }
  //search_UCAC4_with_vizquery( stars, number_of_stars_in_wcs_catalog);
  if( 0!=search_UCAC4_with_vizquery( stars, number_of_stars_in_wcs_catalog, &catalog_search_parameters) ){fprintf(stderr,"ERROR running vizquery\n");return 1;}
  correct_measured_positions( stars, number_of_stars_in_wcs_catalog, REFERENCE_LOCAL_SOLUTION_RADIUS_DEG, 1, &catalog_search_parameters);
  // remove outliers
  for(i=0;i<number_of_stars_in_wcs_catalog;i++)
   if(stars[i].matched_with_catalog==1)
    if( compute_distance_on_sphere(stars[i].catalog_ra,stars[i].catalog_dec,stars[i].corrected_ra_local,stars[i].corrected_dec_local)>catalog_search_parameters.search_radius_second_step_deg )stars[i].matched_with_catalog=0;
  // re-compute the corrections, now applying them to all stars
  correct_measured_positions( stars, number_of_stars_in_wcs_catalog, REFERENCE_LOCAL_SOLUTION_RADIUS_DEG, 0, &catalog_search_parameters);
 
  // compute the number of matched stars
  stars_matched_at_previous_iteration=stars_matched_at_this_iteration;
  for(stars_matched_at_this_iteration=0,i=0;i<number_of_stars_in_wcs_catalog;i++){
   if(stars[i].matched_with_catalog==1)
    stars_matched_at_this_iteration++;
  }


  // Report the current status
  fprintf(stderr,"Excluding outliers - %d stars left (next iteration limit %d)\n", stars_matched_at_this_iteration,stars_matched_at_previous_iteration+(int)(0.1*stars_matched_at_previous_iteration) );
 
  // check if there was a noticable improvement in the solution
  if( stars_matched_at_this_iteration<stars_matched_at_previous_iteration+(int)(0.1*stars_matched_at_previous_iteration) ){
   fprintf(stderr,"Stop iterations.\n");
   if( stars_matched_at_this_iteration<MIN_NUMBER_OF_STARS_FOR_UCAC4_MATCH ){
    fprintf(stderr,"\n\n The number of stars matched with the catalog is suspiciously low!\n Something is not right here... :(\n\n");
   }
   else{
    fprintf(stderr,"\n\n The field is successfully solved and matched with the catalog! :)\n\n");
   }
   break;
  }
 
 }


/*
 // Inspect the output
 FILE * solve_plate_debug;
 solve_plate_debug=fopen("solve_plate_debug.txt","w");
 for(i=0;i<number_of_stars_in_wcs_catalog;i++){
  if(stars[i].matched_with_catalog==1)
   fprintf(solve_plate_debug,"%10lf %10lf  %+10lf %+10lf  %+10lf %+10lf  %+10lf %+10lf %+10lf %+10lf  %+10lf %+10lf   %+10lf %+10lf %+10lf %+10lf  %+10lf  %+10lf %+10lf\n",\
// 1 2
   stars[i].x_pix,stars[i].y_pix,\
// 3 4
   stars[i].d_ra*3600,stars[i].d_dec*3600,\
// 5 6
   stars[i].local_correction_ra*3600, stars[i].local_correction_dec*3600,\
// 7 8
   (stars[i].catalog_ra-stars[i].corrected_ra_local)*3600, (stars[i].catalog_dec-stars[i].corrected_dec_local)*3600, \
// 9 10
   (stars[i].catalog_ra-stars[i].ra_deg_measured_orig)*3600, (stars[i].catalog_dec-stars[i].dec_deg_measured_orig)*3600, \
// 11 12
   (stars[i].catalog_ra-stars[i].corrected_ra)*3600, (stars[i].catalog_dec-stars[i].corrected_dec)*3600, \
// 13 14
   stars[i].catalog_ra, stars[i].catalog_dec, \
// 15 16
   stars[i].computed_d_ra*3600, stars[i].computed_d_dec*3600, \
// 17
   stars[i].mag, \
   (stars[i].catalog_ra-stars[i].corrected_mag_ra)*3600, (stars[i].catalog_dec-stars[i].corrected_mag_dec)*3600 );
 }
 fclose(solve_plate_debug);
*/
 
 // Write output
 write_wcs_catalog( fits_image_filename, stars, number_of_stars_in_wcs_catalog);
 write_matched_stars_to_ds9_region( fits_image_filename, stars, number_of_stars_in_wcs_catalog);
 
 free(stars);
 return 0;
}