ComboCode.cc.data.Sed

1 # -*- coding: utf-8 -*- 2 3 """ 4 Toolbox for reading spectra and spectrum manipulation. 5 6 Author: R. Lombaert 7 8 """ 9 10 import os 11 from numpy import argsort,array 12 from scipy.integrate import trapz 13 from scipy.interpolate import interp1d 14 from numpy.core.defchararray import rfind,ljust 15 import numpy.lib.recfunctions as recfunc 16 import numpy as np 17 from glob import glob 18 import operator 19 20 from astropy import units as u 21 22 from cc.ivs.sed import builder, filters 23 import cc.ivs.sed.reddening as ivs_red 24 from cc.ivs.sed.model import synthetic_flux 25 #from ivs.units import conversions 26 27 import cc.path 28 from cc.tools.io import DataIO 29 from cc.modeling.tools import Reddening 30 from cc.modeling.codes import MCMax 31 32

33 -def getCFlux(wav,seds=[],star_grid=[],nans=1,deredden=[],\ 34 law='Fitz2004Chiar2006',lawtype='ism',map='marshall'):

35 36 ''' 37 Retrieve the continuum flux at a given wavelength from either a model 38 spectrum or an observation. If both seds and star_grid are given, 39 values from the seds are returned first in the array, then the models. 40 41 For now assumes the observation is either an ISO SWS spectrum OR that the 42 continuum point is given in a dictionary that is property of the Sed() 43 object (sed.cflux) with wavelengths as keys. 44 45 star_grid are all models! Works also when seds or star_grid are empty. 46 47 Reddening is taken into account when requested in the models and parameters 48 are taken from the model objects. However, this is only allowed if only one 49 data object is given (otherwise model reddening doesn't make sense) 50 51 Dereddening of data is also possible (and extra arguments can be passed to 52 the reddening law), in which case distances have to be given for the seds. 53 If any model reddening is requested and only sed is given, sed dereddening 54 is always turned off. 55 56 @param wav: The continuum wavelength point 57 @type wav: float 58 59 @keyword seds: The SEDs of the data objects. Number of SEDs sets the amount 60 of Star() objects represent data versus number of models. 61 62 (default: []) 63 @type seds: list(Sed()) 64 @keyword star_grid: The data + model objects 65 66 (default: []) 67 @type star_grid: list(Star()) 68 @keyword nans: Set undefined line strengths as nans. Errors are set as a 69 nan if it concerns mode==dint. Otherwise, they are not set. 70 71 (default: 1) 72 @type nans: bool 73 @keyword deredden: Deredden the SEDs with distances given here. This option 74 is turned off automatically if any reddening is requested 75 in the models and only one sed is given to avoid double 76 correction. Number of distances given must be equal to 77 number of SEDs given. If not, it is also turned off. 78 79 (default: []) 80 @type deredden: list 81 @keyword law: The reddening law for DEREDDENING 82 83 (default: 'Fitz2004Chiar2006') 84 @type law: str 85 @keyword lawtype: The type of Chiar & Tielens reddening law (either ism or 86 gc) for DEREDDENING 87 88 (default: 'ism') 89 @type lawtype: str 90 @keyword map: The galactic 3d extinction model for DEREDDENING. 91 92 (default: 'marshall') 93 @type map: str 94 95 @return: The continuum fluxes in W/m2/Hz with length that of star_grid, as 96 well as errors if applicable. If a combo mode is requested, errors 97 are given when available, and listed as None/nan if not available 98 (Plotting2 module knows how to deal with this). 99 @rtype: (array[float],array[float]) 100 101 ''' 102 103 all_cflux = [] 104 all_eflux = [] 105 rlaw = [] 106 107 dtype = '' 108 if seds: 109 ddict = dict([('SWS',(2.4,45.0)),('PACS',(55.1,189.))]) 110 for k,v in ddict.items(): 111 if v[0] <= wav and wav <= v[1]: 112 dtype = k 113 break 114 115 #-- Only allow model reddening if one sed is given. Multiple: makes no sense 116 # to redden models. None: No galactic coordinates given to redden. 117 redden = [s['REDDENING'] for s in star_grid] if len(seds) == 1 else [] 118 119 #-- Only allow dereddening if enough distances are given, and if model 120 # reddening for one sed is not requested. 121 if len(deredden) != len(seds) or np.any(redden): 122 deredden = [] 123 124 #-- Interpolate the reddening law once. 125 if deredden or redden: 126 wave_arr,rlaw = ivs_red.get_law(name=law,wave=wav,curve=lawtype,\ 127 norm='Ak',wave_units='micron') 128 129 #-- First all data objects 130 for ised,sed in enumerate(seds): 131 dtypes = [] 132 if dtype: 133 dtypes = [dt for dt in sed.data.keys() if dtype in dt[0].upper()] 134 #-- If no SWS spectrum find, check if the flux is available in sed.flux 135 if not dtypes: 136 if sed.cflux.has_key(wav): 137 tflux = sed.cflux[wav] 138 if deredden: 139 ak = sed.getAk(deredden[ised],map=map,law=law) 140 #-- deredden so increase flux, as opposed to redden 141 tflux = tflux * 10**(rlaw*ak/2.5) 142 all_cflux.append(tflux) 143 all_eflux.append(sed.eflux[wav]) 144 else: 145 all_cflux.append(nans and float('nan') or None) 146 all_eflux.append(nans and float('nan') or None) 147 continue 148 #-- At least one dtype spectrum found, take the first one. 149 dt = dtypes[0] 150 abs_err = sed.abs_err[dt[0]] 151 dwave = sed.data[dt][0] 152 dflux = sed.data[dt][1] 153 154 #-- Check if the data object gives the standard deviation 155 if len(sed.data[dt]) > 2: 156 i = np.argmin(abs(dwave-wav)) 157 ilow = i if wav>dwave[i] else i-1 158 iup = i if wav<dwave[i] else i+1 159 errs = sed.data[dt][2] 160 deflux = np.sqrt((errs/dflux)[ilow]**2+(errs/dflux)[iup]**2) 161 else: 162 deflux = 0.0 163 164 #-- Interpolate for flux, and set the error taking into account abs flux 165 # calib uncert. 166 interp = interp1d(dwave,dflux) 167 tflux = interp(wav) 168 if deredden: 169 ak = sed.getAk(deredden[ised],map=map,law=law) 170 #-- deredden so increase flux 171 tflux = tflux * 10**(rlaw*ak/2.5) 172 all_cflux.append(tflux) 173 all_eflux.append(np.sqrt(deflux**2+abs_err**2)) 174 175 #-- Then all model objects 176 all_eflux.extend([nans and float('nan') or None]*len(star_grid)) 177 for s in star_grid: 178 if not s['LAST_MCMAX_MODEL']: 179 all_cflux.append(nans and float('nan') or None) 180 continue 181 cc.path.mout = os.path.join(cc.path.mcmax,s.path_mcmax) 182 dpath = os.path.join(cc.path.mcmax,s.path_mcmax,'models',\ 183 s['LAST_MCMAX_MODEL']) 184 w,f = MCMax.readModelSpectrum(dpath,rt_sed=1) 185 interp = interp1d(w,f) 186 tflux = interp(wav) 187 if s['REDDENING'] and seds: 188 #-- Only one SED is supposed to be given. 189 ak = seds[0].getAk(s['DISTANCE'],map=s['REDDENING_MAP'],\ 190 law=s['REDDENING_LAW']) 191 #-- redden so decrease flux 192 tflux = tflux / 10**(rlaw*ak/2.5) 193 all_cflux.append(tflux) 194 195 #-- All fluxes for SED type data or models are given in Jy. Convert to 196 # W/m2/Hz. Errors are given in relative numbers, so no conversion needed. 197 all_cflux = array(all_cflux)*1e-26 198 all_eflux = array(all_eflux) 199 200 return (all_cflux,all_eflux)

