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  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"
  73
  74 #define PROGRAM "UVS_Zodiacal_Model"
  75
  76 /* ------------------------------------------------------------------------- */
  77
  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 */
  86
  87 /* ------------------------------------------------------------------------- */
  88
  89 void usage()
  90 {
  91
  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 */
 100
 101 /* ------------------------------------------------------------------------- */
 102
 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;
 112
 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);
 121
 122     /* calculate sun's ecliptic coordinates */
 123     ln_get_solar_ecl_coords(jd, &sun_ecl);
 124
 125     eqcoords.ra  = (double) ra;
 126     eqcoords.dec = (double) dec;
 127
 128     ln_get_ecl_from_equ(&eqcoords, jd, &eclcoords);
 129
 130     /* things are symmetric around 0 latitude here..  */
 131     f_temp    = fabs(eclcoords.lng - sun_ecl.lng);
 132
 133     /* wrap coordinates around 180 deg.*/
 134     if (f_temp > 180) f_temp = 360 - f_temp;
 135
 136     *sun_hecl = (float) f_temp;
 137     *sun_beta = (float) fabs(eclcoords.lat);
 138 }/*End calculate solar coordinates*/
 139
 140 /* ------------------------------------------------------------------------- */
 141
 142 int main (int argc, char *argv[])
 143 {
 144
 145     char param_file[MAX_TEXT] = "", out_file[MAX_TEXT] = "";
 146     char zod_file[MAX_TEXT] = "", spec_file[MAX_TEXT] = "";
 147
 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;
 159
 160     /* --------------------------------------------------------------------- */
 161     /*                 Before we begin... set the signals and print a banner */
 162
 163     set_signals();                            /* set signals before starting */
 164     PRINT_BANNER();
 165
 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);
 171
 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);
 176
 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);
 182
 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);
 189
 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);
 193
 194     /* Get helio-ecliptic coordinates of Sun */
 195     get_sun_hecl(hour, day, month, year, ra, dec, &helio_ecl, &helio_beta);
 196
 197     /* Print ecliptic coordinates*/
 198     printf("Info   : Helioecliptic coordinates are: %3.4f %3.4f\n",
 199            helio_ecl, helio_beta);
 200
 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;
 208
 209     zod_intensity = zod_dist[hecl_index + hbeta_index*ZOD_XSIZE];
 210
 211     /*Scale zodiacal light spectrum*/
 212     wave_index = 0;
 213     while (zod_wave[wave_index] < 5000)
 214         ++wave_index;
 215
 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;
 220
 221
 222
 223     /*Print out zodiacal light*/
 224
 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     }
 232
 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");
 252
 253     for (i = 0; i < ndat; ++i)
 254         fprintf(fout, "%f %f\n", zod_wave[i], zod_spec[i]);
 255     close(fout);
 256
 257     printf("Info   : Output written to %s\n", out_file);
 258     return(EXIT_SUCCESS);
 259 }