e1f8724bf59128ad8a936146b477a69c9b40b33b
1 /***************************************************************************
2 * This File is a part of the CADS/UV Zodiacal light model software *
3 * Copyright (C) 2012 by CADS/UV Software Team, *
4 * Indian Institute of Astrophysics *
5 * Bangalore 560034 *
6 * cads_AT_iiap.res.in *
7 * *
8 * This program is free software; you can redistribute it and/or modify *
9 * it under the terms of the GNU General Public License as published by *
10 * the Free Software Foundation; either version 2 of the License, or *
11 * (at your option) any later version. *
12 * *
13 * This program is distributed in the hope that it will be useful, *
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16 * GNU General Public License for more details. *
17 * *
18 * You should have received a copy of the GNU General Public License *
19 * along with this program; if not, write to the *
20 * Free Software Foundation, Inc., *
21 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
22 ***************************************************************************/
23 /*-------------------------------------------------------------------------*
24 Program : cads_zodiacal_model.c
25 This program will write a 2 dimensional file with the zodiacal light
26 spectrum. The spectrum is based on the solar spectrum from Colina et al
27 while the distribution is from Leinert et al. Documentation is in the
28 cads website at http://cads.iiap.res.in/software The Sun position is
29 calculated using libnova (http://libnova.sf.net)
31 The output is in units of photons cm-2 s-1 sr-1 A-1
33 Date : Sep. 8, 2007
34 Author: Jayant Murthy
36 Modification history:
38 Reks, 17th Sept 2007:
39 1. Improved formatting and overall beautification :-)
40 2. configure, make and added into the simulation project
41 3. Removed parse_par. Inputs are now read from a parameter file,
42 including the location of leinart model and zodiacal spectrum files.
43 4. output written to a file instead of terminal.
45 JM: Sep. 20, 2007:
46 1. helioecliptic longitude should be absolute value.
47 2. latitude should also be fabs because we are symmetric around the equator.
49 Reks: 25th Sept 2007: v1.0 ready for release!
50 1. Put back all the options for sending in input parameters - paramfile
51 or commandline or interactive (in the order of preference).
53 JM: Sep. 27, 2007
54 1. Changed input methods. Will read either from a parameter file or
55 from a command line but not interactively.
56 2. Corrected an error where I need to look only at the difference
57 between solar longitude and look direction.
58 3. Header of output file was wrong.
60 Reks: 23 June 2012
61 1. Rebranded to cads software team
62 2. All coordinate transformations and calculations using libnova
63 3. Removed parse_par => no command line options
64 4. Can handle comments in spectrum file and removed the 'rather strict'
65 ZOD_SPEC number of elements in spectrum :)
66 5. More gnu-autotools friendly. Path to input data files are taken care of
67 *--------------------------------------------------------------------------*/
68 #include <libnova/solar.h>
69 #include <libnova/transform.h>
70 #include <libnova/julian_day.h>
71 #include <libnova/utility.h>
72 #include "zodiacal_model.h"
74 #define PROGRAM "UVS_Zodiacal_Model"
76 /* ------------------------------------------------------------------------- */
78 void PRINT_BANNER()
79 {
80 printf("-----------------------------------------------------------\n");
81 printf("%s %s\n", PROGRAM, VERSION);
82 printf(" Calculates the Zodiacal light spectrum towards a given\n");
83 printf(" direction, on a given date and time.\n");
84 printf("-----------------------------------------------------------\n");
85 }/* End PRINT_BANNER */
87 /* ------------------------------------------------------------------------- */
89 void usage()
90 {
92 printf (" Usage: uvs_zodiacal_model \n\n");
93 printf (" If no parameter file exists, new one will be created\n");
94 printf (" which may be edited as required\n\n");
95 printf ("-----------------------------------------------------------\n");
96 printf ("Copyright (c) 2012 The CADS Team, IIAp. [GPL v3.0 or later]\n");
97 printf ("Report Bugs to: http://cads.iiap.res.in/bugzilla \n");
98 printf ("Documentation : http://cads.iiap.res.in/software\n\n");
99 }/*End usage */
101 /* ------------------------------------------------------------------------- */
103 void get_sun_hecl(float t_ut, int day, int month, int year,
104 float ra, float dec, float *sun_hecl, float *sun_beta)
105 {
106 double jd = 0.0;
107 double f_temp = 0.0;
108 struct ln_lnlat_posn sun_ecl;
109 struct ln_date date;
110 struct ln_equ_posn eqcoords;
111 struct ln_lnlat_posn eclcoords;
113 /* calculate Julian day */
114 date.years = year;
115 date.months = month;
116 date.days = day;
117 date.hours = (int) t_ut;
118 date.minutes = 0;
119 date.seconds = 0;
120 jd = ln_get_julian_day(&date);
122 /* calculate sun's ecliptic coordinates */
123 ln_get_solar_ecl_coords(jd, &sun_ecl);
125 eqcoords.ra = (double) ra;
126 eqcoords.dec = (double) dec;
128 ln_get_ecl_from_equ(&eqcoords, jd, &eclcoords);
130 /* things are symmetric around 0 latitude here.. */
131 f_temp = fabs(eclcoords.lng - sun_ecl.lng);
133 /* wrap coordinates around 180 deg.*/
134 if (f_temp > 180) f_temp = 360 - f_temp;
136 *sun_hecl = (float) f_temp;
137 *sun_beta = (float) fabs(eclcoords.lat);
138 }/*End calculate solar coordinates*/
140 /* ------------------------------------------------------------------------- */
142 int main (int argc, char *argv[])
143 {
145 char param_file[MAX_TEXT] = "", out_file[MAX_TEXT] = "";
146 char zod_file[MAX_TEXT] = "", spec_file[MAX_TEXT] = "";
148 float *zod_dist, *zod_spec, *zod_wave;
149 float *table_ecl, *table_beta;
150 float ra = 0.0, dec = 0.0, hour = 0.0;
151 int day = 0, month = 0, year = 0;
152 float helio_beta = 0.0, helio_ecl = 0.0;
153 int i_arg = 0;
154 long hecl_index = 0, hbeta_index = 0, i = 0, wave_index = 0;
155 float zod_intensity = 0.0;
156 float zod_scale = 252; /*1 unit in the table in phot units at 5000 A*/
157 int ndat = 0;
158 FILE *fout = NULL;
160 /* --------------------------------------------------------------------- */
161 /* Before we begin... set the signals and print a banner */
163 set_signals(); /* set signals before starting */
164 PRINT_BANNER();
166 /* start with some default values. Use it
167 if nobody re-defined them in paramfile */
168 strcpy(spec_file, DEFAULT_SPEC_FILE);
169 strcpy(zod_file, DEFAULT_ZOD_FILE);
170 strcpy(out_file, DEFAULT_OUT_FILE);
172 /*Read parameters if present*/
173 /* Get all the input parameters from parameter file */
174 READ_PARAMS(&hour, &day, &month, &year, &ra, &dec, zod_file, spec_file,
175 out_file);
177 /* Tell the user about the input values */
178 printf("Info : Time & date : %04.2f Hrs %02d/%02d/%04d \n",
179 hour, day, month, year);
180 printf("Info : Sky coordinates : (RA = %3.4f, Decl. = %3.4f)\n",
181 ra, dec);
183 /*Zodiacal light related variables*/
184 zod_dist = (float *) malloc(sizeof(float) * ZOD_XSIZE * ZOD_YSIZE);
185 zod_spec = (float *) malloc(sizeof(float) * ZOD_SPEC);
186 zod_wave = (float *) malloc(sizeof(float) * ZOD_SPEC);
187 table_ecl = (float *) malloc(sizeof(float) * ZOD_XSIZE);
188 table_beta = (float *) malloc(sizeof(float) * ZOD_YSIZE);
190 /*Read zodiacal light data*/
191 ZOD_DIST_READ(zod_dist, table_ecl, table_beta, zod_file);
192 ZOD_SPECT_READ(zod_wave, zod_spec, &ndat, spec_file);
194 /* Get helio-ecliptic coordinates of Sun */
195 get_sun_hecl(hour, day, month, year, ra, dec, &helio_ecl, &helio_beta);
197 /* Print ecliptic coordinates*/
198 printf("Info : Helioecliptic coordinates are: %3.4f %3.4f\n",
199 helio_ecl, helio_beta);
201 /*Now lookup zodiacal intensity*/
202 hecl_index = 0;
203 while (table_ecl[hecl_index] < helio_ecl)
204 ++hecl_index;
205 hbeta_index = 0;
206 while (table_beta[hbeta_index] < fabs(helio_beta))
207 ++hbeta_index;
209 zod_intensity = zod_dist[hecl_index + hbeta_index*ZOD_XSIZE];
211 /*Scale zodiacal light spectrum*/
212 wave_index = 0;
213 while (zod_wave[wave_index] < 5000)
214 ++wave_index;
216 /*Scale factor for zodiacal light*/
217 zod_scale *= (zod_intensity / zod_spec[wave_index]);
218 for (i = 0; i < ndat; ++i)
219 zod_spec[i] *= zod_scale;
223 /*Print out zodiacal light*/
225 fout = fopen(out_file, "w");
226 if (fout == NULL)
227 {
228 fprintf(stderr, "ERROR : Unable to create output file, %s\n",
229 out_file);
230 exit(EXIT_FAILURE);
231 }
233 /* Write a header with some info on developers and all input values */
234 fprintf(fout, "#------------------------------------------------#\n");
235 fprintf(fout, "# Zodiacal light Spectrum generated by\n");
236 fprintf(fout, "# %s %s\n", PROGRAM, VERSION);
237 fprintf(fout, "# (based on Leinert's model) \n");
238 fprintf(fout, "# Copyright (c) by the CADS/UV Software Team\n");
239 fprintf(fout, "# URL: http://cads.iiap.res.in/software\n");
240 fprintf(fout, "#------------------------------------------------#\n");
241 fprintf(fout, "# user inputs\n");
242 fprintf(fout, "# Time : %f Hrs %02d/%02d/%04d \n",
243 hour, day, month, year);
244 fprintf(fout, "# Sky Coord : (RA = %3.4f, Decl. = %3.4f)\n", ra, dec);
245 fprintf(fout, "#\n");
246 fprintf(fout, "# Solar spectrum : %s\n", spec_file);
247 fprintf(fout, "# Zodiacal model : %s\n", zod_file);
248 fprintf(fout, "#------------------------------------------------#\n");
249 fprintf(fout, "# Column 1: Wavelength (Angstroms) \n");
250 fprintf(fout, "# Column 2: Flux density (photons/cm^2/s/sr/A) \n");
251 fprintf(fout, "#------------------------------------------------#\n");
253 for (i = 0; i < ndat; ++i)
254 fprintf(fout, "%f %f\n", zod_wave[i], zod_spec[i]);
255 close(fout);
257 printf("Info : Output written to %s\n", out_file);
258 return(EXIT_SUCCESS);
259 }