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 *--------------------------------------------------------------------------*/
67 #ifdef HAVE_CONFIG_H
68 #include <config.h>
69 #endif
70 #include <libnova/solar.h>
71 #include <libnova/transform.h>
72 #include <libnova/julian_day.h>
73 #include <libnova/utility.h>
74 #include "zodiacal_model.h"
76 #define PROGRAM "UVS_Zodiacal_Model"
78 /* ------------------------------------------------------------------------- */
80 void PRINT_BANNER()
81 {
82 printf("-----------------------------------------------------------\n");
83 printf("%s %s\n", PROGRAM, VERSION);
84 printf(" Calculates the Zodiacal light spectrum towards a given\n");
85 printf(" direction, on a given date and time.\n");
86 printf("-----------------------------------------------------------\n");
87 }/* End PRINT_BANNER */
89 /* ------------------------------------------------------------------------- */
91 void usage()
92 {
94 printf (" Usage: uvs_zodiacal_model \n\n");
95 printf (" If no parameter file exists, new one will be created\n");
96 printf (" which may be edited as required\n\n");
97 printf ("-----------------------------------------------------------\n");
98 printf ("Copyright (c) 2012, The CADS Team. [GPL v2.0 or later]\n");
99 printf ("Bug Reports : cads.iiap.res.in/dropbux \n");
100 printf ("Documentation : cads.iiap.res.in/software\n\n");
101 }/*End usage */
103 /* ------------------------------------------------------------------------- */
105 void get_sun_hecl(float t_ut, int day, int month, int year,
106 float ra, float dec, float *sun_hecl, float *sun_beta)
107 {
108 double jd = 0.0;
109 double f_temp = 0.0;
110 struct ln_lnlat_posn sun_ecl;
111 struct ln_date date;
112 struct ln_equ_posn eqcoords;
113 struct ln_lnlat_posn eclcoords;
115 /* calculate Julian day */
116 date.years = year;
117 date.months = month;
118 date.days = day;
119 date.hours = (int) t_ut;
120 date.minutes = 0;
121 date.seconds = 0;
122 jd = ln_get_julian_day(&date);
124 /* calculate sun's ecliptic coordinates */
125 ln_get_solar_ecl_coords(jd, &sun_ecl);
127 eqcoords.ra = (double) ra;
128 eqcoords.dec = (double) dec;
130 ln_get_ecl_from_equ(&eqcoords, jd, &eclcoords);
132 /* things are symmetric around 0 latitude here.. */
133 f_temp = fabs(eclcoords.lng - sun_ecl.lng);
135 /* wrap coordinates around 180 deg.*/
136 if (f_temp > 180) f_temp = 360 - f_temp;
138 *sun_hecl = (float) f_temp;
139 *sun_beta = (float) fabs(eclcoords.lat);
140 }/*End calculate solar coordinates*/
142 /* ------------------------------------------------------------------------- */
144 int main (int argc, char *argv[])
145 {
147 char param_file[MAX_TEXT] = "", out_file[MAX_TEXT] = "";
148 char zod_file[MAX_TEXT] = "", spec_file[MAX_TEXT] = "";
150 float *zod_dist, *zod_spec, *zod_wave;
151 float *table_ecl, *table_beta;
152 float ra = 0.0, dec = 0.0, hour = 0.0;
153 int day = 0, month = 0, year = 0;
154 float helio_beta = 0.0, helio_ecl = 0.0;
155 int i_arg = 0;
156 long hecl_index = 0, hbeta_index = 0, i = 0, wave_index = 0;
157 float zod_intensity = 0.0;
158 float zod_scale = 252; /*1 unit in the table in phot units at 5000 A*/
159 int ndat = 0;
160 FILE *fout = NULL;
162 /* --------------------------------------------------------------------- */
163 /* Before we begin... set the signals and print a banner */
165 set_signals(); /* set signals before starting */
166 PRINT_BANNER();
168 /* start with some default values. Use it
169 if nobody re-defined them in paramfile */
170 strcpy(spec_file, DEFAULT_SPEC_FILE);
171 strcpy(zod_file, DEFAULT_ZOD_FILE);
172 strcpy(out_file, DEFAULT_OUT_FILE);
174 /*Read parameters if present*/
175 /* Get all the input parameters from parameter file */
176 READ_PARAMS(&hour, &day, &month, &year, &ra, &dec, zod_file, spec_file,
177 out_file);
179 /* Tell the user about the input values */
180 printf("Info : Time & date : %04.2f Hrs %02d/%02d/%04d \n",
181 hour, day, month, year);
182 printf("Info : Sky coordinates : (RA = %3.4f, Decl. = %3.4f)\n",
183 ra, dec);
185 /*Zodiacal light related variables*/
186 zod_dist = (float *) malloc(sizeof(float) * ZOD_XSIZE * ZOD_YSIZE);
187 zod_spec = (float *) malloc(sizeof(float) * ZOD_SPEC);
188 zod_wave = (float *) malloc(sizeof(float) * ZOD_SPEC);
189 table_ecl = (float *) malloc(sizeof(float) * ZOD_XSIZE);
190 table_beta = (float *) malloc(sizeof(float) * ZOD_YSIZE);
192 /*Read zodiacal light data*/
193 ZOD_DIST_READ(zod_dist, table_ecl, table_beta, zod_file);
194 ZOD_SPECT_READ(zod_wave, zod_spec, &ndat, spec_file);
196 /* Get helio-ecliptic coordinates of Sun */
197 get_sun_hecl(hour, day, month, year, ra, dec, &helio_ecl, &helio_beta);
199 /* Print ecliptic coordinates*/
200 printf("Info : Helioecliptic coordinates are: %3.4f %3.4f\n",
201 helio_ecl, helio_beta);
203 /*Now lookup zodiacal intensity*/
204 hecl_index = 0;
205 while (table_ecl[hecl_index] < helio_ecl)
206 ++hecl_index;
207 hbeta_index = 0;
208 while (table_beta[hbeta_index] < fabs(helio_beta))
209 ++hbeta_index;
211 zod_intensity = zod_dist[hecl_index + hbeta_index*ZOD_XSIZE];
213 /*Scale zodiacal light spectrum*/
214 wave_index = 0;
215 while (zod_wave[wave_index] < 5000)
216 ++wave_index;
218 /*Scale factor for zodiacal light*/
219 zod_scale *= (zod_intensity / zod_spec[wave_index]);
220 for (i = 0; i < ndat; ++i)
221 zod_spec[i] *= zod_scale;
225 /*Print out zodiacal light*/
227 fout = fopen(out_file, "w");
228 if (fout == NULL)
229 {
230 fprintf(stderr, "ERROR : Unable to create output file, %s\n",
231 out_file);
232 exit(EXIT_FAILURE);
233 }
235 /* Write a header with some info on developers and all input values */
236 fprintf(fout, "#------------------------------------------------#\n");
237 fprintf(fout, "# Zodiacal light Spectrum generated by\n");
238 fprintf(fout, "# %s %s\n", PROGRAM, VERSION);
239 fprintf(fout, "# (based on Leinert's model) \n");
240 fprintf(fout, "# Copyright (c) by the CADS/UV Software Team\n");
241 fprintf(fout, "# URL: http://cads.iiap.res.in/software\n");
242 fprintf(fout, "#------------------------------------------------#\n");
243 fprintf(fout, "# user inputs\n");
244 fprintf(fout, "# Time : %f Hrs %02d/%02d/%04d \n",
245 hour, day, month, year);
246 fprintf(fout, "# Sky Coord : (RA = %3.4f, Decl. = %3.4f)\n", ra, dec);
247 fprintf(fout, "#\n");
248 fprintf(fout, "# Solar spectrum : %s\n", spec_file);
249 fprintf(fout, "# Zodiacal model : %s\n", zod_file);
250 fprintf(fout, "#------------------------------------------------#\n");
251 fprintf(fout, "# Column 1: Wavelength (Angstroms) \n");
252 fprintf(fout, "# Column 2: Flux density (photons/cm^2/s/sr/A) \n");
253 fprintf(fout, "#------------------------------------------------#\n");
255 for (i = 0; i < ndat; ++i)
256 fprintf(fout, "%f %f\n", zod_wave[i], zod_spec[i]);
257 close(fout);
259 printf("Info : Output written to %s\n", out_file);
260 return(EXIT_SUCCESS);
261 }