1 /*
2 gal_bkgd.c
3 sky_model
5 Created by Jayant Murthy on 07/10/12.
6 Copyright (c) 2012 Jayant Murthy. All rights reserved.
7 */
8 #include "sky_model.h"
10 struct BKGD_DATA READ_GAL_BKGD(char bkgd_file[FLEN_FILENAME]) {
12 FILE *file_ptr;
13 int data_index;
14 double val, eval, t1, t2;
15 struct BKGD_DATA gal_bkgd;
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 }
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);
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);
56 }
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 {
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;
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;
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));
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);
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*/
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 }
152 return(EXIT_SUCCESS);
154 }