Source Code for Module ComboCode.cc.tools.io.Atmosphere

  1  # -*- coding: utf-8 -*- 
  2   
  3  """ 
  4  Reading stellar model atmospheres for use with MCMax and GASTRoNOoM. 
  5   
  6  Author: R. Lombaert 
  7   
  8  """ 
  9   
 10  from glob import glob 
 11  from astropy.io import fits as pyfits 
 12  from scipy import rec,array,argmin 
 13  import os 
 14   
 15  import cc.path 
 16   
 17 -class Atmosphere(object): 
 18       
 19      """ 
 20      Reading stellar model atmospheres. 
 21       
 22      """ 
 23       
 24 -    def __init__(self,modeltype,filename=None): 
 25           
 26          """ 
 27          Initializing an Atmosphere() object.  
 28           
 29          The stellar model atmospheres are read from a given folder and the  
 30          preferred model is chosen after initialisation. 
 31           
 32          The location of the model atmospheres is taken from Path.dat. 
 33   
 34          @param modeltype: Type of model atmosphere (any from IVSDATA folder) 
 35          @type modeltype: string 
 36           
 37          @keyword filename: The filename of the model atmosphere in the  
 38                             modelgrids folder. If None, the filename can be set 
 39                             later. 
 40                              
 41                             (default: None) 
 42          @type filename: string 
 43                   
 44          """ 
 45           
 46          self.modeltype = modeltype 
 47          self.filename = filename 
 48          if not self.filename is None: 
 49              self.filepath = os.path.join(cc.path.atm,self.filename) 
 50          self.modellist = glob(os.path.join(cc.path.atm,modeltype+'*')) 
 51          self.modelgrid = None 
 52          self.teff_actual = None 
 53          self.logg_actual = None 
 54   
 55   
 56   
 57 -    def getModelList(self): 
 58           
 59          """ 
 60          Retrieve a list of models available in the modelgrids folder 
 61          for the requested model type. 
 62           
 63          When setting a model, the model is identified by its index in this list. 
 64           
 65          @return: The filepaths available for given modeltype.  
 66          @rtype: list[string] 
 67           
 68          """ 
 69           
 70          return self.modellist 
 71           
 72       
 73       
 74 -    def setModelGrid(self,index): 
 75           
 76          """ 
 77          Set the model atmosphere to be read, identified by the index in the  
 78          model list. 
 79           
 80          Once set here, or if given upon initialisation of the Atmosphere()  
 81          object, the model cannot be changed for this instance. 
 82           
 83          """ 
 84           
 85          if self.filename is None: 
 86              self.filename = os.path.split(self.modellist[index])[1] 
 87              self.filepath = os.path.join(cc.path.atm,self.filename) 
 88   
 89   
 90       
 91 -    def readModelGrid(self): 
 92           
 93          """ 
 94          Read the model atmosphere fits file. 
 95           
 96          """ 
 97           
 98          self.ff = pyfits.open(self.filepath) 
 99          self.header = self.ff[0].header 
100          teffs = array([self.ff[i].header['TEFF']  
101                         for i in xrange(1,len(self.ff))]) 
102          loggs = array([self.ff[i].header['LOGG']  
103                         for i in xrange(1,len(self.ff))]) 
104          self.modelgrid = rec.fromarrays([array(range(1,len(self.ff))),teffs,\ 
105                                           loggs],\ 
106                                          names=['INDEX','TEFF','LOGG']) 
107                                           
108       
109       
110 -    def getHeader(self): 
111           
112          """ 
113          Return the main header of the fits file. 
114           
115          @return: The header information 
116          @rtype: pyfits.NP_pyfits.Header() 
117           
118          """ 
119           
120          return self.header 
121           
122   
123   
124 -    def getModelGrid(self): 
125           
126          """ 
127          Return the TEFF and LOGG grid for this model.  
128           
129          @return: the grid values of the teff and logg parameters 
130          @rtype: recarray 
131           
132          """ 
133           
134          if self.modelgrid is None: self.readModelGrid() 
135          return self.modelgrid 
136           
137           
138           
139 -    def getModel(self,teff,logg): 
140           
141          """ 
142          Return the model atmosphere for given effective temperature and log g. 
143           
144          Not yet scaled to the distance! 
145           
146          Units returned are (micron,Jy) 
147           
148          @param teff: the stellar effective temperature 
149          @type teff: float 
150          @param logg: the log g value 
151          @type logg: float 
152           
153          @return: The model spectrum in (micron,Jy) 
154          @rtype: recarray 
155           
156          """ 
157           
158          c = 2.99792458e18          #in angstrom/s 
159          if self.modelgrid is None: 
160              self.readModelGrid() 
161          mg = self.modelgrid 
162          #- Find the closest temperature in the grid 
163          teff_prox = mg['TEFF'][argmin(abs(mg['TEFF']-teff))] 
164          #- Select all models with that temperature 
165          mgsel = mg[mg['TEFF']==teff_prox] 
166          #- Select the closest log g in the selection 
167          logg_prox = mgsel['LOGG'][argmin(abs(mgsel['LOGG']-logg))] 
168          #- Get the index of the model closest to teff and logg 
169          imodel = mgsel[mgsel['LOGG']==logg_prox]['INDEX'][0] 
170           
171          self.teff_actual = teff_prox 
172          self.logg_actual = logg_prox         
173           
174          wave = self.ff[imodel].data.field('wavelength') 
175          flux = self.ff[imodel].data.field('flux') 
176          if self.header['FLXUNIT'] == 'erg/s/cm2/A': 
177              #- Go to erg/s/cm2/Hz, lFl = nFn, then to Jy (factor 10**(23)) 
178              flux = flux * wave**2 / c * 10**(23) 
179          else: 
180              raise Error('Flux unit unknown in atmosphere model fits file.') 
181          if self.header['WAVUNIT'] == 'angstrom': 
182              wave = wave * 10**(-4) 
183          else: 
184              raise Error('Wavelength unit unknown in atmosphere model fits file.') 
185           
186          model = rec.fromarrays([wave,flux],names=['wave','flux'])         
187          return model 
188