src/gal_bkgd.c

   1 /*
   2   gal_bkgd.c
   3   sky_model
   4
   5   Created by Jayant Murthy on 07/10/12.
   6   Copyright (c) 2012 Jayant Murthy. All rights reserved.
   7 */
   8 #include "sky_model.h"
   9
  10 struct BKGD_DATA READ_GAL_BKGD(char bkgd_file[FLEN_FILENAME]) {
  11
  12     FILE *file_ptr;
  13     int data_index;
  14     double val, eval, t1, t2;
  15     struct BKGD_DATA gal_bkgd;
  16
  17     file_ptr = fopen(bkgd_file, "r"); /*Open input text file*/
  18     if (file_ptr == NULL){
  19         printf("Could not find the file %s\n", bkgd_file);
  20         exit(EXIT_FAILURE);
  21     }
  22
  23         /*Read the data. The maximum number of elements is set by MAX_BKGD_ELEM*/
  24         data_index = 0;
  25         gal_bkgd.ra   = malloc(sizeof(double) * MAX_BKGD_ELEM);
  26         gal_bkgd.dec  = malloc(sizeof(double) * MAX_BKGD_ELEM);
  27         gal_bkgd.x       = malloc(sizeof(double) * MAX_BKGD_ELEM);
  28         gal_bkgd.y       = malloc(sizeof(double) * MAX_BKGD_ELEM);
  29         gal_bkgd.z       = malloc(sizeof(double) * MAX_BKGD_ELEM);
  30         gal_bkgd.bkgd = malloc(sizeof(double) * MAX_BKGD_ELEM);
  31         gal_bkgd.err  = malloc(sizeof(double) * MAX_BKGD_ELEM);
  32
  33     /*Read UV background file*/
  34     while (feof(file_ptr) == 0){ /*Continue until end of input file*/
  35                 fscanf(file_ptr, "%lf %lf %lf %lf", &t1, &t2, &val, &eval);
  36         /*Galactic background is in galactic coordinates*/
  37                 gal_bkgd.ra[data_index] = t1;
  38                 gal_bkgd.dec[data_index] = t2;
  39         /*Convert to cartesian coordinates*/
  40                 gal_bkgd.x[data_index] =
  41         cos(DEGRAD(gal_bkgd.ra[data_index]))*
  42         cos(DEGRAD(gal_bkgd.dec[data_index]));
  43                 gal_bkgd.y[data_index] =
  44         sin(DEGRAD(gal_bkgd.ra[data_index]))*
  45         cos(DEGRAD(gal_bkgd.dec[data_index]));
  46                 gal_bkgd.z[data_index] =
  47         sin(DEGRAD(gal_bkgd.dec[data_index]));
  48         /*Background value at that position and the error*/
  49                 gal_bkgd.bkgd[data_index] = val;
  50                 gal_bkgd.err[data_index] = eval;
  51                 ++data_index;
  52         }/*End read file*/
  53     gal_bkgd.npoints = data_index;
  54     return(gal_bkgd);
  55
  56 }
  57
  58 int CALC_BKGD_FLUX(struct WCS_HEAD wcs_in, long naxes[2], float *data, float *var,
  59                                    float ang_limit, float csc_law,
  60                    struct BKGD_DATA bkgd)
  61 {
  62
  63         double ra0, dec0, x0, y0, z0;
  64         double val, max_data=10000;
  65         double cang_limit, dist, cm_angle;
  66         int status = 0, data_index, index;
  67         int i, j;
  68     int *inc_list, inc_index = 0;
  69
  70
  71         /*ang_limit is the maximum angle for binning data but it is easier to work with
  72          the cosine*/
  73         cang_limit = cos(DEGRAD(ang_limit));
  74     data_index = bkgd.npoints;
  75
  76     /*I don't have to search the huge dataset each time
  77       so I'll restrict the search range.*/
  78     inc_list = (int *) malloc(sizeof(int)*data_index);
  79     /*Go overboard on the limiting size*/
  80     cm_angle = cos(DEGRAD(wcs_in.angsize*1.5 + ang_limit));
  81
  82     for (index = 0; index < data_index; ++index) {
  83         ra0 = wcs_in.xrefval;
  84         dec0 = wcs_in.yrefval;
  85         x0 = cos(DEGRAD(ra0))*cos(DEGRAD(dec0));
  86         y0 = sin(DEGRAD(ra0))*cos(DEGRAD(dec0));
  87         z0 = sin(DEGRAD(dec0));
  88         dist = bkgd.x[index]*x0 + bkgd.y[index]*y0 + bkgd.z[index]*z0;
  89         if (dist > cm_angle) {
  90             inc_list[inc_index] = index;
  91             ++inc_index;
  92         }
  93     }
  94         for (i = 0; i < naxes[0]; ++i){
  95                 for (j = 0; j < naxes[1]; ++j) {
  96                         status = 0;
  97             /*Get the l and b for each pixel in the FITS file*/
  98                         fits_pix_to_world(i + 1, j + 1, wcs_in.xrefval, wcs_in.yrefval,
  99                                                           wcs_in.xrefpix, wcs_in.yrefpix, wcs_in.xinc,
 100                                                           wcs_in.yinc, wcs_in.rot, wcs_in.coordtype,
 101                                                           &ra0, &dec0, &status);
 102
 103                         if (status == 0){/*If the point is valid*/
 104                 /*Convert coordinates into cartesian coordinates*/
 105                                 x0 = cos(DEGRAD(ra0))*cos(DEGRAD(dec0));
 106                                 y0 = sin(DEGRAD(ra0))*cos(DEGRAD(dec0));
 107                                 z0 = sin(DEGRAD(dec0));
 108                 /*Run through all the data points find the angular distance between them*/
 109                                 for (index = 0; index < inc_index; ++index) {
 110                                         dist = bkgd.x[inc_list[index]]*x0 +
 111                            bkgd.y[inc_list[index]]*y0 +
 112                            bkgd.z[inc_list[index]]*z0;
 113                     /*If the distance is less than our limit, we bin weighting with the square
 114                      of the cos(angle) Note that because we want nearer points to be weighted
 115                      higher we actually use 1/cos(angle)*/
 116                                         if (dist > cang_limit){
 117                                                 data[i + j*naxes[0]] += bkgd.bkgd[inc_list[index]]*dist*dist;
 118                                                 var[i + j*naxes[0]] += dist*dist;
 119                                         }/*Endif*/
 120                                 }/*End check through all data values*/
 121                                 if (var[i + j*naxes[0]] > 0){
 122                     /*Divide by the weights*/
 123                                         data[i + j*naxes[0]] /= var[i + j*naxes[0]];
 124                                         var[i + j*naxes[0]] = sqrt(var[i + j*naxes[0]]);
 125                                 } else {/*If there is no valid data*/
 126                                         data[i + j*naxes[0]] = -1;
 127                                         var[i + j*naxes[0]] = 1e9;
 128                                 }
 129                                 if (data[i + j*naxes[0]] > max_data)
 130                     max_data = data[i + j*naxes[0]];
 131                         }/*End add data*/
 132                 }/*End j loop*/
 133         }/*End i loop*/
 134
 135     /*Fill in the empty areas using a cosecant law*/
 136         for (i = 0; i < naxes[0]; ++i){
 137                 for (j = 0; j < naxes[1]; ++j) {
 138                         if (var[i + j*naxes[0]] >= 1e9) {
 139                                 status = 0;
 140                                 fits_pix_to_world(i + 1, j + 1, wcs_in.xrefval, wcs_in.yrefval,
 141                                                                   wcs_in.xrefpix, wcs_in.yrefpix, wcs_in.xinc,
 142                                                                   wcs_in.yinc, wcs_in.rot, wcs_in.coordtype,
 143                                                                   &ra0, &dec0, &status);
 144                                 if (fabs(dec0) > 0) val = csc_law/sin(fabs(dec0)*RADEG);
 145                 if (val > max_data)
 146                     val = max_data;/*Cap the values*/
 147                                 data[i + j*naxes[0]] = max_data;
 148                         }
 149                 }
 150         }
 151
 152         return(EXIT_SUCCESS);
 153
 154 }