201 202 203

204 -def calcPhotometry(w,f,photbands):

205 206 ''' 207 Calculate the (model) photometry for a given set of photometric bands. 208 209 Reddening is assumed to have been done before this. 210 211 @param w: the wavelengths in micron 212 @type w: array() 213 @param f: Flux grid in Jy 214 @type f: array() 215 @param photbands: the photometric bands 216 @type photbands: array(str) 217 218 @return: The photometry in Jy 219 @rtype: array 220 221 ''' 222 223 #-- Convert wavelength from micron to angstrom. astropy conversion module 224 # returns a "Quantity" that has properties "unit" and "value" 225 mlam = (w*u.micron).to(u.AA).value 226 227 #-- Convert Jy flux density to erg/s/cm2/aa for windowed integration 228 # Requires equivalency because conversion depends on the wavelength where 229 # the flux is measured. Uses spectral_density equivalency. 230 mflam = (f*u.Jy).to(u.erg/u.s/u.cm**2/u.AA,\ 231 equivalencies=u.spectral_density(w*u.micron)).value 232 233 # mflam = conversions.convert('Jy','erg/s/cm2/AA',f,wave=(w,'micron')) 234 # mlam = conversions.convert('micron','AA',w) 235 mphot = synthetic_flux(mlam,mflam,photbands,units=['Fnu']*len(photbands)) 236 237 mphot = (mphot*u.erg/u.s/u.Hz/u.cm**2).to(u.Jy).value 238 # mphot = conversions.convert('erg/s/cm2/Hz','Jy',mphot) 239 return mphot

