src/zodiacal_programs.c

   1 /*
   2  *  zodiacal_programs.c
   3  *
   4  *  Created by Jayant Murthy on 21/04/11.
   5  *
   6  */
   7
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <math.h>
  11 #include <string.h>
  12 #include "sky_model.h"
  13 #include "fitsio.h"
  14 #include "solar.h"
  15
  16 /********************** SUN_RA_DEC *****************************/
  17 /*Using libnova*/
  18
  19 void SUN_RA_DEC(float t_ut, float day, float month,
  20                 float year, double *ra, double *dec)
  21 {
  22         double jd;
  23     struct ln_equ_posn sun_eq;
  24
  25         /*Calculate Julian Day*/
  26         if (month<=2) {
  27                 month=month+12;
  28                 year=year-1;
  29         }
  30         jd = (int) (365.25*year) + (int)(30.6001*(month+1)) - 15 + 1720996.5 + day + t_ut/24.0;
  31     ln_get_solar_equ_coords (jd, &sun_eq);
  32         *ra = sun_eq.ra;
  33         *dec= sun_eq.dec;
  34 }/*End calculate solar coordinates*/
  35 /*********************************************************************/
  36
  37 /*************** ZOD_DIST_READ ****************************/
  38 /* Reads zodiacal light from Leinert table ***************/
  39
  40 struct ZOD_DATA ZOD_DIST_READ(char *zod_file)
  41 {
  42         FILE *ZODIAC=NULL;
  43         int ix, iy;
  44         float f_tmp;
  45         char s_tmp[FLEN_FILENAME];
  46     struct ZOD_DATA zodi;
  47
  48 /*Check for existence of zodiacal light file*/
  49         if ((ZODIAC = fopen(zod_file, "r")) == NULL){
  50                 fprintf(stderr, "ERROR  : Unable to open Zodiacal light file %s \n", zod_file);
  51                 exit(EXIT_FAILURE);
  52         }/*Open zodiacal light file if it exists*/
  53
  54 /*Zodiacal light related variables*/
  55         zodi.zod_dist   = (float *) malloc(sizeof(float) * ZOD_XSIZE * ZOD_YSIZE);
  56         zodi.table_ecl  = (float *) malloc(sizeof(float) * ZOD_XSIZE);
  57         zodi.table_beta = (float *) malloc(sizeof(float) * ZOD_YSIZE);
  58
  59
  60         /*The Zodiacal light file has two comment lines*/
  61         fgets(s_tmp, FLEN_FILENAME, ZODIAC);
  62         fgets(s_tmp, FLEN_FILENAME, ZODIAC);
  63
  64         fscanf(ZODIAC, "%f", &f_tmp);         /*Dummy read of first element*/
  65         for (iy = 0; iy < ZOD_YSIZE; ++iy){ /*Now read the remaining latitudes*/
  66                 fscanf(ZODIAC, "%f", &f_tmp);
  67                 zodi.table_beta[iy] = f_tmp;
  68         }
  69
  70         for (ix = 0; ix < ZOD_XSIZE; ++ix){
  71                 fscanf(ZODIAC, "%f", &f_tmp);
  72                 zodi.table_ecl[ix] = f_tmp;/*Reading longitudes*/
  73                 for (iy = 0; iy < ZOD_YSIZE; ++iy){
  74                         fscanf(ZODIAC, "%f", &f_tmp);/*Read actual values*/
  75                         zodi.zod_dist[ix + iy*ZOD_XSIZE] = f_tmp;
  76                 }
  77         }/*Read over zodiacal file*/
  78
  79     return(zodi);
  80         fclose(ZODIAC);
  81 }/*End ZOD_DIST*/
  82 /*********************************************************************/
  83
  84 /*************** CALC_ZOD_SCALE_FACTOR *******************/
  85 /* Calculates the zodiacal light over the entire sky based on
  86    helioecliptic coordinates.
  87  */
  88
  89 /* CALC_ZOD_SCALE_FACTOR*/
  90
  91 int CALC_ZOD_FLUX(struct ZOD_DATA zodi,
  92                                                   float hour, float day, float month, float year,
  93                                                   long *naxes,
  94                                                   float *zod_data, struct WCS_HEAD wcs_in)
  95 {
  96         int hecl_index, hbeta_index, ix, iy;
  97         float zod_scale;
  98         float user_beta, helio_ecl;
  99         int status=0;
 100         double ra, dec, sun_ra, sun_dec, sun_ecl, sun_beta, coord[2];
 101
 102         /*Read zodiacal light data*/
 103         /*Get position of Sun*/
 104         SUN_RA_DEC(hour, day, month, year, &sun_ra, &sun_dec);
 105
 106         /*We need the helioecliptic longitude and the beta angle for our model*/
 107         /*Convert celestial to ecliptic coordinates for Sun*/
 108         coord[0] = sun_ra;
 109         coord[1] = sun_dec;
 110         EQ_to_EC(coord);
 111         sun_ecl  = coord[0];
 112         sun_beta = coord[1];
 113
 114 /*Step through the requested field. The WCS coordinates are in the WCS structure.
 115   The axis size is given by naxes*/
 116
 117         for (ix = 0; ix < naxes[0]; ++ix) {
 118                 for (iy = 0; iy < naxes[1]; ++iy) {
 119                         fits_pix_to_world(ix + 1, iy + 1, wcs_in.xrefval, wcs_in.yrefval,
 120                                                           wcs_in.xrefpix, wcs_in.yrefpix, wcs_in.xinc,
 121                                                           wcs_in.yinc, wcs_in.rot, wcs_in.coordtype,
 122                                                           &ra, &dec, &status);/*Returns galactic coordinates*/
 123                         if (status == 0) {
 124                                 /*Convert galactic coordinates into helioecliptic*/
 125                                 coord[0] = ra;  coord[1] = dec;
 126                                 GA_to_EC(coord);
 127                                 helio_ecl = fabs(coord[0] - sun_ecl);
 128                                 if (helio_ecl > 180) helio_ecl = 360 - helio_ecl;
 129                                 user_beta = fabs(coord[1]);
 130
 131                                 /*Now lookup zodiacal intensity*/
 132                                 hecl_index = 0;
 133                                 while (zodi.table_ecl[hecl_index] < helio_ecl)
 134                                         ++hecl_index;
 135                 if (hecl_index > 0) {
 136                     if ((zodi.table_ecl[hecl_index] - helio_ecl) >
 137                         (helio_ecl - zodi.table_ecl[hecl_index - 1]))
 138                         --hecl_index;
 139                 }
 140                                 hbeta_index = 0;
 141                                 while (zodi.table_beta[hbeta_index] < user_beta)
 142                                         ++hbeta_index;
 143                 if (hbeta_index > 0) {
 144                     if ((zodi.table_beta[hbeta_index] - user_beta) >
 145                         (user_beta - zodi.table_beta[hbeta_index - 1]))
 146                         --hbeta_index;
 147                 }
 148                                 zod_scale = zodi.zod_dist[hecl_index + hbeta_index*ZOD_XSIZE];
 149                                 zod_data[ix + iy*naxes[0]] = zod_scale;
 150                         }
 151                         status = 0;
 152                 }
 153         }
 154         return(EXIT_SUCCESS);
 155 }