src/read_write.c

   1 /*
   2  *  read_write.c
   3  *  diffuse
   4  *
   5  *  Created by Jayant Murthy on 21/04/11.
   6  *  Copyright 2011 __MyCompanyName__. All rights reserved.
   7  *
   8  */
   9 #include "sky_model.h"
  10
  11
  12 /*Function to create FITS header*/
  13 void FITS_WRITE_FILE(char file_out[FLEN_FILENAME], struct WCS_HEAD wcs_out,
  14                                         long naxes[2], float *data, char gal[3])
  15 {
  16     int status=0;
  17     char str_out[FLEN_KEYWORD];
  18     fitsfile *fptr;
  19
  20     if (fits_create_file(&fptr, file_out, &status)){/*Create FITS file*/
  21                 printf("Error in creating %s", file_out);
  22                 exit(EXIT_FAILURE);
  23     }
  24     fits_create_img(fptr, FLOAT_IMG, 2, naxes, &status);
  25     fits_write_key(fptr, TDOUBLE, "CRVAL1", &wcs_out.xrefval,
  26                                    "REF POINT VALUE IN DEGREES", &status);
  27     fits_write_key(fptr, TDOUBLE, "CRPIX1", &wcs_out.xrefpix,
  28                                    "REF POINT PIXEL LOCATION", &status);
  29     fits_write_key(fptr, TDOUBLE, "CDELT1", &wcs_out.xinc,
  30                                    "DEGREES PER PIXEL", &status);
  31     fits_write_key(fptr, TDOUBLE, "CROTA1", &wcs_out.rot,
  32                                    "ROTATION FROM ACTUAL AXIS", &status);
  33     if (strcmp(gal, "g") == 0) strcpy(str_out,"GLON");
  34     else strcpy(str_out,"RA--");
  35     strcat(str_out, wcs_out.coordtype);
  36     fits_write_key(fptr, TSTRING, "CTYPE1", &str_out,
  37                    "COORDINATE TYPE", &status);
  38     fits_write_key(fptr, TDOUBLE, "CRVAL2", &wcs_out.yrefval,
  39                    "REF POINT VALUE IN DEGREES", &status);
  40     fits_write_key(fptr, TDOUBLE, "CRPIX2", &wcs_out.yrefpix,
  41                    "REF POINT PIXEL LOCATION", &status);
  42     fits_write_key(fptr, TDOUBLE, "CDELT2", &wcs_out.yinc,
  43                    "DEGREES PER PIXEL", &status);
  44     fits_write_key(fptr, TDOUBLE, "CROTA2", &wcs_out.rot,
  45                    "ROTATION FROM ACTUAL AXIS", &status);
  46     if (strcmp(gal, "g") == 0) strcpy(str_out,"GLAT");
  47     else strcpy(str_out,"DEC-");
  48     strcat(str_out, wcs_out.coordtype);
  49     fits_write_key(fptr, TSTRING, "CTYPE2", &str_out,
  50                    "COORDINATE TYPE", &status);
  51     fits_write_2d_flt(fptr, 1, naxes[0], naxes[0], naxes[1], data, &status);
  52         fits_close_file(fptr, &status);
  53 }/*End of FITS_WRITE_FILE*/
  54 /*********************************************************************/
  55
  56 void GEN_DEFAULT_PARAMS() {
  57
  58         FILE *fparam = NULL;
  59
  60         fparam = fopen(PARAM_FILE, "w");
  61         if(fparam == NULL){
  62                 fprintf(stderr, "ERROR  : Unable to create %s\n", PARAM_FILE);
  63                 exit(EXIT_FAILURE);
  64         }
  65
  66         fprintf(fparam,
  67                         "#-----------------------------------------------------#\n"
  68                         "# Parameter FILE for sky model                         \n"
  69                         "#                                                      \n"
  70                         "# Purpose: Input parameters for running the Diffuse    \n"
  71                         "#          light modelling software will be read from  \n"
  72                         "#          this file.                                  \n"
  73                         );
  74
  75         fprintf(fparam,
  76                         "# NOTE:    Any line beginning with a hash (\"#\") or   \n"
  77                         "#          those enclosed within \"/* -- */ \" are     \n"
  78                         "#          comments and will be ignored by the software\n"
  79                         "#                                                      \n"
  80                         "# format of this file is: KEY = VALUE                  \n"
  81                         "#                                                      \n"
  82                         "#-----------------------------------------------------#\n\n"
  83                         );
  84
  85         fprintf(fparam,
  86                         "# TIME & DATE Information:\n"
  87                         "# -----------------------\n"
  88                         "# 1. Time of the day (UT, hours. range: 0.0 - 24.0)\n"
  89                         "# 2. Day of the month (unitless, range: 0 - 31)\n"
  90                         "# 3. Month of the year (unitless, range: 1 - 12) \n"
  91                         "# 4. Year (unitless, range: 1900 - 2100 for max. accuracy) \n\n"
  92                         "TIME_UT      = 12.00\n"
  93                         "DAY_OF_MONTH = 1.0\n"
  94                         "MONTH        = 8\n"
  95                         "YEAR         = 2007\n\n"
  96                         );
  97
  98         /*The files */
  99         fprintf(fparam,
 100                 "# List of files:\n"
 101                 "# -----------------------\n"
 102                );
 103         fprintf(fparam,
 104                 "OUTPUT_FILE      = diffuse_output.fits\n"
 105                );
 106         fprintf(fparam,
 107                 "ZODIAC_MODEL     = %szodiacal_data/leinert_dist.txt\n",
 108                 CADSDATA);
 109         fprintf(fparam,
 110                 "ZODIAC_SPECTRUM  = %szodiacal_data/zodiacal_spec.txt\n",
 111                 CADSDATA);
 112         fprintf(fparam,
 113                 "BKGD_MODEL       = %sdiffuse_data/nuv_allsky.txt\n",
 114                 CADSDATA);
 115         fprintf(fparam,
 116                 "BKGD_SPECTRUM    = %sdiffuse_data/bkgd_spectrum.txt\n",
 117                 CADSDATA);
 118         fprintf(fparam,
 119                 "STELLAR_SPECTRUM = %sstellar_data/castelli\n",
 120                 CADSDATA);
 121         fprintf(fparam,
 122                 "CROSS_SEC        = %sstellar_data/cross_sec.txt\n",
 123                 CADSDATA);
 124         fprintf(fparam,
 125                 "STAR_CAT         = %sstellar_data/hip_main.txt\n",
 126                 CADSDATA);
 127         fprintf(fparam,
 128                 "FILTER_SPECTRUM  = %ssample_filters/galex_FUV.txt\n\n"
 129                );
 130
 131         fprintf(fparam,
 132                         "# WCS Parameters for output FITS file\n"
 133             "#GAL_COORD = 1 => output is in galactic coordinates, othrwise equatorial"
 134                         "# ------------------------\n"
 135                         "XREFVAL = 0\n"
 136                         "YREFVAL = 0\n"
 137                         "XINC = -0.5\n"
 138                         "YINC = 0.5\n"
 139             "ANGSIZE = 360\n"
 140                         "COORDTYPE = -AIT\n"
 141             "GAL_COORD = 1\n\n"
 142                         );
 143
 144         fprintf(fparam,
 145                         "# Background parameters\n"
 146                         "# ---------------------\n"
 147                         "ANG_LIMIT = 3\n"
 148                         "CSC_LAW = 450.0\n\n"
 149             "DC_LEVEL = 0\n\n"
 150                         );
 151
 152         fprintf(fparam,
 153                         "# Which elements to include\n"
 154                         "# -----------------------\n"
 155                         "ZOD_INC = 1\n"
 156                         "GAL_INC = 1\n"
 157                         "STAR_INC = 1\n\n"
 158                         );
 159
 160         printf("Default Parameter file has been created. Please check and modify, if necessary.\n");
 161         fclose(fparam);
 162         exit(EXIT_FAILURE);
 163
 164 }/* END GEN_DEFAULT_PARAMS */
 165 /*********************************************************************/
 166
 167 int READ_PARAMS(struct INP_PAR *inp_par, struct WCS_HEAD *wcs_in)
 168 /*
 169                 float *hour, float *day, float *month, float *year,
 170                                 float *wave_ref, float *ang_limit, float *csc_law,
 171                                 int *zod_inc, int *gal_inc,
 172                 char *zod_file, char *spec_file, char *bkgd_file, char *out_file,
 173                                 struct WCS_HEAD *wcs_in) {
 174 */
 175                 {
 176         FILE    *fp;
 177     char    line[FLEN_FILENAME], name[FLEN_FILENAME];
 178     char    equal[FLEN_FILENAME], data[FLEN_FILENAME];
 179
 180         if((fp = fopen(PARAM_FILE,"r")) == NULL){
 181                 GEN_DEFAULT_PARAMS();
 182                 if((fp = fopen(PARAM_FILE,"r")) == NULL)
 183                         exit(EXIT_FAILURE);
 184         }
 185         /*Inititalize parameters*/
 186     inp_par->hour = 0;
 187     inp_par->day = 0;
 188                     inp_par->month = 0;
 189                     inp_par->csc_law=0;
 190                     inp_par->ang_limit = 5;
 191                     inp_par->dc_level= 0;
 192                     inp_par->gal_inc = 1;
 193                     inp_par->star_inc = 1;
 194                     inp_par->zod_inc = 1;
 195
 196         /* Read input parameters */
 197         while (!feof(fp)) {
 198                 fgets(line, FLEN_FILENAME, fp);
 199                 if ((line[0] != '#') && (line[0] != '\n')) {
 200                         sscanf(line, "%s %s %s ", name, equal, data);
 201                         if(!strcmp("TIME_UT", name))        inp_par->hour = atof(data);
 202                         if(!strcmp("DAY_OF_MONTH", name))   inp_par->day = atof(data);
 203                         if(!strcmp("MONTH", name))          inp_par->month = atof(data);
 204                         if(!strcmp("YEAR", name))           inp_par->year = atof(data);
 205                         if(!strcmp("OUTPUT_FILE", name))    strcpy(inp_par->out_file, data);
 206                         if(!strcmp("ZODIAC_MODEL", name))      strcpy(inp_par->zod_file, data);
 207                         if(!strcmp("ZODIAC_SPECTRUM", name))   strcpy(inp_par->zod_spec, data);
 208                         if(!strcmp("FILTER_SPECTRUM", name))   strcpy(inp_par->filter_file, data);
 209             if(!strcmp("BKGD_SPECTRUM", name))   strcpy(inp_par->bkgd_spec, data);
 210                         if(!strcmp("BKGD_MODEL", name))         strcpy(inp_par->bkgd_model, data);
 211             if(!strcmp("CROSS_SEC", name))              strcpy(inp_par->sigma_file, data);
 212             if(!strcmp("STELLAR_SPECTRUM", name))
 213                                         strcpy(inp_par->castelli_file, data);
 214             if(!strcmp("STAR_CAT", name))
 215                 strcpy(inp_par->hipparcos_file, data);
 216                         if(!strcmp("ANG_LIMIT", name))          inp_par->ang_limit = atof(data);
 217                         if(!strcmp("CSC_LAW", name))            inp_par->csc_law = atof(data);
 218             if(!strcmp("DC_LEVEL", name))               inp_par->dc_level = atof(data);
 219             if (!strcmp("ZOD_INC", name))       inp_par->zod_inc = atoi(data);
 220             if (!strcmp("GAL_INC", name))       inp_par->gal_inc = atoi(data);
 221             if (!strcmp("STAR_INC", name))      inp_par->star_inc = atoi(data);
 222                         if(!strcmp("XREFVAL", name))        wcs_in->xrefval = atof(data);
 223                         if(!strcmp("YREFVAL", name))            wcs_in->yrefval = atof(data);
 224                         if(!strcmp("XINC", name))                       wcs_in->xinc = atof(data);
 225                         if(!strcmp("YINC", name))           wcs_in->yinc = atof(data);
 226                         if(!strcmp("ANGSIZE", name))        wcs_in->angsize = atof(data);
 227             if(!strcmp("COORDTYPE", name))      strcpy(wcs_in->coordtype, data);
 228
 229                 }
 230         }
 231     wcs_in->rot = 0;
 232         strcpy(line, "End of File");
 233         /*This is to take care of the end of the file*/
 234         fclose(fp);
 235         return(EXIT_SUCCESS);
 236 }/* END READ_PARAMS */
 237
 238 struct INP_PAR READ_INPUT_PARAMETERS(int argc, char *argv[])
 239 {
 240     int i_arg = 2;
 241     FILE *para_file=NULL;
 242     struct INP_PAR inp_par;
 243
 244     if (PARSE_PAR_CHAR(argc, argv, inp_par.par_file, &i_arg)){
 245         printf("Parameter file name is required on command line\n.");
 246         exit(EXIT_FAILURE);
 247     }/*File which contains names of other files*/
 248
 249     if((para_file = fopen(inp_par.par_file,"r")) == NULL){
 250         printf("Could not open file list:\n");
 251         exit(EXIT_FAILURE);
 252     }
 253     strcpy(inp_par.out_dir, "");
 254     /*Read list of files*/
 255     fscanf(para_file, "%s",inp_par.castelli_file); /*Location of Castelli data*/
 256     fscanf(para_file, "%s",inp_par.sigma_file); /* Optical constants */
 257     fscanf(para_file, "%s",inp_par.hipparcos_file); /*Location of Hipparcos data*/
 258     fscanf(para_file, "%s",inp_par.location_file); /* Locations for calculations */
 259     fscanf(para_file, "%s", inp_par.out_dir);   /*Writes all output files in this directory*/
 260     fclose(para_file);
 261
 262     /*Output file name is derived from parameters*/
 263     strcpy(inp_par.out_file,"total_scattered.dat");
 264
 265     return (inp_par);
 266 }
 267
 268 /**************** SPECT_READ**************************/
 269 /*Reads the spectral (2 column) files*/
 270 int SPECT_READ(float *wave, float *spec, char *spec_file)
 271 {
 272         FILE *spec_file_ptr = NULL;
 273         int ispec;
 274         float f_tmp1, f_tmp2;
 275         char s_tmp[FLEN_FILENAME];
 276
 277         if ((spec_file_ptr = fopen(spec_file, "r")) == NULL){
 278                 fprintf(stderr, "ERROR  : Unable to open %s \n", spec_file);
 279                 exit(EXIT_FAILURE);
 280         }/*Open spectrum if it is exists*/
 281
 282     ispec =0;
 283     while (!feof(spec_file_ptr)) {
 284         if (ispec >  MAX_SPEC_ELEM){
 285             printf("%s %s\n", "Max. size exceeded for file:",spec_file);
 286             exit(EXIT_FAILURE);
 287         }
 288         /*Read line*/
 289         fgets(s_tmp, FLEN_FILENAME, spec_file_ptr);
 290                 if (sscanf(s_tmp, "%f %f", &f_tmp1, &f_tmp2) == 2){
 291             *(wave + ispec) = f_tmp1;
 292             *(spec + ispec) = f_tmp2;
 293             ++ispec;
 294         }
 295     }
 296         fclose(spec_file_ptr);
 297     return (ispec);
 298 }/*END SPECT_READ*/
 299 /*********************************************************************/
 300
 301