240 241 242

243 -def buildPhotometry(star_name,fn='Photometric_IvS',remove=[]):

244 ''' 245 Retrieve the photometry of a star through the IvS repo's SED builder. 246 247 Save the photometry in a target location. 248 249 @param star_name: Name of the star (cc name from Star.dat) 250 @type star_name: str 251 252 @keyword fn: Output filename of the photometry file. Always appends 253 '_STAR.dat' where STAR is star_name. The file is saved in dp. 254 255 (default: 'Photometric_IvS') 256 @type fn: str 257 @keyword remove: Photometry to be removed from the output file. 258 e.g. ['WISE','DENIS'] 259 260 (default: []) 261 @type remove: list[str] 262 263 ''' 264 265 #-- Get the SIMBAD name of the star 266 si = DataIO.getInputData(filename='Star.dat').index(star_name) 267 sn_sim = DataIO.getInputData(keyword='STAR_NAME_PLOTS')[si] 268 sn_sim = sn_sim.replace('$','').replace('\\','') 269 270 #-- Define the outputfolder of the Raw data 271 ofn_raw = os.path.join(cc.path.dphot,'_'.join([star_name,'phot.txt'])) 272 273 ## Retrieve p = builder.SED() 274 star = builder.SED(sn_sim,photfile=ofn_raw) 275 star.get_photometry() 276 277 ## Read in photometric data 278 columns = ['wave', 'meas', 'emeas', 'photband', 'bibcode', 'comments'] 279 data = np.genfromtxt(star.photfile,usecols = (9,10,11,4,15,16),\ 280 names=columns,dtype = None) 281 282 #-- Find and remove WISE and DENIS data 283 for photband in remove: 284 selection = rfind(data['photband'],photband) 285 data = data[np.where(selection==-1)] 286 287 #-- Remove colour measurements (NaN in wavelength) 288 data = data[np.where(np.isfinite(data['wave']))] 289 290 #-- Fill photbands with whitespace for formatting 291 data['photband'] = ljust(data['photband'],15) 292 293 #-- Flux: erg/s/cm**2/AA -> Jy (10**-23 erg/s/cm**2/Hz) 294 # From http://astro.wku.edu/strolger/UNITS.txt (non-linear!) 295 c = 2.99792458e18 # in AA/s 296 data['meas'] = data['meas']*data['wave']**2/c*1e23 297 data['emeas'] = data['emeas']*data['wave']**2/c*1e23 298 299 #-- Wavelength: angstrom -> micron 300 data['wave'] = data['wave']*10**-4 301 302 #-- Sort the data on wavelength 303 data.sort() 304 305 ofn_final = os.path.join(cc.path.dsed,'_'.join([fn,star_name+'.dat'])) 306 hdr = 'Photometry extracted with ivs.sed.builder. \nWave (micron) Flux '+\ 307 '(Jy) Error Flux (Jy) Photband Bibcode Comments' 308 np.savetxt(ofn_final,data,fmt=['%.8e']*3+['%s']*3,header=hdr)

