1 /****************************************************************************
2 This program will calculate the count rate per pixel
3 in a given area of the sky.
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
10 **************************************************************************/
11 #include "sky_model.h"
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;
21 /*Location*/
22 double ra, dec, x, y;
24 /*Input*/
25 struct INP_PAR inp_par;
27 /*Output*/
28 struct WCS_HEAD wcs_in;
29 long naxes[2];
30 int ix, iy;
31 float *diffuse_data;
33 /*Zodiacal variables*/
34 float *zod_data, zod_scale, tot_zod;
35 struct ZOD_DATA zodi;
36 struct SPECTRA zod_spec;
38 /*Filters*/
39 struct SPECTRA filter_spec;
41 /*Stars*/
42 struct STARS *hipstars;
43 struct SPECTRA *stellar_spectra;
44 int istars;
45 float tot_stars;
47 /*Unclassified variables*/
48 int status;
49 int i, tst;
50 float *err;
53 /*====================== READS OUT HEADERLINES===============================*/
55 /*Begin Initialization*/
56 /*Read parameters if present*/
57 tst = READ_PARAMS(&inp_par, &wcs_in);
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));
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]);
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);
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);
80 /*Filters*/
81 filter_spec.wavelength = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
82 filter_spec.spectrum = (float *) malloc(sizeof(float) * MAX_SPEC_ELEM);
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);
94 /*Read the filter*/
95 filter_spec.nelements = SPECT_READ(filter_spec.wavelength,
96 filter_spec.spectrum, inp_par.filter_file);
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 }
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) {
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 }
151 }/*Read all the stars into memory*/
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);
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 */