# common import statements import os import numpy as np from scipy.io.idl import readsav print "This package contains the following functions:" print "readWhelanModels reads models by Whelan et al 2011 and stores them in a record" ############################ # This routine reads all the models from D. Whelan ############################ def readWhelanModels(help=False,jy=True,location='none'): """ This routine reads all the models produced by Dave Whelan The parameters describing the models are: - The star formation efficiency - the size of the outer radius - The clumpyness of the shell - The size of the inner radius - The viewing angle The first four parameters define the name of a file (i.e. the number of files is equal to the product of the number of values for these parameters) while the viewing angle is varied inside each file. The viewing angles is defined by 10 values of theta and 20 values of phi. The spectrum is defined by 200 values of lambda (mic). In the models, the wavelength is tabulated in microns, while the output flux is tabulated in erg.s-1.cm-2 and the flux is computed assuming the source is at 10 Mpc. I want to output the model flux in Jy which means that I need to multiply the tabulated fluxes by (10^9)*lam[mic]/3.0. """ # HISTORY # 03-Oct-2011 MS Convert the model flux to Jy # 25-Oct-2011 MS In the conversion to Jy I have established that the flux density (in Jy) is # (1e9/3)*lam[mic]*(lamflam[erg.s-1.cm-2]. Note that there is a flux discrepancy between # the models as published and as available on the site. Dave Whelan is aware of this. if (help): print "result = readWhelanModels(help=T/F,jy=T/F,location=location)" print "" print "Reads all the Whelan models and returns them" print "" print "help if True prints this" print "jy if True (default) converts the spectra in Jy, otherwise" print " leaves them in their original units (erg.s-1.cm-2)" print "location a string to indicate the location of the model files" return # check whether location is specified if (location == 'none'): # the location was not provided so we use the default print "There is little chance that the default location is meaningful for you..." location = "/Users/sauvage/DangerZone/Temporary_Data/Whelan/" # first define the arrays of the variables are are constant sfe = np.array([5,10,15,25,50]) sfeStr = ['05','10','15','25','50'] rout = np.array([25,50,100]) clump = np.array([0.01,0.1,0.5,0.9,1.0]) # defines the dimensions of the "inner" variables: # number of lines of sight nLos = 200 # number of sampled wavelengths nLam = 200 # Create the datatype that will hold the models it consists of 4 float fields that will contain # the values of the SFE, Rout, Clump and Rin that index a model uniquely as well as two 2-D array # that contains the values of lambda and lamflam for each node of the line of sight and lambda grid # I may reduce the lambda array as all of the spectra share the same array. modDType = np.dtype([('sfe',float),('rout',float),('rin',float),('fclump',float),\ ('specgrid',float,(nLos,nLam)),('lamgrid',float,(nLam))]) firstFile = True for iSfe in range(sfe.size): for iRout in range(rout.size): for iClump in range(clump.size): # first build the directory name from the main three variables dirName = 'sfe'+sfeStr[iSfe]+'_rout'+str(rout[iRout])+'_f'+str(clump[iClump])+'/' # check that the directory exists and counts its countent if (os.path.exists(location+dirName)): # there are files in this directory theFiles = os.listdir(location+dirName) # read the files for iFile in range(len(theFiles)): # instantiate an object to hold the data thisModel = np.zeros(1,dtype=modDType) thisModel['sfe'] = sfe[iSfe] thisModel['rout'] = rout[iRout] thisModel['fclump'] = clump[iClump] # first I need to extract the value of rin from the name name = 'cl'+sfeStr[iSfe]+'_r'+str(rout[iRout])+'_f'+str(clump[iClump])+'.' rInStr = theFiles[iFile].replace(name,'').replace('_spectra.txt','') # I now have the value of rIn for this spectral family rin = float(rInStr) print rin thisModel['rin'] = rin # now read the whole file in one go (thanks Pierre Chanial for the trick!) data = np.loadtxt(location+dirName+theFiles[iFile], dtype=[('lam',float), ('lamflam',float)]) # now data is a 1D array of pairs of values data['lam'] contains the value of lambda in microns # data['lamflam'] contains the flux values in units of erg s-1 cm-2 #specGrid = np.ndarray(shape=(nLos,nLam),dtype=float) #lamGrid = np.ndarray(shape=(nLos,nLam),dtype=float) #use the first lambda grid, they are all the same thisModel['lamgrid'][0][:] = data['lam'][0:nLam] for iLos in range(nLos): if (jy): # on the fly I convert the model to Jy thisModel['specgrid'][0][iLos,:] = data['lam'][0:nLam]*data['lamflam'][iLos*nLam:(iLos+1)*nLam]*(1e9)/3.0 else: thisModel['specgrid'][0][iLos,:] = data['lamflam'][iLos*nLam:(iLos+1)*nLam] # now see if we need to concatenate this object to already read ones if (firstFile): # this is the first time we read a model result = thisModel.copy() firstFile = False else: result = np.concatenate([result,thisModel]) else: print dirName+' does not exist' # Print some succint help print "Models are stored in a numpy 1D array, with elements that are of composite datatype" print "For each element you will find the value of the model parameters, the wavelength grid" print "and the 200 computed line of sight SEDs" print "To start navigating this array type myArray.dtype (where myArray is the output variable" print "of this function)" return result