309 310 311

312 -def normalizeSed(sed):

313 314 ''' 315 Normalize an SED by division by integrated energy. 316 317 Input is given by F_nu, integration of the SED is done for nu*F_nu. 318 319 @param sed: The input SED wavelenth and F_nu. 320 @type sed: (list,list) 321 322 ''' 323 324 c = 2.99792458e14 #in micron/s 325 326 return [sed[0],sed[1]/trapz(x=sed[0],y=sed[1]*(c/sed[0]))]

327 328 329

330 -class Sed(object):

331 332 ''' 333 Tools for: 334 - reading data 335 - downloading photometry 336 337 ''' 338

339 - def __init__(self,star_name,remove=[]):

340 341 ''' 342 Initializing an Sed instance. 343 344 Setting starting parameters from a star object. 345 346 @param star_name: The star name of the object 347 @type star_name: string 348 349 @keyword remove: Photometry to be removed from the output file. 350 e.g. ['WISE','DENIS'] 351 352 (default: []) 353 @type remove: list[str] 354 355 ''' 356 357 self.instrument = 'SED' 358 self.photbands = np.empty(0,dtype=str) 359 self.photwave = np.empty(0) 360 self.star_name = star_name 361 self.data = dict() 362 self.setStarPars() 363 self.setData(remove=remove) 364 self.readData() 365 self.readPhotInfo() 366 self.ak = dict() 367 368 #-- Extra property defining continuum flux points that can be set 369 # externally. Assumed to give (wavelength,flux) as (key,value) pair 370 # Used by, e.g., the Sed.getCFlux() method. An error must also be 371 # given (in relative number). 372 self.cflux = dict() 373 self.eflux = dict()

374 375

376 - def setStarPars(self):

377 378 """ 379 Set some standard stellar parameters such as Ak and galactic position. 380 381 """ 382 383 self.star_index = DataIO.getInputData().index(self.star_name) 384 self.ll = DataIO.getInputData(keyword='LONG',rindex=self.star_index) 385 self.bb = DataIO.getInputData(keyword='LAT',rindex=self.star_index) 386 snp = DataIO.getInputData(keyword='STAR_NAME_PLOTS',\ 387 remove_underscore=1,rindex=self.star_index) 388 self.star_name_plots = snp

389 390 391

392 - def setData(self,**kwargs):

393 394 ''' 395 Select available data. 396 397 Based on the data file types in Sed.dat and the available data files. 398 399 Also calls the buildPhotometry method to create a photometry file from 400 the IvS Sed builder tool 401 402 Any keywords required for buildPhotometry can be passed here. 403 404 ''' 405 406 data_types = DataIO.getInputData(keyword='DATA_TYPES',\ 407 filename='Sed.dat') 408 abs_errs = DataIO.getInputData(keyword='ABS_ERR',filename='Sed.dat') 409 410 if 'Photometric_IvS' in data_types: 411 buildPhotometry(self.star_name,**kwargs) 412 413 self.data_types = [] 414 self.data_filenames = [] 415 self.abs_err = dict() 416 for dt,ierr in zip(data_types,abs_errs): 417 searchpath = os.path.join(cc.path.dsed,'%s_*%s*.dat'\ 418 %(dt,self.star_name)) 419 add_files = glob(searchpath) 420 for ff in add_files: 421 if ff not in self.data_filenames: 422 self.data_filenames.append(ff) 423 self.data_types.append(dt) 424 self.abs_err[dt] = ierr

