src/stars.c

   1 /*
   2   stars.c
   3   sky_model
   4
   5   Created by Jayant Murthy on 07/10/12.
   6   Copyright (c) 2012 Jayant Murthy. All rights reserved.
   7 */
   8
   9 #include "sky_model.h"
  10 /*Function to read single line from current position in stars file*/
  11
  12 int HIP_READ_LINE(FILE *star_file, struct STARS *line)
  13 {
  14     int i, tst;
  15     char hip_line[HIP_MAIN_REC_LEN], *token, *read_res, s_tmp[13], s_tmp1[13];
  16     double  gl, gb;
  17
  18     read_res = fgets(hip_line, HIP_MAIN_REC_LEN, star_file);/*Read first bracket*/
  19     if (read_res == NULL)
  20         return(EXIT_FAILURE);
  21     token = strtok(hip_line, HIP_DELIM); /*Set up tokens*/
  22     for (i = 1; i < 77; ++i) {
  23         token = strtok(NULL, HIP_DELIM); /*read token*/
  24         switch (i) {
  25             case 1:  line->HIP_NO   = atoi(token);
  26                 break;
  27             case 8:  line->ra       = atof(token);
  28                 break;
  29             case 9:  line->dec      = atof(token);
  30                 break;
  31             case 11: if (atof(token) > 0)
  32                 line->distance = 1000./atof(token);
  33             else line->distance = 1.e6;
  34                 break;
  35             case 32: line->B_mag    = atof(token);
  36                 break;
  37             case 34: line->V_mag    = atof(token);
  38                 break;
  39             case 71: line->HD_NO    = atoi(token);
  40                 break;
  41             case 76: strcpy(s_tmp, token);
  42                 sscanf(s_tmp, "%s", s_tmp1);
  43                 strcpy(line->sp_type, s_tmp1);
  44                 break;
  45         }/*End switch*/
  46         tst = CONV_EQ_TO_GAL(line->ra, line->dec, &gl, &gb);
  47         line->gl = gl;
  48         line->gb = gb;
  49
  50         line->x = cos(DEGRAD(line->gl))*cos(DEGRAD(line->gb))*line->distance;
  51         line->y = sin(DEGRAD(line->gl))*cos(DEGRAD(line->gb))*line->distance;
  52         line->z = sin(DEGRAD(line->gb))*line->distance;
  53     }/*End for*/
  54     return(0);
  55 }/*End READ_LINE*/
  56 /*********************************************************************/
  57 int READ_CASTELLI_SPECTRA(char *spec_dir, struct SPECTRA
  58                           *stellar_spectra)
  59 {
  60     int status = 0, anynull;
  61     fitsfile *fptr;
  62     char filename[MAX_FILE_LENGTH];
  63     int i;
  64     int temper[N_CASTELLI_MODELS], gindex[N_CASTELLI_MODELS];
  65     float nullval;
  66     char stemper[6];
  67
  68     /*Lookup table for temperatures and g*/
  69     for (i = 0; i <= 37; ++i)
  70         temper[i] = 50000 - i*1000;
  71     for (i = 38; i < 76; ++i)
  72         temper[i] = 13000 - (i - 37)*250;
  73     for (i = 0; i <= 4; ++i)
  74         gindex[i] = 12;
  75     for (i = 5; i <= 10; ++i)
  76         gindex[i] = 11;
  77     for (i = 11; i <= 63; ++i)
  78         gindex[i] = 10;
  79     for (i = 64; i <= 75; ++i)
  80         gindex[i] = 11;
  81
  82     /*Read stellar spectra*/
  83     for (i = 0; i < N_CASTELLI_MODELS; ++i) {
  84         strcpy(filename, spec_dir);
  85         strcat(filename, "/ckp00_");
  86         sprintf(stemper, "%i", temper[i]);
  87         strcat(filename, stemper);
  88         strcat(filename, ".fits");
  89         sprintf(stellar_spectra[i].filename, "%i", temper[i]);
  90         fits_open_file(&fptr, filename, READONLY, &status);
  91         fits_movabs_hdu(fptr, 2, NULL, &status);
  92         fits_read_col(fptr, TFLOAT, gindex[i], 1, 1, N_CASTELLI_SPECTRA,
  93                       &nullval, stellar_spectra[i].spectrum, &anynull,
  94                       &status);
  95         fits_read_col(fptr, TFLOAT, 1, 1, 1, N_CASTELLI_SPECTRA,
  96                       &nullval, stellar_spectra[i].wavelength, &anynull,
  97                       &status);
  98
  99         fits_close_file(fptr, &status);
 100     }
 101
 102     return(EXIT_SUCCESS);
 103 }
 104
 105 int GET_STAR_TEMP(struct STARS *hipline)
 106 {
 107     char sptype[10];
 108     char *str=sptype;
 109
 110     if (strncmp(hipline->sp_type, "sd", 2) == 0) {
 111         strcpy(sptype, hipline->sp_type);
 112         str = str + 2;
 113         strcpy(hipline->sp_type, str);
 114     }
 115     if (strncmp(hipline->sp_type, "O3", 2) == 0) {
 116         hipline->temperature = 5;
 117     }
 118     else if (strncmp(hipline->sp_type, "O4", 2) == 0) {
 119         hipline->temperature = 7;
 120     }
 121     else if (strncmp(hipline->sp_type, "O5", 2) == 0) {
 122         hipline->temperature = 10;
 123     }
 124     else if (strncmp(hipline->sp_type, "O6", 2) == 0) {
 125         hipline->temperature = 11;
 126     }
 127     else if (strncmp(hipline->sp_type, "O7", 2) == 0) {
 128         hipline->temperature = 13;
 129     }
 130     else if (strncmp(hipline->sp_type, "O8", 2) == 0) {
 131         hipline->temperature = 15;
 132     }
 133     else if (strncmp(hipline->sp_type, "O9", 2) == 0) {
 134         hipline->temperature = 18;
 135     }
 136     else if (strncmp(hipline->sp_type, "B0", 2) == 0) {
 137         hipline->temperature = 20;
 138     }
 139     else if (strncmp(hipline->sp_type, "B1", 2) == 0) {
 140         hipline->temperature = 25;
 141     }
 142     else if (strncmp(hipline->sp_type, "B2", 2) == 0) {
 143         hipline->temperature = 28;
 144     }
 145     else if (strncmp(hipline->sp_type, "B3", 2) == 0) {
 146         hipline->temperature = 31;
 147     }
 148     else if (strncmp(hipline->sp_type, "B4", 2) == 0) {
 149         hipline->temperature = 33;
 150     }
 151     else if (strncmp(hipline->sp_type, "B5", 2) == 0) {
 152         hipline->temperature = 35;
 153     }
 154     else if (strncmp(hipline->sp_type, "B6", 2) == 0) {
 155         hipline->temperature = 36;
 156     }
 157     else if (strncmp(hipline->sp_type, "B7", 2) == 0) {
 158         hipline->temperature = 37;
 159     }
 160     else if (strncmp(hipline->sp_type, "B8", 2) == 0) {
 161         hipline->temperature = 41;
 162     }
 163     else if (strncmp(hipline->sp_type, "B9", 2) == 0) {
 164         hipline->temperature = 49;
 165     }
 166     else if (strncmp(hipline->sp_type, "A0", 2) == 0) {
 167         hipline->temperature = 51;
 168     }
 169     else if (strncmp(hipline->sp_type, "A1", 2) == 0) {
 170         hipline->temperature = 52;
 171     }
 172     else if (strncmp(hipline->sp_type, "A2", 2) == 0) {
 173         hipline->temperature = 53;
 174     }
 175     else if (strncmp(hipline->sp_type, "A3", 2) == 0) {
 176         hipline->temperature = 56;
 177     }
 178     else if (strncmp(hipline->sp_type, "A4", 2) == 0) {
 179         hipline->temperature = 56;
 180     }
 181     else if (strncmp(hipline->sp_type, "A5", 2) == 0) {
 182         hipline->temperature = 56;
 183     }
 184     else if (strncmp(hipline->sp_type, "A6", 2) == 0) {
 185         hipline->temperature = 57;
 186     }
 187     else if (strncmp(hipline->sp_type, "A7", 2) == 0) {
 188         hipline->temperature = 58;
 189     }
 190     else if (strncmp(hipline->sp_type, "A8", 2) == 0) {
 191         hipline->temperature = 59;
 192     }
 193     else if (strncmp(hipline->sp_type, "A9", 2) == 0) {
 194         hipline->temperature = 60;
 195     }
 196     else if (strncmp(hipline->sp_type, "F0", 2) == 0) {
 197         hipline->temperature = 60;
 198     }
 199     else if (strncmp(hipline->sp_type, "F1", 2) == 0) {
 200         hipline->temperature = 63;
 201     }
 202     else if (strncmp(hipline->sp_type, "F2", 2) == 0) {
 203         hipline->temperature = 61;
 204     }
 205     else if (strncmp(hipline->sp_type, "F3", 2) == 0) {
 206         hipline->temperature = 62;
 207     }
 208     else if (strncmp(hipline->sp_type, "F4", 2) == 0) {
 209         hipline->temperature = 62;
 210     }
 211     else if (strncmp(hipline->sp_type, "F5", 2) == 0) {
 212         hipline->temperature = 63;
 213     }
 214     else if (strncmp(hipline->sp_type, "F6", 2) == 0) {
 215         hipline->temperature = 63;
 216     }
 217     else if (strncmp(hipline->sp_type, "F7", 2) == 0) {
 218         hipline->temperature = 63;
 219     }
 220     else if (strncmp(hipline->sp_type, "F8", 2) == 0) {
 221         hipline->temperature = 64;
 222     }
 223     else if (strncmp(hipline->sp_type, "F9", 2) == 0) {
 224         hipline->temperature = 65;
 225     }
 226     else if (strncmp(hipline->sp_type, "G0", 2) == 0) {
 227         hipline->temperature = 65;
 228     }
 229     else if (strncmp(hipline->sp_type, "G1", 2) == 0) {
 230         hipline->temperature = 66;
 231     }
 232     else if (strncmp(hipline->sp_type, "G2", 2) == 0) {
 233         hipline->temperature = 66;
 234     }
 235     else if (strncmp(hipline->sp_type, "G3", 2) == 0) {
 236         hipline->temperature = 66;
 237     }
 238     else if (strncmp(hipline->sp_type, "G4", 2) == 0) {
 239         hipline->temperature = 66;
 240     }
 241     else if (strncmp(hipline->sp_type, "G5", 2) == 0) {
 242         hipline->temperature = 66;
 243     }
 244     else if (strncmp(hipline->sp_type, "G6", 2) == 0) {
 245         hipline->temperature = 66;
 246     }
 247     else if (strncmp(hipline->sp_type, "G7", 2) == 0) {
 248         hipline->temperature = 67;
 249     }
 250     else if (strncmp(hipline->sp_type, "G8", 2) == 0) {
 251         hipline->temperature = 67;
 252     }
 253     else if (strncmp(hipline->sp_type, "G9", 2) == 0) {
 254         hipline->temperature = 67;
 255     }
 256     else if (strncmp(hipline->sp_type, "K0", 2) == 0) {
 257         hipline->temperature = 68;
 258     }
 259     else if (strncmp(hipline->sp_type, "K1", 2) == 0) {
 260         hipline->temperature = 69;
 261     }
 262     else if (strncmp(hipline->sp_type, "K2", 2) == 0) {
 263         hipline->temperature = 70;
 264     }
 265     else if (strncmp(hipline->sp_type, "K3", 2) == 0) {
 266         hipline->temperature = 70;
 267     }
 268     else if (strncmp(hipline->sp_type, "K4", 2) == 0) {
 269         hipline->temperature = 71;
 270     }
 271     else if (strncmp(hipline->sp_type, "K5", 2) == 0) {
 272         hipline->temperature = 72;
 273     }
 274     else if (strncmp(hipline->sp_type, "K6", 2) == 0) {
 275         hipline->temperature = 72;
 276     }
 277     else if (strncmp(hipline->sp_type, "K7", 2) == 0) {
 278         hipline->temperature = 73;
 279     }
 280     else if (strncmp(hipline->sp_type, "K8", 2) == 0) {
 281         hipline->temperature = 73;
 282     }
 283     else if (strncmp(hipline->sp_type, "K9", 2) == 0) {
 284         hipline->temperature = 73;
 285     }
 286     else if (strncmp(hipline->sp_type, "M0", 2) == 0) {
 287         hipline->temperature = 74;
 288     }
 289     else if (strncmp(hipline->sp_type, "M1", 2) == 0) {
 290         hipline->temperature = 74;
 291     }
 292     else if (strncmp(hipline->sp_type, "M2", 2) == 0) {
 293         hipline->temperature = 75;
 294     }
 295     else if (strncmp(hipline->sp_type, "M3", 2) == 0) {
 296         hipline->temperature = 75;
 297     }
 298     else if (strncmp(hipline->sp_type, "M4", 2) == 0) {
 299         hipline->temperature = 75;
 300     }
 301     else if (strncmp(hipline->sp_type, "M5", 2) == 0) {
 302         hipline->temperature = 75;
 303     }
 304     else if (strncmp(hipline->sp_type, "M6", 2) == 0) {
 305         hipline->temperature = 75;
 306     }
 307     else if (strncmp(hipline->sp_type, "M7", 2) == 0) {
 308         hipline->temperature = 75;
 309     }
 310     else if (strncmp(hipline->sp_type, "M8", 2) == 0) {
 311         hipline->temperature = 75;
 312     }
 313     else if (strncmp(hipline->sp_type, "M9", 2) == 0) {
 314         hipline->temperature = 75;
 315     }
 316     else if (strncmp(hipline->sp_type, "O", 1) == 0) {
 317         hipline->temperature = 13;
 318     }
 319     else if (strncmp(hipline->sp_type, "B", 1) == 0) {
 320         hipline->temperature = 35;
 321     }
 322     else if (strncmp(hipline->sp_type, "M", 1) == 0) {
 323         hipline->temperature = 75;
 324     }
 325     else if (strncmp(hipline->sp_type, "C", 1) == 0) {
 326         hipline->temperature = 75;
 327     }
 328     else if (strncmp(hipline->sp_type, "A", 1) == 0) {
 329         hipline->temperature = 56;
 330     }
 331     else if (strncmp(hipline->sp_type, "R", 1) == 0) {
 332         hipline->temperature = 75;
 333     }
 334     else if (strncmp(hipline->sp_type, "G", 1) == 0) {
 335         hipline->temperature = 66;
 336     }
 337     else if (strncmp(hipline->sp_type, "W", 1) == 0) {
 338         hipline->temperature = 35;
 339     }
 340     else if (strncmp(hipline->sp_type, "K", 1) == 0) {
 341         hipline->temperature = 72;
 342     }
 343     else if (strncmp(hipline->sp_type, "N", 1) == 0) {
 344         hipline->temperature = 72;
 345     }
 346     else if (strncmp(hipline->sp_type, "S", 1) == 0) {
 347         hipline->temperature = 72;
 348     }
 349     else if (strncmp(hipline->sp_type, "F", 1) == 0) {
 350         hipline->temperature = 63;
 351     }
 352     else if (strncmp(hipline->sp_type, "DA", 2) == 0) {
 353         hipline->temperature = 35;
 354     }
 355     /*Special Cases*/
 356     else {
 357         printf("Unimplemented spectral type: %s\n", hipline->sp_type);
 358         hipline->temperature = 66;
 359     }
 360     return (EXIT_SUCCESS);
 361 }
 362 /**************************************************************************/
 363
 364 int READ_CROSS_SEC(struct SPECTRA line, struct INP_PAR file_list)
 365 {
 366     FILE    *fp;
 367     float   wave, cross_sec, t1, t2, t3;
 368     char    str[132], *read_res;
 369     int i = 0, j;
 370     struct SPECTRA temp;
 371
 372     /*  The following steps are to find out the Absorption CrossSection from
 373      the cross-section file corresponding to our wavelength of interest. This
 374      file was taken directly from Draine's page:
 375      ftp://ftp.astro.princeton.edu/draine/dust/mix/kext_albedo_WD_MW_3.1_60_D03.all
 376      */
 377
 378     if ((fp = fopen(file_list.sigma_file, "r")) == NULL){
 379         printf("Error opening Cross Section File:::\n");
 380         exit(EXIT_FAILURE);
 381     }
 382     temp.spectrum = malloc(sizeof(float) * MAX_SPEC_ELEM);
 383     temp.wavelength = malloc(sizeof(float) * MAX_SPEC_ELEM);
 384
 385     /*Skip all lines until we reach -*/
 386     strcpy(str, "          ");
 387     while ((strncmp("-", str, 1) != 0)) {
 388         read_res = fgets(str, 132, fp);
 389     }
 390     /*Find the cross section at the nearest wavelength*/
 391     i = 0;
 392     while (feof(fp) == 0){
 393         read_res = fgets(str, 132, fp);/*Read first bracket*/
 394         sscanf(str, "%f %f %f %e %e",
 395                &wave, &t1, &t2, &cross_sec, &t3);
 396         temp.wavelength[i] = wave*10000;/*convert to Å*/
 397         temp.spectrum[i]   = cross_sec;
 398         i = i + 1;
 399     }
 400     /*We have to reverse the wavelengths because of Draine's format*/
 401     for (j=0; j < i; ++j) {
 402         line.wavelength[j] = temp.wavelength[i-j-1];
 403         line.spectrum[j] = temp.spectrum[i-j-1];
 404     }
 405     fclose(fp);
 406     free(temp.spectrum);
 407     free(temp.wavelength);
 408     return (i);
 409 }
 410 /**************************************************************************/
 411
 412 float CALC_FLUX(struct STARS *hipstars, struct INP_PAR files,
 413                 struct SPECTRA cross_sec,
 414                 struct SPECTRA *stellar_spectra,
 415                 struct SPECTRA filters,
 416                 float gl, float gb, float dist, float lambda,
 417                 float ave_sigma)
 418 {
 419     int istars, windex, sindex, cindex;
 420     float x, y, z;
 421     float dstarsq, tau, flux = 0;
 422     struct SPECTRA star_flux;
 423
 424     /*Distance from star to point of interest*/
 425     x = cos(DEGRAD(gl))*cos(DEGRAD(gb))*dist;
 426     y = sin(DEGRAD(gl))*cos(DEGRAD(gb))*dist;
 427     z = sin(DEGRAD(gb))*dist;
 428
 429     star_flux.wavelength = (float *) malloc(sizeof(float) * N_CASTELLI_SPECTRA);
 430     star_flux.spectrum   = (float *) malloc(sizeof(float) * N_CASTELLI_SPECTRA);
 431     star_flux.nelements  = N_CASTELLI_SPECTRA;
 432
 433     istars = 0;
 434     dstarsq = pow((hipstars[istars].x - x),2) +
 435     pow((hipstars[istars].y - y),2) +
 436     pow((hipstars[istars].z - z),2);
 437     sindex = hipstars[istars].temperature;
 438     cindex = 0;
 439     for (windex = 0; windex < star_flux.nelements; ++windex) {
 440         star_flux.wavelength[windex] = stellar_spectra[sindex].wavelength[windex];
 441         while ((cross_sec.wavelength[cindex] < star_flux.wavelength[windex])
 442                && (cindex < cross_sec.nelements))
 443             ++cindex;
 444         tau = cross_sec.spectrum[cindex]*sqrt(dstarsq)*ave_sigma*PCinCM(1.);
 445         star_flux.spectrum[windex] =
 446         stellar_spectra[sindex].spectrum[windex] * expf(-tau) *
 447         hipstars[istars].scale/dstarsq*
 448         1e9*stellar_spectra[sindex].wavelength[windex]/6.626/3.;
 449     }
 450     flux = CALC_SCALE_FACTOR(filters, star_flux);
 451     return(flux);
 452 }
 453
 454 int GET_SCALE_FACTOR(struct STARS *hipstars,
 455                      struct SPECTRA *stellar_spectra)
 456 {
 457     float scale;
 458     int bindex, vindex, sindex, doprint;
 459     float ebv, b_mag, v_mag, bflux, vflux, bv;
 460
 461     if (hipstars->V_mag == 0){
 462         hipstars->scale = 0;/* A few stars are bad*/
 463     } else {
 464         sindex = hipstars->temperature;
 465         /*We scale at 4400 and 5500 Å*/
 466         vindex = 0;
 467         while (stellar_spectra[sindex].wavelength[vindex] < 5500)
 468             ++vindex;
 469         bindex = 0;
 470         while (stellar_spectra[sindex].wavelength[bindex] < 4400)
 471             ++bindex;
 472
 473         bflux = stellar_spectra[sindex].spectrum[bindex];
 474         vflux = stellar_spectra[sindex].spectrum[vindex];
 475         b_mag = -2.5 * log10(bflux/6.61);/*Convert flux into magnitudes*/
 476         v_mag = -2.5 * log10(vflux/3.64);
 477         bv = hipstars->B_mag - hipstars->V_mag;
 478         ebv = bv - (b_mag - v_mag);
 479         if (ebv < 0) ebv = 0;
 480         hipstars->ebv = ebv;
 481         /*Scale factor to convert model to Earth with no extinction*/
 482         scale = 3.64e-9*pow(10, -0.4 * (hipstars->V_mag - 3.2*ebv))/
 483         vflux;
 484         /*Scale to a distance of 1 pc*/
 485         hipstars->scale = scale*hipstars->distance*hipstars->distance;
 486
 487         doprint = 0;
 488         if (doprint > 0)
 489             printf("%li %s %s %i %10.3e %10.3e\n", hipstars->HIP_NO, hipstars->sp_type,
 490                    stellar_spectra[sindex].filename, sindex, bflux, vflux);
 491     }
 492
 493     return(EXIT_SUCCESS);
 494 }
 495
 496 int READ_HIPPARCOS_CAT(struct INP_PAR inp_par, struct SPECTRA *stellar_spectra,
 497                        struct STARS *hipstars, struct SPECTRA filter)
 498 {
 499     FILE *hipfile_ptr=NULL;
 500     int i, tst;
 501     struct STARS hipline;
 502     struct SPECTRA cross_sec;
 503
 504     tst = READ_CASTELLI_SPECTRA(inp_par.castelli_file, stellar_spectra);
 505     cross_sec.spectrum = malloc(sizeof(float) * MAX_SPEC_ELEM);
 506     cross_sec.wavelength = malloc(sizeof(float) * MAX_SPEC_ELEM);
 507     cross_sec.nelements = READ_CROSS_SEC(cross_sec, inp_par);
 508
 509     hipfile_ptr = fopen(inp_par.hipparcos_file, "r");
 510     if (hipfile_ptr == NULL){
 511         printf("%s %s\n", "Can't find stellar file:", inp_par.hipparcos_file);
 512         exit(EXIT_FAILURE);
 513     }
 514     i = 0;
 515     while (feof(hipfile_ptr) == 0) {
 516         tst   = HIP_READ_LINE(hipfile_ptr, &hipline);
 517         tst = GET_STAR_TEMP(&hipline);
 518         tst = GET_SCALE_FACTOR(&hipline, stellar_spectra);
 519         hipline.flux_at_earth = CALC_FLUX(&hipline, inp_par, cross_sec, stellar_spectra, filter,
 520                                           0, 0, 0, 2300, 1);
 521         hipstars[i] = hipline;
 522         ++i;
 523     }
 524     hipstars[0].nstars = i-1;
 525     fclose(hipfile_ptr);
 526     return (EXIT_SUCCESS);
 527
 528 }