1 /*
2 * zodiacal_programs.c
3 *
4 * Created by Jayant Murthy on 21/04/11.
5 *
6 */
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <string.h>
12 #include "sky_model.h"
13 #include "fitsio.h"
14 #include "solar.h"
16 /********************** SUN_RA_DEC *****************************/
17 /*Using libnova*/
19 void SUN_RA_DEC(float t_ut, float day, float month,
20 float year, double *ra, double *dec)
21 {
22 double jd;
23 struct ln_equ_posn sun_eq;
25 /*Calculate Julian Day*/
26 if (month<=2) {
27 month=month+12;
28 year=year-1;
29 }
30 jd = (int) (365.25*year) + (int)(30.6001*(month+1)) - 15 + 1720996.5 + day + t_ut/24.0;
31 ln_get_solar_equ_coords (jd, &sun_eq);
32 *ra = sun_eq.ra;
33 *dec= sun_eq.dec;
34 }/*End calculate solar coordinates*/
35 /*********************************************************************/
37 /*************** ZOD_DIST_READ ****************************/
38 /* Reads zodiacal light from Leinert table ***************/
40 struct ZOD_DATA ZOD_DIST_READ(char *zod_file)
41 {
42 FILE *ZODIAC=NULL;
43 int ix, iy;
44 float f_tmp;
45 char s_tmp[FLEN_FILENAME];
46 struct ZOD_DATA zodi;
48 /*Check for existence of zodiacal light file*/
49 if ((ZODIAC = fopen(zod_file, "r")) == NULL){
50 fprintf(stderr, "ERROR : Unable to open Zodiacal light file %s \n", zod_file);
51 exit(EXIT_FAILURE);
52 }/*Open zodiacal light file if it exists*/
54 /*Zodiacal light related variables*/
55 zodi.zod_dist = (float *) malloc(sizeof(float) * ZOD_XSIZE * ZOD_YSIZE);
56 zodi.table_ecl = (float *) malloc(sizeof(float) * ZOD_XSIZE);
57 zodi.table_beta = (float *) malloc(sizeof(float) * ZOD_YSIZE);
60 /*The Zodiacal light file has two comment lines*/
61 fgets(s_tmp, FLEN_FILENAME, ZODIAC);
62 fgets(s_tmp, FLEN_FILENAME, ZODIAC);
64 fscanf(ZODIAC, "%f", &f_tmp); /*Dummy read of first element*/
65 for (iy = 0; iy < ZOD_YSIZE; ++iy){ /*Now read the remaining latitudes*/
66 fscanf(ZODIAC, "%f", &f_tmp);
67 zodi.table_beta[iy] = f_tmp;
68 }
70 for (ix = 0; ix < ZOD_XSIZE; ++ix){
71 fscanf(ZODIAC, "%f", &f_tmp);
72 zodi.table_ecl[ix] = f_tmp;/*Reading longitudes*/
73 for (iy = 0; iy < ZOD_YSIZE; ++iy){
74 fscanf(ZODIAC, "%f", &f_tmp);/*Read actual values*/
75 zodi.zod_dist[ix + iy*ZOD_XSIZE] = f_tmp;
76 }
77 }/*Read over zodiacal file*/
79 return(zodi);
80 fclose(ZODIAC);
81 }/*End ZOD_DIST*/
82 /*********************************************************************/
84 /*************** CALC_ZOD_SCALE_FACTOR *******************/
85 /* Calculates the zodiacal light over the entire sky based on
86 helioecliptic coordinates.
87 */
89 /* CALC_ZOD_SCALE_FACTOR*/
91 int CALC_ZOD_FLUX(struct ZOD_DATA zodi,
92 float hour, float day, float month, float year,
93 long *naxes,
94 float *zod_data, struct WCS_HEAD wcs_in)
95 {
96 int hecl_index, hbeta_index, ix, iy;
97 float zod_scale;
98 float user_beta, helio_ecl;
99 int status=0;
100 double ra, dec, sun_ra, sun_dec, sun_ecl, sun_beta, coord[2];
102 /*Read zodiacal light data*/
103 /*Get position of Sun*/
104 SUN_RA_DEC(hour, day, month, year, &sun_ra, &sun_dec);
106 /*We need the helioecliptic longitude and the beta angle for our model*/
107 /*Convert celestial to ecliptic coordinates for Sun*/
108 coord[0] = sun_ra;
109 coord[1] = sun_dec;
110 EQ_to_EC(coord);
111 sun_ecl = coord[0];
112 sun_beta = coord[1];
114 /*Step through the requested field. The WCS coordinates are in the WCS structure.
115 The axis size is given by naxes*/
117 for (ix = 0; ix < naxes[0]; ++ix) {
118 for (iy = 0; iy < naxes[1]; ++iy) {
119 fits_pix_to_world(ix + 1, iy + 1, wcs_in.xrefval, wcs_in.yrefval,
120 wcs_in.xrefpix, wcs_in.yrefpix, wcs_in.xinc,
121 wcs_in.yinc, wcs_in.rot, wcs_in.coordtype,
122 &ra, &dec, &status);/*Returns galactic coordinates*/
123 if (status == 0) {
124 /*Convert galactic coordinates into helioecliptic*/
125 coord[0] = ra; coord[1] = dec;
126 GA_to_EC(coord);
127 helio_ecl = fabs(coord[0] - sun_ecl);
128 if (helio_ecl > 180) helio_ecl = 360 - helio_ecl;
129 user_beta = fabs(coord[1]);
131 /*Now lookup zodiacal intensity*/
132 hecl_index = 0;
133 while (zodi.table_ecl[hecl_index] < helio_ecl)
134 ++hecl_index;
135 if (hecl_index > 0) {
136 if ((zodi.table_ecl[hecl_index] - helio_ecl) >
137 (helio_ecl - zodi.table_ecl[hecl_index - 1]))
138 --hecl_index;
139 }
140 hbeta_index = 0;
141 while (zodi.table_beta[hbeta_index] < user_beta)
142 ++hbeta_index;
143 if (hbeta_index > 0) {
144 if ((zodi.table_beta[hbeta_index] - user_beta) >
145 (user_beta - zodi.table_beta[hbeta_index - 1]))
146 --hbeta_index;
147 }
148 zod_scale = zodi.zod_dist[hecl_index + hbeta_index*ZOD_XSIZE];
149 zod_data[ix + iy*naxes[0]] = zod_scale;
150 }
151 status = 0;
152 }
153 }
154 return(EXIT_SUCCESS);
155 }