src/zodiacal_model.c

   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)
  30
  31  The output is in units of photons cm-2 s-1 sr-1 A-1
  32
  33  Date  : Sep. 8, 2007
  34  Author: Jayant Murthy
  35
  36   Modification history:
  37
  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.
  44
  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.
  48
  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).
  52
  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.
  59
  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"
  75
  76 #define PROGRAM "UVS_Zodiacal_Model"
  77
  78 /* ------------------------------------------------------------------------- */
  79
  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 */
  88
  89 /* ------------------------------------------------------------------------- */
  90
  91 void usage()
  92 {
  93
  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 */
 102
 103 /* ------------------------------------------------------------------------- */
 104
 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;
 114
 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);
 123
 124     /* calculate sun's ecliptic coordinates */
 125     ln_get_solar_ecl_coords(jd, &sun_ecl);
 126
 127     eqcoords.ra  = (double) ra;
 128     eqcoords.dec = (double) dec;
 129
 130     ln_get_ecl_from_equ(&eqcoords, jd, &eclcoords);
 131
 132     /* things are symmetric around 0 latitude here..  */
 133     f_temp    = fabs(eclcoords.lng - sun_ecl.lng);
 134
 135     /* wrap coordinates around 180 deg.*/
 136     if (f_temp > 180) f_temp = 360 - f_temp;
 137
 138     *sun_hecl = (float) f_temp;
 139     *sun_beta = (float) fabs(eclcoords.lat);
 140 }/*End calculate solar coordinates*/
 141
 142 /* ------------------------------------------------------------------------- */
 143
 144 int main (int argc, char *argv[])
 145 {
 146
 147     char param_file[MAX_TEXT] = "", out_file[MAX_TEXT] = "";
 148     char zod_file[MAX_TEXT] = "", spec_file[MAX_TEXT] = "";
 149
 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;
 161
 162     /* --------------------------------------------------------------------- */
 163     /*                 Before we begin... set the signals and print a banner */
 164
 165     set_signals();                            /* set signals before starting */
 166     PRINT_BANNER();
 167
 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);
 173
 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);
 178
 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);
 184
 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);
 191
 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);
 195
 196     /* Get helio-ecliptic coordinates of Sun */
 197     get_sun_hecl(hour, day, month, year, ra, dec, &helio_ecl, &helio_beta);
 198
 199     /* Print ecliptic coordinates*/
 200     printf("Info   : Helioecliptic coordinates are: %3.4f %3.4f\n",
 201            helio_ecl, helio_beta);
 202
 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;
 210
 211     zod_intensity = zod_dist[hecl_index + hbeta_index*ZOD_XSIZE];
 212
 213     /*Scale zodiacal light spectrum*/
 214     wave_index = 0;
 215     while (zod_wave[wave_index] < 5000)
 216         ++wave_index;
 217
 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;
 222
 223
 224
 225     /*Print out zodiacal light*/
 226
 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     }
 234
 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");
 254
 255     for (i = 0; i < ndat; ++i)
 256         fprintf(fout, "%f %f\n", zod_wave[i], zod_spec[i]);
 257     close(fout);
 258
 259     printf("Info   : Output written to %s\n", out_file);
 260     return(EXIT_SUCCESS);
 261 }