src/sky_model.c

   1 /****************************************************************************
   2  This program will calculate the count rate per pixel
   3  in a given area of the sky.
   4
   5  Author: Jayant Murthy (jmurthy@yahoo.com)
   6  Program Documentation: https://docs.google.com/document/d/1BIT1I2IkJfBDlnYJI9a-E5dtHK-9ZWWv5k4p8vAcmb4/edit
   7  Date: Oct. 7, 2012
   8  Oct. 18, 2012 (JM): Initialized variables, added integrated counts, DC level
   9
  10  **************************************************************************/
  11 #include "sky_model.h"
  12
  13 int main(int argc, char *argv[], char *envp[])
  14 {
  15     /*Diffuse variables*/
  16     struct SPECTRA bkgd_spec;
  17     struct  BKGD_DATA bkgd;
  18     float bkgd_scale, solid_angle, tot_bkgd;
  19     float *data;
  20
  21     /*Location*/
  22     double ra, dec, x, y;
  23
  24     /*Input*/
  25     struct  INP_PAR  inp_par;
  26
  27     /*Output*/
  28     struct  WCS_HEAD wcs_in;
  29     long    naxes[2];
  30         int ix, iy;
  31     float *diffuse_data;
  32
  33     /*Zodiacal variables*/
  34         float   *zod_data, zod_scale, tot_zod;
  35     struct  ZOD_DATA zodi;
  36     struct SPECTRA zod_spec;
  37
  38     /*Filters*/
  39     struct SPECTRA filter_spec;
  40
  41     /*Stars*/
  42         struct  STARS   *hipstars;
  43     struct  SPECTRA *stellar_spectra;
  44     int istars;
  45     float tot_stars;
  46
  47     /*Unclassified variables*/
  48     int status;
  49         int     i, tst;
  50         float   *err;
  51
  52
  53 /*====================== READS OUT HEADERLINES===============================*/
  54
  55     /*Begin Initialization*/
  56     /*Read parameters if present*/
  57     tst = READ_PARAMS(&inp_par, &wcs_in);
  58
  59         naxes[0] = wcs_in.angsize/fabs(wcs_in.xinc) + 1;
  60         naxes[1] = wcs_in.angsize/wcs_in.yinc + 1;
  61     wcs_in.xrefpix = naxes[0]/2.;
  62     wcs_in.yrefpix = naxes[1]/2.;
  63     /*We have to define the solid angle for diffuse sources*/
  64     solid_angle = fabs(DEGRAD(wcs_in.xinc)*DEGRAD(wcs_in.yinc));
  65
  66         /*Allocate arrays*/
  67     data  = calloc(sizeof(float), naxes[0]*naxes[1]);
  68     err   = calloc(sizeof(float), naxes[0]*naxes[1]);
  69         diffuse_data = calloc(sizeof(float), naxes[0]*naxes[1]);
  70
  71     /*Zodiacal light*/
  72         zod_data   = calloc(sizeof(float), naxes[0]*naxes[1]);
  73     zod_spec.wavelength = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
  74         zod_spec.spectrum   = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
  75
  76     /*Diffuse background*/
  77     bkgd_spec.wavelength = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
  78     bkgd_spec.spectrum = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
  79
  80     /*Filters*/
  81     filter_spec.wavelength = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
  82     filter_spec.spectrum = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
  83
  84     /*Stars*/
  85     stellar_spectra = (struct SPECTRA *)
  86                     malloc(sizeof(struct STARS)*N_CASTELLI_MODELS);
  87     for (i = 0; i < N_CASTELLI_MODELS; ++i){
  88         stellar_spectra[i].spectrum = malloc(sizeof(float) * N_CASTELLI_SPECTRA);
  89         stellar_spectra[i].wavelength = malloc(sizeof(float) * N_CASTELLI_SPECTRA);
  90         stellar_spectra[i].nelements = N_CASTELLI_MODELS;
  91     }
  92     hipstars = (struct STARS *) malloc(sizeof(struct STARS)*NSTARS);
  93
  94     /*Read the filter*/
  95     filter_spec.nelements = SPECT_READ(filter_spec.wavelength,
  96                                        filter_spec.spectrum, inp_par.filter_file);
  97
  98     /*Read Galactic background if needed*/
  99         if (inp_par.gal_inc > 0){
 100         bkgd = READ_GAL_BKGD(inp_par.bkgd_model);
 101         bkgd_spec.nelements = SPECT_READ(bkgd_spec.wavelength,
 102                                          bkgd_spec.spectrum, inp_par.bkgd_spec);
 103         bkgd_scale = CALC_SCALE_FACTOR(filter_spec, bkgd_spec);
 104                 CALC_BKGD_FLUX(wcs_in, naxes, data, err, inp_par.ang_limit,
 105                        inp_par.csc_law, bkgd);
 106     }
 107         if (inp_par.zod_inc > 0){
 108         zodi = ZOD_DIST_READ(inp_par.zod_file);
 109         zod_spec.nelements = SPECT_READ(zod_spec.wavelength,
 110                                         zod_spec.spectrum, inp_par.zod_spec);
 111                 tst = CALC_ZOD_FLUX(zodi, inp_par.hour, inp_par.day, inp_par.month,
 112                               inp_par.year, naxes, zod_data, wcs_in);
 113         zod_scale = CALC_SCALE_FACTOR(filter_spec, zod_spec);
 114     }
 115         tot_bkgd  =0;
 116     tot_zod = 0;
 117         for (iy = 0; iy < naxes[1]; ++iy) {
 118                 for (ix = 0; ix < naxes[0]; ++ix) {
 119             diffuse_data[ix + iy*naxes[0]] =
 120             data[ix + iy*naxes[0]]*bkgd_scale*solid_angle +
 121                         zod_data[ix + iy*naxes[0]]*zod_scale*solid_angle +
 122             inp_par.dc_level;
 123             tot_bkgd += data[ix + iy*naxes[0]]*bkgd_scale*solid_angle;
 124             tot_zod += zod_data[ix + iy*naxes[0]]*zod_scale*solid_angle;
 125                 }
 126         }
 127
 128     /*Do we include stars?*/
 129     if (inp_par.star_inc > 0) {
 130         tst = READ_HIPPARCOS_CAT(inp_par, stellar_spectra, hipstars,
 131                                  filter_spec);
 132         tot_stars = 0;
 133         for (istars = 0; istars < NSTARS; ++istars) {
 134
 135             status = 0;
 136             ra = hipstars[istars].gl;
 137             dec = hipstars[istars].gb;
 138             status = 0;
 139             fits_world_to_pix(ra, dec, wcs_in.xrefval, wcs_in.yrefval,
 140                           wcs_in.xrefpix, wcs_in.yrefpix, wcs_in.xinc,
 141                           wcs_in.yinc, wcs_in.rot, wcs_in.coordtype,
 142                           &x, &y, &status);
 143             ix = (int) x - 1;
 144             iy = (int) y - 1;
 145             if ((status == 0) && (ix >= 0) && (iy >= 0) && (ix < naxes[0]) && (iy < naxes[1])){
 146                 diffuse_data[ix + iy*naxes[0]] += hipstars[istars].flux_at_earth;
 147                 tot_stars += hipstars[istars].flux_at_earth;
 148             }
 149         }
 150
 151     }/*Read all the stars into memory*/
 152
 153     /*Write FITS file*/
 154     remove(inp_par.out_file);
 155     FITS_WRITE_FILE(inp_par.out_file, wcs_in, naxes, diffuse_data, "g");
 156     printf("Bkgd: %f, Zodiacal: %f, Stars: %f\n", tot_bkgd, tot_zod, tot_stars);
 157
 158 /*Free variables*/
 159     free(data);
 160     free(err);
 161     free(diffuse_data);
 162     free(zod_data);
 163     free(zod_spec.wavelength);
 164     free(zod_spec.spectrum);
 165     free(bkgd_spec.wavelength);
 166     free(bkgd_spec.spectrum);
 167     free(filter_spec.wavelength);
 168     free(filter_spec.spectrum);
 169     for (i = 0; i < N_CASTELLI_MODELS; ++i){
 170         free(stellar_spectra[i].spectrum);
 171         free(stellar_spectra[i].wavelength);
 172     }
 173     free(stellar_spectra);
 174     free(hipstars);
 175 }   /* End of MAIN */