425 426 427

428 - def readData(self):

429 430 ''' 431 Read the raw SED data. 432 433 ''' 434 435 for dt,fn in zip(self.data_types,self.data_filenames): 436 data = DataIO.readCols(fn,nans=0) 437 #-- Currently, error bars only available for these types of data. 438 if 'Photometric' in dt or 'MIDI' in dt or 'Sacha' in dt: 439 #-- Sort MIDI data 440 if 'MIDI' in dt: 441 cdat = [dd[(data[0]<=13.)*(data[0]>=8.)] for dd in data] 442 i = argsort(cdat[0]) 443 self.data[(dt,fn)] = (cdat[0][i],cdat[1][i],cdat[2][i]) 444 else: 445 self.data[(dt,fn)] = (data[0],data[1],data[2]) 446 if dt == 'Photometric_IvS': 447 self.photbands = data[3] 448 self.photwave = data[0] 449 else: 450 #-- Still sorting for PACS. Obsolete when separate bands for 451 # PACS are available. 452 i = argsort(data[0]) 453 self.data[(dt,fn)] = (data[0][i],data[1][i])

454 455 456

457 - def readPhotInfo(self,level=.5):

458 459 ''' 460 Read the photometry band information associated with photometry of this 461 SED. 462 463 @keyword level: The level at which the cut off for significant 464 transmission of the photometric bands is placed. 465 466 (default: 0.5) 467 @type level: float 468 469 ''' 470 471 #-- Get photometry bands info from IvS repo. recarray structure same as 472 # self.photbands_ivs 473 filter_info = filters.get_info() 474 keep = np.searchsorted(filter_info['photband'],self.photbands) 475 self.filter_info = filter_info[keep] 476 self.filter_info.eff_wave = self.filter_info.eff_wave/1e4 477 478 response = [filters.get_response(photband) 479 for photband in self.photbands] 480 selection = [waver[transr/max(transr)>level]/1e4 481 for waver,transr in response] 482 wlower = [sel[0] for sel in selection] 483 wupper = [sel[-1] for sel in selection] 484 self.filter_info = recfunc.append_fields(self.filter_info,\ 485 ['wlower','wupper'],\ 486 [wlower,wupper],usemask=0,\ 487 asrecarray=1)

488

489 - def getAk(self,distance=None,map='marshall',law='Fitz2004Chiar2006'):

490 491 ''' 492 Helper method to retrieve the Ak extinction magnitude for a given 493 distance. 494 495 The Ak values are saved in the sed object, to cut down the overhead in 496 subsequent calls. 497 498 The law and map have to be passed as well. Defaults are marshall and 499 Fitz2004Chiar2006 respectively. 500 501 @param distance: The distance to the source. If the default, the total 502 extinction is calculated in given direction. 503 504 (default: None) 505 @type distance: float 506 @keyword map: The galactic 3d extinction model. 507 508 (default: 'marshall') 509 @type map: str 510 @keyword law: The reddening law 511 512 (default: 'Fitz2004Chiar2006') 513 @type law: str 514 515 @return: The extinction magnitude in K-band for the requested distance. 516 @rtype: float 517 518 ''' 519 520 distance = float(distance) 521 if not self.ak.has_key((distance,map,law)): 522 aki = Reddening.getAk(self.ll,self.bb,distance,map,law) 523 self.ak[(distance,map,law)] = aki 524 525 return self.ak[(distance,map,law)]

526

Source Code for Module ComboCode.cc.data.Sed