diff --git a/config.py b/config.py index a3546060f..e94f6c9b7 100644 --- a/config.py +++ b/config.py @@ -5,30 +5,30 @@ Usage: $ python config.py compile """ -import os -import collections -import ast +#-- Make sure the script doesn't confuse the IvS io module with the standard +# Python io module. import sys +sys.path = sys.path[1:] +import os import glob as glob_module #-- You can add directories here, but be sure that the relative paths within # those directories are correct! -data_dirs = [os.getenv('ivsdata'),'/STER/pieterd/IVSDATA/', '/STER/kristofs/IVSdata','/STER/jorisv/IVSDATA/', - '/STER/kenneth/Python_repository/','/home/ben/public_html/opacities','/STER/michelh/IVSDATA/'] - +data_dirs = [os.getenv('ivsdata'),'/STER/jorisv/IVSDATA/', + '/STER/kenneth/Python_repository/','/STER/michelh/IVSDATA/', + os.path.join(os.path.expanduser('~'),'ComboCode','aux','ivsdata')] + ivs_dirs = dict(coralie='/STER/coralie/', hermes='/STER/mercator/hermes/') - - def get_datafile(relative_path,basename): """ Reconstruct the location of a data file and check whether it exists. - + If the file exists, it will return the absolute path of the filename. - + If the file does not exist, it will raise an IOError. - + @param relative_path: relative path starting from main data directory tree @type relative_path: str @param basename: filename @@ -38,16 +38,16 @@ def get_datafile(relative_path,basename): """ for data_dir in data_dirs: if data_dir is None: continue - + filename = os.path.join(data_dir,relative_path,basename) - + if os.path.isfile(filename): break else: str_data_dirs = ", ".join([idir for idir in data_dirs if idir is not None]) relative_file = os.path.join(relative_path,basename) raise IOError, "File %s not found in any of the specified data directories %s"%(relative_file,str_data_dirs) - + return filename @@ -56,7 +56,7 @@ def get_datafile(relative_path,basename): def glob(relative_path,arg='*'): """ Glob the files in a relative path. - + @param relative_path: relative path starting from main data directory tree @type relative_path: str @param arg: argument to use in glob @@ -70,17 +70,15 @@ def glob(relative_path,arg='*'): files += glob_module.glob(os.path.join(data_dir,relative_path,arg)) files.sort() return files - -if __name__=="__main__": + +if __name__=="__main__": import subprocess from ivs.aux import loggers import shutil - import time logger = loggers.get_basic_logger() - to_install = ['spectra/pyrotin4', 'timeseries/deeming','timeseries/eebls','timeseries/multih', 'timeseries/pyclean','timeseries/pyKEP','timeseries/pyscargle', diff --git a/readme.txt b/readme.txt index 1a0ee1609..59a954d14 100644 --- a/readme.txt +++ b/readme.txt @@ -10,12 +10,30 @@ $ cd ivs $ git pull +* Make sure that your python path points to the ivs folder, so that you can simply import +using, for example: + + >>> from ivs.statistics import linearregression + + If your python repository is in, e.g., ~/python/ivs/ , you can put in your .bash_profile: + + export PYTHONPATH=/home/YOURNAME/python:$PYTHONPATH + + Warning: don't put ~/python/ivs in your Python path, but just ~/python. + +* For the next part of the installation, temporarily move the io folder elsewhere, e.g. + + $ cd ivs + $ mv io io_backup + +* In the config.py file in the ivs folder, add the path where the IvS data catalogs (variable: data_dir) can be found. + * The IvS Python repository contains mostly python routines. Some of the time-critical -functions, however, are written in fortran. To compile them you can run +functions, however, are written in fortran. This assumes you have gfortran installed on your system. To compile them you can run $ python config.py compile -If you want to specify your own fortran compiler you can do so with +Alternatively, if you want to specify your own fortran compiler you can do so with $ python config.py compile f77 @@ -26,21 +44,12 @@ Note: sometimes the compilation process fails. If so, try to compile spectra/pyr $ cd ../ $ python config.py compile -* In the config file you may also change the paths where the data catalogs (variable: data_dir) can be found, if you are not using the default locations. - -* Make sure that your python path points to the ivs folder, so that you can simply import -using, for example: - - >>> from ivs.statistics import linearregression - - If your python repository is in, e.g., ~/python/ivs/ , you can put in your .bash_profile: - - export PYTHONPATH=/home/YOURNAME/python:$PYTHONPATH - -Warning: don't put ~/python/ivs in your Python path, but just ~/python. +* With the installation finished, move the io folder back: + $ cd ivs + $ mv io_backup io -* To generate the documentation, simply run the script +* To generate the documentation (which, as a user, you typically do not have to do), simply run the script $ python makedoc.py @@ -66,3 +75,6 @@ Encountered errors and their solutions: numpy.distutils.fcompiler.CompilerNotFound: gnu95: f90 nor f77 A: Install gfortran. +2. Q: When I run "python config.py compile", I get the following error: +/bin/sh: f2py: command not found +A: Install f2py. If you do have f2py installed, but under a different name for the executable (for instance, f2py-2.7), replace f2py with f2py-2.7 in config.py (search for cmd = ). diff --git a/sed/builder.py b/sed/builder.py index b004dcab6..bcfee84b2 100644 --- a/sed/builder.py +++ b/sed/builder.py @@ -7,7 +7,7 @@ outside of the SED class if you know what you're doing. Table of contents: - + 1. Retrieving and plotting photometry of a target 2. Where/what is my target? 3. SED fitting using a grid based approach @@ -133,14 +133,14 @@ target, the script will B{not retrieve the photometry from the internet again}, but will use the contents of the file instead. The purpose is minimizing network traffic and maximizing speed. If you want to refresh the search, simply manually -delete the .phot file or set C{force=True} when calling L{get_photometry}. +delete the .phot file or set C{force=True} when calling L{get_photometry}. The content of the .phot file is most easily read using the L{ivs.io.ascii.read2recarray} function. Be careful, as it contains both absolute fluxes as flux ratios. >>> data = ascii.read2recarray('HD180642.phot') -Notice that in the C{.phot} files, also a C{comment} column is added. You can +Notice that in the C{.phot} files, also a C{comment} column is added. You can find translation of some of the flags here (i.e. upper limit, extended source etc..), or sometimes just additional remarks on variability etc. Not all catalogs have this feature implemented, so you are still responsible yourself for checking @@ -229,7 +229,7 @@ You can L{include}/L{exclude} photoemtry based on name, wavelength range, source and index, and only select absolute photometry or colors (L{include_abs},L{include_colors}). When working in interactive mode, in particular the index is useful. Print the -current set of photometry to the screen with +current set of photometry to the screen with >>> print(photometry2str(mysed.master,color=True,index=True)) @@ -280,9 +280,9 @@ as an argument so you can make e.g. the 'color' SED and 'absolute value' SED): 1. the grid (see L{SED.plot_grid}) - 2. the SED (see L{SED.plot_sed}) + 2. the SED (see L{SED.plot_sed}) 3. the fit statistics (see L{SED.plot_chi2}) - + To change the grid, load the L{ivs.sed.model} module and call L{ivs.sed.model.set_defaults} with appropriate arguments. See that module for conventions on the grid structure when adding custom grids. @@ -298,17 +298,17 @@ underestimated, it could be that the rescaling factor is also wrong, which means that the true probability region can be larger or smaller! -Subsection 3.2 Binary star +Subsection 3.2 Binary star -------------------------- -The SED class can create SEDs for multiple stars as well. There are 2 options +The SED class can create SEDs for multiple stars as well. There are 2 options available, the multiple SED fit which in theory can handle any number of stars, -and the binary SED fit which is for binaries only, and uses the mass of both +and the binary SED fit which is for binaries only, and uses the mass of both components to restrict the radii when combining two model SEDs. As an example we take the system PG1104+243, which consists of a subdwarf B star, and a G2 type mainsequence star. The photometry from the standard catalogues that -are build in this class, is of to low quality, so we use photometry obtained from +are build in this class, is of to low quality, so we use photometry obtained from U{the subdwarf database}. >>> mysed = SED('PG1104+243') @@ -367,86 +367,86 @@ the following information (we print part of each header):: EXTNAME = 'DATA ' - XTENSION= 'BINTABLE' / binary table extension - BITPIX = 8 / array data type - NAXIS = 2 / number of array dimensions - NAXIS1 = 270 / length of dimension 1 - NAXIS2 = 67 / length of dimension 2 - TTYPE1 = 'meas ' - TTYPE2 = 'e_meas ' - TTYPE3 = 'flag ' - TTYPE4 = 'unit ' - TTYPE5 = 'photband' - TTYPE6 = 'source ' - TTYPE7 = '_r ' - TTYPE8 = '_RAJ2000' - TTYPE9 = '_DEJ2000' - TTYPE10 = 'cwave ' - TTYPE11 = 'cmeas ' - TTYPE12 = 'e_cmeas ' - TTYPE13 = 'cunit ' - TTYPE14 = 'color ' - TTYPE15 = 'include ' - TTYPE16 = 'synflux ' - TTYPE17 = 'mod_eff_wave' - TTYPE18 = 'chi2 ' - - EXTNAME = 'MODEL ' - XTENSION= 'BINTABLE' / binary table extension - BITPIX = 8 / array data type - NAXIS = 2 / number of array dimensions - NAXIS1 = 24 / length of dimension 1 - NAXIS2 = 1221 / length of dimension 2 - TTYPE1 = 'wave ' - TUNIT1 = 'A ' - TTYPE2 = 'flux ' - TUNIT2 = 'erg/s/cm2/AA' - TTYPE3 = 'dered_flux' - TUNIT3 = 'erg/s/cm2/AA' - TEFFL = 23000.68377498454 - TEFF = 23644.49138963689 - TEFFU = 29999.33189058337 - LOGGL = 3.014445328877565 - LOGG = 4.984788506855546 - LOGGU = 4.996095055525759 - EBVL = 0.4703171900728142 - EBV = 0.4933185398871652 - EBVU = 0.5645144879211454 - ZL = -2.499929434601112 - Z = 0.4332336811329144 - ZU = 0.4999776537652627 - SCALEL = 2.028305798068863E-20 - SCALE = 2.444606991671813E-20 - SCALEU = 2.830281842143698E-20 - LABSL = 250.9613352757437 - LABS = 281.771013453664 - LABSU = 745.3149766772975 - CHISQL = 77.07958742733673 - CHISQ = 77.07958742733673 - CHISQU = 118.8587169011471 - CI_RAWL = 0.9999999513255379 - CI_RAW = 0.9999999513255379 - CI_RAWU = 0.9999999999999972 - CI_RED = 0.5401112973063139 - CI_REDL = 0.5401112973063139 - CI_REDU = 0.9500015229597392 - - EXTNAME = 'IGRID_SEARCH' - XTENSION= 'BINTABLE' / binary table extension - BITPIX = 8 / array data type - NAXIS = 2 / number of array dimensions - NAXIS1 = 80 / length of dimension 1 - NAXIS2 = 99996 / length of dimension 2 - TTYPE1 = 'teff ' - TTYPE2 = 'logg ' - TTYPE3 = 'ebv ' - TTYPE4 = 'z ' - TTYPE5 = 'chisq ' - TTYPE6 = 'scale ' - TTYPE7 = 'escale ' - TTYPE8 = 'Labs ' - TTYPE9 = 'CI_raw ' - TTYPE10 = 'CI_red ' + XTENSION= 'BINTABLE' / binary table extension + BITPIX = 8 / array data type + NAXIS = 2 / number of array dimensions + NAXIS1 = 270 / length of dimension 1 + NAXIS2 = 67 / length of dimension 2 + TTYPE1 = 'meas ' + TTYPE2 = 'e_meas ' + TTYPE3 = 'flag ' + TTYPE4 = 'unit ' + TTYPE5 = 'photband' + TTYPE6 = 'source ' + TTYPE7 = '_r ' + TTYPE8 = '_RAJ2000' + TTYPE9 = '_DEJ2000' + TTYPE10 = 'cwave ' + TTYPE11 = 'cmeas ' + TTYPE12 = 'e_cmeas ' + TTYPE13 = 'cunit ' + TTYPE14 = 'color ' + TTYPE15 = 'include ' + TTYPE16 = 'synflux ' + TTYPE17 = 'mod_eff_wave' + TTYPE18 = 'chi2 ' + + EXTNAME = 'MODEL ' + XTENSION= 'BINTABLE' / binary table extension + BITPIX = 8 / array data type + NAXIS = 2 / number of array dimensions + NAXIS1 = 24 / length of dimension 1 + NAXIS2 = 1221 / length of dimension 2 + TTYPE1 = 'wave ' + TUNIT1 = 'A ' + TTYPE2 = 'flux ' + TUNIT2 = 'erg/s/cm2/AA' + TTYPE3 = 'dered_flux' + TUNIT3 = 'erg/s/cm2/AA' + TEFFL = 23000.68377498454 + TEFF = 23644.49138963689 + TEFFU = 29999.33189058337 + LOGGL = 3.014445328877565 + LOGG = 4.984788506855546 + LOGGU = 4.996095055525759 + EBVL = 0.4703171900728142 + EBV = 0.4933185398871652 + EBVU = 0.5645144879211454 + ZL = -2.499929434601112 + Z = 0.4332336811329144 + ZU = 0.4999776537652627 + SCALEL = 2.028305798068863E-20 + SCALE = 2.444606991671813E-20 + SCALEU = 2.830281842143698E-20 + LABSL = 250.9613352757437 + LABS = 281.771013453664 + LABSU = 745.3149766772975 + CHISQL = 77.07958742733673 + CHISQ = 77.07958742733673 + CHISQU = 118.8587169011471 + CI_RAWL = 0.9999999513255379 + CI_RAW = 0.9999999513255379 + CI_RAWU = 0.9999999999999972 + CI_RED = 0.5401112973063139 + CI_REDL = 0.5401112973063139 + CI_REDU = 0.9500015229597392 + + EXTNAME = 'IGRID_SEARCH' + XTENSION= 'BINTABLE' / binary table extension + BITPIX = 8 / array data type + NAXIS = 2 / number of array dimensions + NAXIS1 = 80 / length of dimension 1 + NAXIS2 = 99996 / length of dimension 2 + TTYPE1 = 'teff ' + TTYPE2 = 'logg ' + TTYPE3 = 'ebv ' + TTYPE4 = 'z ' + TTYPE5 = 'chisq ' + TTYPE6 = 'scale ' + TTYPE7 = 'escale ' + TTYPE8 = 'Labs ' + TTYPE9 = 'CI_raw ' + TTYPE10 = 'CI_red ' Subsection 3.4 Loading SED fits @@ -495,7 +495,7 @@ Where C{radius} and C{distance} have equal units. The true absolute luminosity (solar units) is related to the absolute luminosity from the SED (solar units) via the radius (solar units):: - + >> L_abs_true = L_abs * radius**2 Finally, the angular diameter can be computed via:: @@ -509,7 +509,7 @@ independent constraint on the surface gravity, given the effective temperature of the model (you are free to give errors or not, just keep in mind that in the former case, you will get an array of Uncertainties rather than floats):: - + >> teff = 4260.,'K' >> nu_max = 38.90,0.86,'muHz' >> logg_slo = conversions.derive_logg_slo(teff,nu_max,unit='[cm/s2']) @@ -519,7 +519,7 @@ >> Deltanu0 = 4.80,0.02,'muHz' >> R_slo = conversions.derive_radius_slo(numax,Deltanu0,teff,unit='Rsol') - + Then from the radius, you can get both the absolute luminosity and distance to your star. @@ -530,7 +530,7 @@ dependent on some prior assumptions such as the shape of the galaxy and the distribution of stars. To estimate the probability density value due to a measured parallax of a star at particular distance, you can call L{distance.distprob}:: - + >> gal_latitude = 0.5 >> plx = 3.14,0.5 >> d_prob = distance.distprob(d,gal_latitude,plx) @@ -565,7 +565,10 @@ import pylab as pl from matplotlib import mlab -import Image +try: + import image +except ImportError: + import Image as image import numpy as np import scipy.stats from scipy.interpolate import Rbf @@ -607,9 +610,9 @@ def fix_master(master,e_default=None): """ Clean/extend/fix record array received from C{get_photometry}. - + This function does a couple of things: - + 1. Adds common but uncatalogized colors like 2MASS.J-H if not already present. WARNING: these colors can mix values from the same system but from different catalogs! @@ -624,7 +627,7 @@ def fix_master(master,e_default=None): value are untrostworthy because the calibration error is larger than that. 7. USNOB1 and ANS photometry are set the have a minimum error of 30% due to uncertainties in response curves. - + @param master: record array containing photometry. This should have the fields 'meas','e_meas','unit','photband','source','_r','_RAJ2000','DEJ2000', 'cmeas','e_cmeas','cwave','cunit' @@ -648,7 +651,7 @@ def fix_master(master,e_default=None): add_rows = [] for cat in set_cats: master_ = master[cats==cat] - + for color in ['2MASS.J-H','2MASS.KS-H','TD1.1565-1965','TD1.2365-1965','TD1.2365-2740', 'JOHNSON.J-H','JOHNSON.K-H','JOHNSON.B-V','JOHNSON.U-B','JOHNSON.R-V','JOHNSON.I-V', 'GALEX.FUV-NUV','TYCHO2.BT-VT','WISE.W1-W2','WISE.W3-W2','WISE.W4-W3', @@ -660,17 +663,17 @@ def fix_master(master,e_default=None): system,band = color.split('.') band0,band1 = band.split('-') band0,band1 = '%s.%s'%(system,band0),'%s.%s'%(system,band1) - + if band0 in master_['photband'] and band1 in master_['photband'] and not color in master_['photband']: #-- where are the bands located? index0 = list(master_['photband']).index(band0) index1 = list(master_['photband']).index(band1) - + #-- start a new row to add the color row = list(master_[index1]) row[columns.index('photband')] = color row[columns.index('cwave')] = np.nan - + #-- it could be a magnitude difference if master_['unit'][index0]=='mag': row[columns.index('meas')] = master_['meas'][index0]-master_['meas'][index1] @@ -689,7 +692,7 @@ def fix_master(master,e_default=None): row[columns.index('cunit')] = 'flux_ratio' add_rows.append(tuple(row)) master = numpy_ext.recarr_addrows(master,add_rows) - + #-- add an extra column with a flag to distinguish colors from absolute # fluxes, and a column with flags to include/exclude photometry # By default, exclude photometry if the effective wavelength is above @@ -709,17 +712,17 @@ def fix_master(master,e_default=None): master['e_cmeas'][no_error] = np.abs(e_default*master['cmeas'][no_error]) small_error = master['e_cmeas']<(0.01*np.abs(master['cmeas'])) master['e_cmeas'][small_error] = 0.01*np.abs(master['cmeas'][small_error]) - + #-- the measurements from USNOB1 are not that good because the response # curve is approximated by the JOHNSON filter. Set the errors to min 30% # The same holds for ANS: here, the transmission curves are very uncertain for i,photband in enumerate(master['photband']): if 'USNOB1' in photband or 'ANS' in photband: master['e_cmeas'][i] = max(0.30*master['cmeas'][i],master['e_cmeas'][i]) - + #-- remove negative fluxes master = master[master['cmeas']>0] - + return master @@ -727,43 +730,43 @@ def decide_phot(master,names=None,wrange=None,sources=None,indices=None,ptype='a """ Exclude/include photometric passbands containing one of the strings listed in photbands. - + This function will set the flags in the 'include' field of the extended master record array. - + Colours also have a wavelength in this function, and it is defined as the average wavelength of the included bands (doesn't work for stromgren C1 and M1). - + include=False excludes photometry based on the input parameters. include=True includes photometry based on the input parameters. - + Some examples: - + 1. Exclude all measurements:: - + >> decide_phot(master,wrange=(-np.inf,+np.inf),ptype='all',include=False) - + 2. Include all TD1 fluxes and colors:: - + >> decide_phot(master,names=['TD1'],ptype='all',include=True) - + 3. Include all V band measurements from all systems (but not the colors):: - + >> decide_phot(master,names=['.V'],ptype='abs',include=True) - + 4. Include all Geneva colors and exclude Geneva magnitudes:: - + >> decide_phot(master,names=['GENEVA'],ptype='col',include=True) >> decide_phot(master,names=['GENEVA'],ptype='abs',include=False) - + 5. Exclude all infrared measurements beyond 1 micron:: - + >> decide_phot(master,wrange=(1e4,np.inf),ptype='all',include=False) - + 6. Include all AKARI measurements below 10 micron:: - + >> decide_phot(master,names=['AKARI'],wrange=(-np.inf,1e5),ptype='all',include=True) - + @param master: record array containing all photometry @type master: numpy record array @param names: strings excerpts to match filters @@ -781,7 +784,7 @@ def decide_phot(master,names=None,wrange=None,sources=None,indices=None,ptype='a @type include: boolean @return: master record array with adapted 'include' flags @rtype: numpy record array - """ + """ #-- exclude/include passbands based on their names if names is not None: logger.info('%s photometry based on photband containining one of %s'%((include and 'Include' or "Exclude"),names)) @@ -805,7 +808,7 @@ def decide_phot(master,names=None,wrange=None,sources=None,indices=None,ptype='a master['include'][index] = include #-- exclude/include passbands based on their wavelength if not ptype=='col': - master['include'][(wrange[0]>> mysed = SED('HD129929') - + B{2. If previous data exists:} - + >>> #mysed = SED(photfile='HD129929.phot') # will set ID with 'oname' field from SIMBAD >>> #mysed = SED(ID='bla', photfile='HD129929.phot') # Sets custom ID - + The C{ID} variable is used internally to look up data, so it should be something SIMBAD understands and that designates the target. - + @param plx: parallax (and error) of the object @type plx: tuple (plx,e_plx) """ @@ -1002,7 +1005,7 @@ def __init__(self,ID=None,photfile=None,plx=None,load_fits=True,load_hdf5=True,l # coordinates, which are not available from SESAME. else: self.photfile = photfile - + #-- load information from the photometry file if it exists if not os.path.isfile(self.photfile): try: @@ -1020,11 +1023,11 @@ def __init__(self,ID=None,photfile=None,plx=None,load_fits=True,load_hdf5=True,l try: self.info['jradeg'],self.info['jdedeg'] = ra,dec = corot.resolve(corot_id) gal = conversions.convert('equatorial','galactic',(ra,dec),epoch='2000') - self.info['galpos'] = float(gal[0])/np.pi*180,float(gal[1])/np.pi*180 + self.info['galpos'] = float(gal[0])/np.pi*180,float(gal[1])/np.pi*180 logger.info('Resolved position via CoRoT EXO catalog') except: logger.warning("Star %s not recognised by CoRoT"%(os.path.basename(ID))) - + if plx is not None: if not 'plx' in self.info: self.info['plx'] = {} @@ -1042,16 +1045,16 @@ def __init__(self,ID=None,photfile=None,plx=None,load_fits=True,load_hdf5=True,l self.load_fits() if load_hdf5: self.load_hdf5() - + #-- prepare for information on fitting processes self.CI_limit = 0.95 - + def __repr__(self): """ Machine readable string representation of an SED object. """ return "SED('{}')".format(self.ID) - + def __str__(self): """ Human readable string representation of an SED object. @@ -1079,16 +1082,16 @@ def __str__(self): if hasattr(self,'master'): txt.append(photometry2str(self.master,color=True)) return "\n".join(txt) - + #{ Handling photometric data def get_photometry(self,radius=None,ra=None,dec=None, include=None,exclude=None,force=False, units='erg/s/cm2/AA'): """ Search photometry on the net or from the phot file if it exists. - + For bright stars, you can set radius a bit higher... - + @param radius: search radius (arcseconds) @type radius: float. """ @@ -1111,29 +1114,29 @@ def get_photometry(self,radius=None,ra=None,dec=None, if 'jradeg' in self.info: master['_RAJ2000'] -= self.info['jradeg'] master['_DEJ2000'] -= self.info['jdedeg'] - + #-- fix the photometry: set default errors to 2% and print it to the # screen self.master = fix_master(master,e_default=0.1) logger.info('\n'+photometry2str(master)) - + #-- write to file self.save_photometry() - + def get_spectrophotometry(self,directory=None,force_download=False): """ Retrieve and combine spectra. - + B{WARNING:} this function creates FUSE and DIR directories! """ if directory is None and os.path.dirname(self.ID) == '': directory = os.getcwd() elif os.path.dirname(self.ID) != '': directory = os.path.dirname(self.ID) - + if not os.path.isdir(directory): os.mkdir(directory) - + #-- add spectrophotometric filters to the set photbands = filters.add_spectrophotometric_filters(R=200,lambda0=950,lambdan=3350) if hasattr(self,'master') and self.master is not None and not force_download: @@ -1148,14 +1151,14 @@ def get_spectrophotometry(self,directory=None,force_download=False): out1 = [] else: out1 = glob.glob(fuse_direc+'/*') - #-- IUE spectra + #-- IUE spectra if not os.path.isdir(iue_direc) or force_download: out2 = vizier.get_IUE_spectra(ID=os.path.basename(self.ID),directory=iue_direc,select='lo') if out2 is None: out2 = [] else: out2 = glob.glob(iue_direc+'/*') - #-- read them in to combine + #-- read them in to combine list_of_spectra = [fits.read_fuse(ff)[:3] for ff in out1] list_of_spectra+= [fits.read_iue(ff)[:3] for ff in out2] #-- combine @@ -1165,25 +1168,25 @@ def get_spectrophotometry(self,directory=None,force_download=False): units = ['erg/s/cm2/AA']*N source = ['specphot']*N self.add_photometry_fromarrays(flux,err,units,photbands,source,flags=nspec) - + #-- write to file self.master = fix_master(self.master,e_default=0.1) logger.info('\n'+photometry2str(self.master)) self.save_photometry() - - - - + + + + def exclude(self,names=None,wrange=None,sources=None,indices=None): """ Exclude (any) photometry from fitting process. - + If called without arguments, all photometry will be excluded. """ if names is None and wrange is None and sources is None and indices is None: wrange = (-np.inf,np.inf) decide_phot(self.master,names=names,wrange=wrange,sources=sources,indices=indices,include=False,ptype='all') - + def exclude_colors(self,names=None,wrange=None,sources=None,indices=None): """ Exclude (color) photometry from fitting process. @@ -1191,14 +1194,14 @@ def exclude_colors(self,names=None,wrange=None,sources=None,indices=None): if names is None and wrange is None and sources is None and indices is None: wrange = (-np.inf,0) decide_phot(self.master,names=names,wrange=wrange,sources=sources,indices=indices,include=False,ptype='col') - + def exclude_abs(self,names=None,wrange=None,sources=None,indices=None): """ Exclude (absolute) photometry from fitting process. """ decide_phot(self.master,names=names,wrange=wrange,sources=sources,indices=indices,include=False,ptype='abs') - - + + def include(self,names=None,wrange=None,sources=None,indices=None): """ Include (any) photometry in fitting process. @@ -1206,32 +1209,32 @@ def include(self,names=None,wrange=None,sources=None,indices=None): if names is None and wrange is None and sources is None and indices is None: wrange = (-np.inf,np.inf) decide_phot(self.master,names=names,wrange=wrange,sources=sources,indices=indices,include=True,ptype='all') - + def include_colors(self,names=None,wrange=None,sources=None,indices=None): """ Include (color) photometry in fitting process. """ decide_phot(self.master,names=names,wrange=wrange,sources=sources,indices=indices,include=True,ptype='col') - + def include_abs(self,names=None,wrange=None,sources=None,indices=None): """ Include (absolute) photometry in fitting process. """ decide_phot(self.master,names=names,wrange=wrange,sources=sources,indices=indices,include=True,ptype='abs') - + def set_photometry_scheme(self,scheme,infrared=(1,'micron')): """ Set a default scheme of colors/absolute values to fit the SED. - + Possible values: - + 1. scheme = 'abs': means excluding all colors, including all absolute values 2. scheme = 'color': means including all colors, excluding all absolute values 3. scheme = 'combo': means inculding all colors, and one absolute value per system (the one with the smallest relative error) 4. scheme = 'irfm': means mimic infrared flux method: choose absolute values in the infrared (define with C{infrared}), and colours in the optical - + @param infrared: definition of start of infrared for infrared flux method @type infrared: tuple (value , unit ) """ @@ -1284,15 +1287,15 @@ def set_photometry_scheme(self,scheme,infrared=(1,'micron')): use_colors = sum(self.master['include'] & self.master['color']) use_abs = sum(self.master['include'] & -self.master['color']) logger.info('Fitting procedure will use colors (%d) < %.3g%s < absolute values (%d)'%(use_colors,infrared[0],infrared[1],use_abs)) - - + + def add_photometry_fromarrays(self,meas,e_meas,units,photbands,source,flags=None): """ Add photometry from numpy arrays - + By default unflagged, and at the ra and dec of the star. No color and included. - + @param meas: original measurements (fluxes, magnitudes...) @type meas: array @param e_meas: error on original measurements in same units as C{meas} @@ -1347,7 +1350,7 @@ def add_photometry_fromarrays(self,meas,e_meas,units,photbands,source,flags=None extra_master['bibcode'][i] = '-' if 'comments' in extra_master.dtype.names: extra_master['comments'][i] = '-' - + logger.info('Original measurements:\n%s'%(photometry2str(self.master))) logger.info('Appending:\n%s'%(photometry2str(extra_master))) self.master = fix_master(np.hstack([self.master,extra_master]),e_default=0.1) @@ -1356,14 +1359,14 @@ def add_photometry_fromarrays(self,meas,e_meas,units,photbands,source,flags=None #} #{ Additional information - + def is_target(self,name): """ Check if this SED represents the object `name'. - + Purpose: solve alias problems. Maybe the ID is 'HD129929', and you are checking for "V* V836 Cen", which is the same target. - + @param name: object name @type name: str @return: True if this object instance represent the target "name". @@ -1377,34 +1380,34 @@ def is_target(self,name): return False if oname==self.info['oname']: return True - + def has_photfile(self): """ Check if this SED has a phot file. - + @return: True if this object instance has a photfile @rtype: bool """ return os.path.isfile(self.photfile) - + def get_distance_from_plx(self,plx=None,lutz_kelker=True,unit='pc'): """ Get the probability density function for the distance, given a parallax. - + If no parallax is given, the catalogue value will be used. If that one is also not available, a ValueError will be raised. - + Parallax must be given in mas. - + If the parallax is a float without uncertainty, an uncertainty of 0 will be assumed. - + If C{lutz_kelker=True}, a probability density function will be given for the distance, otherwise the distance is computed from simply inverting the parallax. - + Distance is returned in parsec (pc). - + @return: distance @rtype: (float,float) """ @@ -1429,7 +1432,7 @@ def get_distance_from_plx(self,plx=None,lutz_kelker=True,unit='pc'): dist = (conversions.Unit(plx,'kpc-1')**-1).convert(unit) logger.info('Computed distance to target via parallax inversion: {:s}'.format(dist)) return dist.get_value() - + def get_interstellar_reddening(self,distance=None, Rv=3.1): """ Under construction. @@ -1443,24 +1446,24 @@ def get_interstellar_reddening(self,distance=None, Rv=3.1): if ext is not None: output[model] = ext/Rv return output - + def get_angular_diameter(self): """ Under construction. """ raise NotImplementedError - + def compute_distance(self,mtype='igrid_search'): """ Compute distance from radius and angular diameter (scale factor). - + The outcome is added to the C{results} attribute:: - + distance = r/sqrt(scale) - + This particularly useful when you added constraints from solar-like oscillations (L{add_constraint_slo}). - + @return: distance,uncertainty (pc) @rtype: (float,float) """ @@ -1469,37 +1472,37 @@ def compute_distance(self,mtype='igrid_search'): e_radius = grid['e_radius'] scale = grid['scale'] e_scale = grid['escale'] - + radius = conversions.unumpy.uarray([radius,e_radius]) scale = conversions.unumpy.uarray([scale,e_scale]) - + distance = radius/scale**0.5 distance = conversions.convert('Rsol','pc',distance) distance,e_distance = conversions.unumpy.nominal_values(distance),\ conversions.unumpy.std_devs(distance) self.results[mtype]['grid'] = pl.mlab.rec_append_fields(grid,\ ['distance','e_distance'],[distance,e_distance]) - d,e_d = distance.mean(),distance.std() + d,e_d = distance.mean(),distance.std() logger.info("Computed distance from R and scale: {0}+/-{1} pc".format(d,e_d)) return d,e_d - - - - + + + + #} - + def generate_ranges(self,type='single',start_from='igrid_search',distribution='uniform',**kwargs): """ Generate sensible search range for each parameter. """ limits = {} exist_previous = (start_from in self.results and 'CI' in self.results[start_from]) - + #-- run over all parnames, and generate the parameters ranges for each # component: # (1) if a previous parameter search is done and the user has not set # the parameters range, derive a sensible parameter range from the - # previous results. + # previous results. # (2) if no previous search is done, use infinite upper and lower bounds # (3) if ranges are given, use those. Cycle over them, and use the # strategy from (1) if no parameter ranges are given for a specific @@ -1520,7 +1523,7 @@ def generate_ranges(self,type='single',start_from='igrid_search',distribution='u #-- else we can derive a better parameter range else: parrange = (self.results[start_from]['CI'][lkey], - self.results[start_from]['CI'][ukey]) + self.results[start_from]['CI'][ukey]) elif parrange is None: parrange = (-np.inf,np.inf) @@ -1535,13 +1538,13 @@ def generate_ranges(self,type='single',start_from='igrid_search',distribution='u # the type if it was given, or gives the type when it needed to be derived logger.info('Parameter ranges calculated starting from {0:s} and using distribution {1:s}.'.format(start_from,distribution)) return limits - + #{ Fitting routines - + def clip_grid(self,mtype='igrid_search',CI_limit=None): """ Clip grid on CI limit, to save memory. - + @param mtype: type or results to clip @type mtype: str @param CI_limit: confidence limit to clip on @@ -1553,27 +1556,27 @@ def clip_grid(self,mtype='igrid_search',CI_limit=None): new_grid = new_grid[new_grid['ci_red']<=CI_limit] self.results[mtype]['grid'] = new_grid logger.info("Clipped grid at {:.6f}%".format(CI_limit*100)) - + def collect_results(self, grid=None, fitresults=None, mtype='igrid_search', selfact='chisq', **kwargs): """creates a record array of all fit results and removes the failurs""" - + gridnames = sorted(grid.keys()) fitnames = sorted(fitresults.keys()) pardtypes = [(name,'f8') for name in gridnames]+[(name,'f8') for name in fitnames] pararrays = [grid[name] for name in gridnames]+[fitresults[name] for name in fitnames] grid_results = np.rec.fromarrays(pararrays,dtype=pardtypes) - + #-- exclude failures failures = np.isnan(grid_results[selfact]) if sum(failures): logger.info('Excluded {0} failed results (nan)'.format(sum(failures))) grid_results = grid_results[-failures] - + #-- make room for chi2 statistics grid_results = mlab.rec_append_fields(grid_results, 'ci_raw', np.zeros(len(grid_results))) grid_results = mlab.rec_append_fields(grid_results, 'ci_red', np.zeros(len(grid_results))) - + #-- take the previous results into account if they exist: if not mtype in self.results: self.results[mtype] = {} @@ -1583,19 +1586,19 @@ def collect_results(self, grid=None, fitresults=None, mtype='igrid_search', self ex_grid = np.rec.fromarrays([self.results[mtype]['grid'][exname] for exname in ex_names], names=ex_names) grid_results = np.hstack([ex_grid,grid_results]) - - #-- inverse sort according to chisq: this means the best models are at the end + + #-- inverse sort according to chisq: this means the best models are at the end # (mainly for plotting reasons, so that the best models are on top). sa = np.argsort(grid_results[selfact])[::-1] grid_results = grid_results[sa] - + self.results[mtype]['grid'] = grid_results - + logger.info('Total of %d grid points, best chisq=%s'%(len(grid_results), grid_results[-1])) - + def calculateDF(self, **ranges): """ Calculates the degrees of freedom from the given ranges""" - + df, df_info = 1, ['theta'] for range_name in ranges: if re.search('ebv\d?range$', range_name): @@ -1609,22 +1612,22 @@ def calculateDF(self, **ranges): df_info.append(range_name[0:-5]) df_info.sort() logger.info('Degrees of freedom = {} ({})'.format(df,', '.join(df_info))) - - return df, df_info - + + return df, df_info + def calculate_statistics(self, df=None, ranges=None, mtype='igrid_search', selfact='chisq'): """ Calculates the Chi2 and reduced Chi2 based on nr of observations and degrees of freedom. If df is not provided, tries to calculate df from the fitting ranges. """ - + #-- If nessessary calculate degrees of freedom from the ranges if df == None and ranges != None: df,df_info = self.calculateDF(**ranges) elif df == None: logger.warning('Cannot compute degrees of freedom!!! CHI2 might not make sense. (using df=5)') df = 5 - + #-- Do the statistics N = sum(self.master['include']) k = N-df @@ -1632,22 +1635,22 @@ def calculate_statistics(self, df=None, ranges=None, mtype='igrid_search', selfa logger.warning('Not enough data to compute CHI2: it will not make sense') k = 1 logger.info('Calculated statistics based on df={0} and Nobs={1}'.format(df,N)) - + #-- Rescale if needed and compute confidence intervals results = self.results[mtype]['grid'] factor = max(results[selfact][-1]/k,1) logger.warning('CHI2 rescaling factor equals %g'%(factor)) results['ci_raw'] = scipy.stats.distributions.chi2.cdf(results[selfact],k) results['ci_red'] = scipy.stats.distributions.chi2.cdf(results[selfact]/factor,k) - + #-- Store the results self.results[mtype]['grid'] = results self.results[mtype]['factor'] = factor - + def calculate_confidence_intervals(self,mtype='igrid_search',chi2_type='red',CI_limit=None): """ Compute confidence interval of all columns in the results grid. - + @param mtype: type of results to compute confidence intervals of @type mtype: str @param chi2_type: type of chi2 (raw or reduced) @@ -1669,17 +1672,17 @@ def calculate_confidence_intervals(self,mtype='igrid_search',chi2_type='red',CI_ cilow.append(grid_results[name][region].min()) value.append(grid_results[name][-1]) cihigh.append(grid_results[name][region].max()) - + return dict(name=grid_results.dtype.names, value=value, cilow=cilow, cihigh=cihigh) - + def store_confidence_intervals(self, mtype='igrid_search', name=None, value=None, cilow=None, \ cihigh=None, **kwargs): """ Saves the provided confidence intervals in the result dictionary of self. - The provided confidence intervals will be saved in a dictionary: + The provided confidence intervals will be saved in a dictionary: self.results[mtype]['CI'] with as key for the value the name, for cilow: name_l and for cihigh: name_u. - + @param mtype: the search type @type mtype: str @param name: names of the parameters @@ -1697,57 +1700,57 @@ def store_confidence_intervals(self, mtype='igrid_search', name=None, value=None self.results[mtype]['CI'][name+'_l'] = cil self.results[mtype]['CI'][name] = val self.results[mtype]['CI'][name+'_u'] = cih - + # Send the stored CI to logger logger.info(self.ci2str(mtype=mtype)) - + def igrid_search(self,points=100000,teffrange=None,loggrange=None,ebvrange=None, zrange=(0,0),rvrange=(3.1,3.1),vradrange=(0,0), df=None,CI_limit=None,set_model=True,**kwargs): """ Fit fundamental parameters using a (pre-integrated) grid search. - + If called consecutively, the ranges will be set to the CI_limit of previous estimations, unless set explicitly. - + If called for the first time, the ranges will be +/- np.inf by defaults, unless set explicitly. """ if CI_limit is None or CI_limit > 1.0: CI_limit = self.CI_limit - + #-- set defaults limits ranges = self.generate_ranges(teffrange=teffrange,\ loggrange=loggrange,ebvrange=ebvrange,\ zrange=zrange,rvrange=rvrange,vradrange=vradrange) - + #-- grid search on all include data: extract the best CHI2 include_grid = self.master['include'] logger.info('The following measurements are included in the fitting process:\n%s'%(photometry2str(self.master[include_grid]))) - + #-- build the grid, run over the grid and calculate the CHI2 - pars = fit.generate_grid_pix(self.master['photband'][include_grid],points=points,**ranges) + pars = fit.generate_grid_pix(self.master['photband'][include_grid],points=points,**ranges) chisqs,scales,e_scales,lumis = fit.igrid_search_pix(self.master['cmeas'][include_grid], self.master['e_cmeas'][include_grid], self.master['photband'][include_grid],**pars) fitres = dict(chisq=chisqs, scale=scales, escale=e_scales, labs=lumis) - + #-- collect all the results in a record array self.collect_results(grid=pars, fitresults=fitres, mtype='igrid_search') - + #-- do the statistics self.calculate_statistics(df=df, ranges=ranges, mtype='igrid_search') - + #-- compute the confidence intervals ci = self.calculate_confidence_intervals(mtype='igrid_search',chi2_type='red',CI_limit=CI_limit) self.store_confidence_intervals(mtype='igrid_search', **ci) - + #-- remember the best model if set_model: self.set_best_model() - + def generate_fit_param(self, start_from='igrid_search', **pars): - """ - generates a dictionary with parameter information that can be handled by fit.iminimize + """ + generates a dictionary with parameter information that can be handled by fit.iminimize """ result = dict() for key in pars.keys(): @@ -1766,11 +1769,11 @@ def generate_fit_param(self, start_from='igrid_search', **pars): result[key+"_value"] = pars[key] elif pars[key] == None and not key+"_value" in result: result[key+"_value"] = self.results[start_from]['CI'][key] - - return result - + + return result + def calculate_iminimize_CI(self, mtype='iminimize', CI_limit=0.66, **kwargs): - + #-- Get the best fit parameters and ranges pars = {} skip = ['scale', 'chisq', 'nfev', 'labs', 'ci_raw', 'ci_red', 'scale', 'escale'] @@ -1781,27 +1784,27 @@ def calculate_iminimize_CI(self, mtype='iminimize', CI_limit=0.66, **kwargs): pars[name+"range"] = [self.results[mtype]['CI'][name+"_l"],\ self.results[mtype]['CI'][name+"_u"]] pars.update(kwargs) - pars = self.generate_fit_param(**pars) - + pars = self.generate_fit_param(**pars) + #-- calculate the confidence intervalls include_grid = self.master['include'] ci = fit.calculate_iminimize_CI(self.master['cmeas'][include_grid], self.master['e_cmeas'][include_grid], self.master['photband'][include_grid], CI_limit=CI_limit,constraints=self.constraints, **pars) - + #-- Add the scale factor ci['name'] = np.append(ci['name'], 'scale') ci['value'] = np.append(ci['value'], self.results[mtype]['grid']['scale'][-1]) ci['cilow'] = np.append(ci['cilow'], min(self.results[mtype]['grid']['scale'])) ci['cihigh'] = np.append(ci['cihigh'], max(self.results[mtype]['grid']['scale'])) - + logger.info('Calculated %s%% confidence intervalls for all parameters'%(CI_limit)) - + self.store_confidence_intervals(mtype='iminimize', **ci) - + def calculate_iminimize_CI2D(self,xpar, ypar, mtype='iminimize', limits=None, res=10, **kwargs): - + #-- get the best fitting parameters pars = {} skip = ['scale', 'chisq', 'nfev', 'labs', 'ci_raw', 'ci_red', 'scale', 'escale'] @@ -1812,24 +1815,24 @@ def calculate_iminimize_CI2D(self,xpar, ypar, mtype='iminimize', limits=None, re pars[name+"range"] = [self.results[mtype]['CI'][name+"_l"],\ self.results[mtype]['CI'][name+"_u"]] pars.update(kwargs) - pars = self.generate_fit_param(**pars) - + pars = self.generate_fit_param(**pars) + logger.info('Calculating 2D confidence intervalls for %s--%s'%(xpar,ypar)) - + #-- calculate the confidence intervalls include_grid = self.master['include'] x,y,chi2, raw, red = fit.calculate_iminimize_CI2D(self.master['cmeas'][include_grid], self.master['e_cmeas'][include_grid], self.master['photband'][include_grid], - xpar, ypar, limits=limits, res=res, + xpar, ypar, limits=limits, res=res, constraints=self.constraints, **pars) - + #-- store the CI if not 'CI2D' in self.results[mtype]: self.results[mtype]['CI2D'] = {} self.results[mtype]['CI2D'][xpar+"-"+ypar] = dict(x=x, y=y, ci_chi2=chi2,\ ci_raw=raw, ci_red=red) - - + + def _get_imin_ci(self, mtype='iminimize',**ranges): """ returns ci information for store_confidence_intervals """ names, values, cil, cih = [],[],[],[] @@ -1843,16 +1846,16 @@ def _get_imin_ci(self, mtype='iminimize',**ranges): values.append(self.results[mtype]['grid']['scale'][-1]) cil.append(min(self.results[mtype]['grid']['scale'])) cih.append(max(self.results[mtype]['grid']['scale'])) - return dict(name=names, value=values, cilow=cil, cihigh=cih) - + return dict(name=names, value=values, cilow=cil, cihigh=cih) + def iminimize(self, teff=None, logg=None, ebv=None, z=0, rv=3.1, vrad=0, teffrange=None, loggrange=None, ebvrange=None, zrange=None, rvrange=None, vradrange=None, - points=None, distance=None, start_from='igrid_search',df=None, CI_limit=None, + points=None, distance=None, start_from='igrid_search',df=None, CI_limit=None, calc_ci=False, set_model=True, **kwargs): - """ + """ Basic minimizer method for SED fitting implemented using the lmfit library from sigproc.fit """ - + #-- set defaults limits and starting values ranges = self.generate_ranges(teffrange=teffrange,loggrange=loggrange,\ ebvrange=ebvrange,zrange=zrange,rvrange=rvrange,\ @@ -1860,43 +1863,43 @@ def iminimize(self, teff=None, logg=None, ebv=None, z=0, rv=3.1, vrad=0, teffran pars = self.generate_fit_param(teff=teff, logg=logg, ebv=ebv, z=z, \ rv=rv, vrad=vrad, start_from=start_from, **ranges) fitkws = dict(distance=distance) if distance != None else dict() - + #-- grid search on all include data: extract the best CHI2 include_grid = self.master['include'] logger.info('The following measurements are included in the fitting process:\n%s'% \ (photometry2str(self.master[include_grid]))) - + #-- pass all ranges and starting values to the fitter grid, chisq, nfev, scale, lumis = fit.iminimize(self.master['cmeas'][include_grid], self.master['e_cmeas'][include_grid], self.master['photband'][include_grid], fitkws=fitkws, points=points,**pars) - + logger.info('Minimizer Succes with startpoints=%s, chi2=%s, nfev=%s'%(len(chisq), chisq[0], nfev[0])) #-- handle the results fitres = dict(chisq=chisq, nfev=nfev, scale=scale, labs=lumis) self.collect_results(grid=grid, fitresults=fitres, mtype='iminimize', selfact='chisq') - + #-- Do the statistics self.calculate_statistics(df=5, ranges=ranges, mtype='iminimize', selfact='chisq') - + #-- Store the confidence intervals ci = self._get_imin_ci(mtype='iminimize',**ranges) self.store_confidence_intervals(mtype='iminimize', **ci) if calc_ci: self.calculate_iminimize_CI(mtype='iminimize', CI_limit=CI_limit) - + #-- remember the best model if set_model: self.set_best_model(mtype='iminimize') - - + + def imc(self,teffrange=None,loggrange=None,ebvrange=None,zrange=None,start_from='imc',\ distribution='uniform',points=None,fitmethod='fmin',disturb=True): - + limits,type = self.generate_ranges(teffrange=teffrange,loggrange=loggrange,\ ebvrange=ebvrange,zrange=zrange,distribution=distribution,\ start_from='imc') - + #-- grid search on all include data: extract the best CHI2 include = self.master['include'] meas = self.master['cmeas'][include] @@ -1906,13 +1909,13 @@ def imc(self,teffrange=None,loggrange=None,ebvrange=None,zrange=None,start_from= emeas = self.master['e_cmeas'][include]*1e-6 photbands = self.master['photband'][include] logger.info('The following measurements are included in the fitting process:\n%s'%(photometry2str(self.master[include]))) - + #-- generate initial guesses - teffs,loggs,ebvs,zs,radii = fit.generate_grid(self.master['photband'][include],type=type,points=points+25,**limits) + teffs,loggs,ebvs,zs,radii = fit.generate_grid(self.master['photband'][include],type=type,points=points+25,**limits) NrPoints = len(teffs)>points and points or len(teffs) firstoutput = np.zeros((len(teffs)-NrPoints,9)) output = np.zeros((NrPoints,9)) - + #-- fit the original data a number of times for i,(teff,logg,ebv,z) in enumerate(zip(teffs[NrPoints:],loggs[NrPoints:],ebvs[NrPoints:],zs[NrPoints:])): try: @@ -1921,20 +1924,20 @@ def imc(self,teffrange=None,loggrange=None,ebvrange=None,zrange=None,start_from= firstoutput[i,0] = warnflag except IOError: firstoutput[i,0] = 3 - + logger.info("{0}/{1} fits on original data failed (max func call)".format(sum(firstoutput[:,0]==1),firstoutput.shape[0])) logger.info("{0}/{1} fits on original failed (max iter)".format(sum(firstoutput[:,0]==2),firstoutput.shape[0])) logger.info("{0}/{1} fits on original data failed (outside of grid)".format(sum(firstoutput[:,0]==3),firstoutput.shape[0])) - + #-- retrieve the best fitting result and make it the first entry of output keep = (firstoutput[:,0]==0) & (firstoutput[:,1]>0) best = firstoutput[keep,-2].argmin() output[-1,:] = firstoutput[keep][best,:] - + # calculate the factor with which to multiply the scale #factor = np.sqrt(output[-1,5]/len(meas)) #print factor - + #-- now do the actual Monte Carlo simulation for i,(teff,logg,ebv,z) in enumerate(zip(teffs[:NrPoints-1],loggs[:NrPoints-1],ebvs[:NrPoints-1],zs[:NrPoints-1])): newmeas = meas + np.random.normal(scale=emeas) #*factor) @@ -1944,18 +1947,18 @@ def imc(self,teffrange=None,loggrange=None,ebvrange=None,zrange=None,start_from= output[i,0] = warnflag except IOError: output[i,0] = 3 - + logger.info("{0}/{1} MC simulations failed (max func call)".format(sum(output[:,0]==1),NrPoints)) logger.info("{0}/{1} MC simulations failed (max iter)".format(sum(output[:,0]==2),NrPoints)) logger.info("{0}/{1} MC simulations failed (outside of grid)".format(sum(output[:,0]==3),NrPoints)) - + #-- remove nonsense results keep = (output[:,0]==0) & (output[:,1]>0) output = output[keep] output = np.rec.fromarrays(output[:,1:-1].T,names=['teff','logg','ebv','z','labs','chisq','scale']) #-- derive confidence intervals and median values #print np.median(output[:,1:],axis=0) - + self.results['imc'] = {} self.results['imc']['CI'] = {} self.results['imc']['grid'] = output @@ -1971,29 +1974,29 @@ def imc(self,teffrange=None,loggrange=None,ebvrange=None,zrange=None,start_from= self.results['imc']['CI'][name+'_u'])) self.set_best_model(mtype='imc') - + #} - + #{ Interfaces - + def clear(self): """ Clear the results. """ self.results = {} - + def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004',**kwargs): """ Get reddenend and unreddened model """ logger.info('Interpolating approximate full SED of best model') - + #-- synthetic flux include = self.master['include'] synflux = np.zeros(len(self.master['photband'])) keep = (self.master['cwave']<1.6e6) | np.isnan(self.master['cwave']) keep = keep & include - + if mtype in ['igrid_search', 'iminimize']: #-- get the metallicity right files = model.get_file(z='*') @@ -2018,11 +2021,11 @@ def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004',**kwargs): ebv=self.results[mtype]['CI']['ebv'], z=self.results[mtype]['CI']['z'], photbands=self.master['photband'][keep]) - + flux,flux_ur = flux*scale,flux_ur*scale - + synflux[keep] = synflux_ - + #synflux,Labs = model.get_itable(teff=self.results[mtype]['CI']['teff'], # logg=self.results[mtype]['CI']['logg'], # ebv=self.results[mtype]['CI']['ebv'], @@ -2035,16 +2038,16 @@ def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004',**kwargs): self.results[mtype]['model'] = wave,flux,flux_ur self.results[mtype]['synflux'] = eff_waves,synflux,self.master['photband'] self.results[mtype]['chi2'] = chi2 - + def set_model(self,wave,flux,label='manual'): """ Manually set the best SED model. - + This is particularly useful when working with calibrators for which there is an observed SED. - + The label will be used to store the model in the C{results} attribute. - + @param wave: wavelength array (angstrom) @type wave: ndarray @param flux: flux array (erg/s/cm2/AA) @@ -2057,24 +2060,24 @@ def set_model(self,wave,flux,label='manual'): is_color = self.master['color'] include = self.master['include'] synflux = np.zeros(len(photbands)) - + #-- compute synthetic photometry synflux[(-is_color) & include] = model.synthetic_flux(wave,flux,photbands[(-is_color) & include]) synflux[is_color & include] = model.synthetic_color(wave,flux,photbands[is_color & include]) chi2 = (self.master['cmeas']-synflux)**2/self.master['e_cmeas']**2 eff_waves = filters.eff_wave(photbands) - + if not label in self.results: self.results[label] = {} self.results[label]['model'] = wave,flux,flux self.results[label]['synflux'] = eff_waves,synflux,photbands self.results[label]['chi2'] = chi2 logger.debug('Stored model SED in {}'.format(label)) - + def get_model(self,label='igrid_search'): """ Retrieve the best SED model. - + B{Warning}: the grid search interpolation is also done with interpolation in metallicity, while this is not the case for the best full SED model. Interpolation of the full SED is only done in teff/logg, and the @@ -2090,27 +2093,27 @@ def get_model(self,label='igrid_search'): """ wave,flux,urflux = self.results[label]['model'] return wave,flux,urflux - + def sample_gridsearch(self,NrSamples=1,df=None,selfact='chisq'): """ Retrieve an element from the results of a grid search according to the derived probability. - + An array of length "NrSamples" containing 'grid-indices' is returned, so the actual parameter values of the corresponding model can be retrieved from the results dictionary. - + @param NrSamples: the number of samples you wish to draw @type NrSamples: int """ ranges = self.generate_ranges(teffrange=None,loggrange=None,ebvrange=None,\ zrange=None,rvrange=None,vradrange=None,start_from='igrid_search') - + #-- If nessessary calculate degrees of freedom from the ranges if df == None and ranges != None: df,df_info = self.calculateDF(**ranges) elif df == None: logger.warning('Cannot compute degrees of freedom!!! CHI2 might not make sense. (using df=5)') df = 5 - + #-- Do the statistics N = sum(self.master['include']) k = N-df @@ -2118,22 +2121,22 @@ def sample_gridsearch(self,NrSamples=1,df=None,selfact='chisq'): logger.warning('Not enough data to compute CHI2: it will not make sense') k = 1 logger.info('Statistics based on df={0} and Nobs={1}'.format(df,N)) - + #-- Compute the pdf and cdf probdensfunc = scipy.stats.distributions.chi2.pdf(self.results['igrid_search']['grid'][selfact],k) cumuldensfunc = probdensfunc.cumsum() - + #-- Uniformly sample the cdf, to get a sampling according to the pdf sample = pl.uniform(cumuldensfunc[0],cumuldensfunc[-1],NrSamples) indices = np.zeros(NrSamples,int) for i in xrange(NrSamples): indices[i] = (abs(cumuldensfunc-sample[i])).argmin() return indices - + def chi2(self,select=None,reduced=False,label='igrid_search'): """ Calculate chi2 of best model. - + TEMPORARY!!! API WILL CHANGE!!! """ #-- select photometry @@ -2147,30 +2150,30 @@ def chi2(self,select=None,reduced=False,label='igrid_search'): photbands = master['photband'] is_color = master['color'] synflux = np.zeros(len(photbands)) - + wave,flux,urflux = self.get_model(label=label) synflux[-is_color] = model.synthetic_flux(wave,flux,photbands[-is_color]) synflux[is_color] = model.synthetic_color(wave,flux,photbands[is_color]) x2 = (master['cmeas']-synflux)**2/master['e_cmeas']**2 x2 = x2.mean() return x2 - + #} #{ Add constraints - + def add_constraint_distance(self,distance=None,mtype='igrid_search',**kwargs): """ Use the distance to compute additional information. - + Compute radii, absolute luminosities and masses, and add them to the results. - + Extra kwargs go to L{get_distance_from_plx}. - + B{Warning:} after calling this function, the C{labs} column in the grid is actually absolutely calibrated and reflects the true absolute luminosity instead of the absolute luminosity assuming 1 solar radius. - + @param distance: distance in solar units and error @type distance: tuple (float,float) @param mtype: type of results to add the information to @@ -2180,7 +2183,7 @@ def add_constraint_distance(self,distance=None,mtype='igrid_search',**kwargs): kwargs['lutz_kelker'] = False # we can't handle asymmetric error bars kwargs['unit'] = 'Rsol' distance = self.get_distance_from_plx(**kwargs) - + #-- compute the radius, absolute luminosity and mass: note that there is # an uncertainty on the scaling factor *and* on the distance! scale = conversions.unumpy.uarray([self.results[mtype]['grid']['scale'],\ @@ -2190,7 +2193,7 @@ def add_constraint_distance(self,distance=None,mtype='igrid_search',**kwargs): labs = self.results[mtype]['grid']['labs']*radius**2 mass = conversions.derive_mass((self.results[mtype]['grid']['logg'],'[cm/s2]'),\ (radius,'Rsol'),unit='Msol') - + #-- store the results in the grid mass,e_mass = conversions.unumpy.nominal_values(mass),\ conversions.unumpy.std_devs(mass) @@ -2208,18 +2211,18 @@ def add_constraint_distance(self,distance=None,mtype='igrid_search',**kwargs): else: self.results[mtype]['grid'] = pl.mlab.rec_append_fields(self.results[mtype]['grid'],\ labels,data) - + #-- update the confidence intervals self.calculate_confidence_intervals(mtype=mtype) logger.info('Added constraint: distance (improved luminosity, added info on radius and mass)') - + def add_constraint_slo(self,numax,Deltanu0,mtype='igrid_search',chi2_type='red'): """ Use diagnostics from solar-like oscillations to put additional constraints on the parameters. - + If the results are constrained with the distance before L{add_constraint_distance}, then these results are combined with the SLO constraints. - + @param numax: frequency of maximum amplitude @type numax: 3-tuple (value,error,unit) @param Deltanu0: large separation (l=0) @@ -2234,17 +2237,17 @@ def add_constraint_slo(self,numax,Deltanu0,mtype='igrid_search',chi2_type='red') logg = grid['logg'] pvalues = [1-grid['ci_'+chi2_type]] names = ['ci_'+chi2_type+'_phot'] - + #-- we *always* have a logg, that's the way SEDs work: logg_slo = conversions.derive_logg_slo((teff,'K'),numax) logg_slo,e_logg_slo = conversions.unumpy.nominal_values(logg_slo),\ conversions.unumpy.std_devs(logg_slo) # compute the probabilities: scale to standard normal logg_prob = scipy.stats.distributions.norm.sf(abs(logg_slo-logg)/e_logg_slo) - + pvalues.append(logg_prob) - names.append('ci_logg_slo') - + names.append('ci_logg_slo') + #-- but we only sometimes have a radius (viz. when the distance is # known). There is an uncertainty on that radius! radius_slo = conversions.derive_radius_slo(numax,Deltanu0,(teff,'K'),unit='Rsol') @@ -2284,41 +2287,41 @@ def add_constraint_slo(self,numax,Deltanu0,mtype='igrid_search',chi2_type='red') grid = pl.mlab.rec_append_fields(grid,['radius','e_radius','e_labs','mass','e_mass','distance','e_distance'],\ [radius_slo,e_radius_slo,e_labs,mass,e_mass,distance,e_distance]) logger.info('Added constraint: {0:s} via slo, improved luminosity'.format(', '.join(['radius','e_radius','e_labs','mass','e_mass']))) - + #-- combine p values using Fisher's method combined_pvals = evaluate.fishers_method(pvalues) - + #-- add the information to the grid self.results[mtype]['grid'] = pl.mlab.rec_append_fields(grid,\ names,pvalues) - + #-- and replace the original confidence intervals, and re-order self.results[mtype]['grid']['ci_'+chi2_type] = 1-combined_pvals sa = np.argsort(self.results[mtype]['grid']['ci_'+chi2_type])[::-1] self.results[mtype]['grid'] = self.results[mtype]['grid'][sa] - + #-- update the confidence intervals self.calculate_confidence_intervals(mtype=mtype) logger.info('Added constraint: {0:s} via slo and replaced ci_{1:s} with combined CI'.format(', '.join(names),chi2_type)) - + def add_constraint_reddening(self,distance=None,ebv=None,e_ebv=0.1,Rv=3.1,\ model=None,mtype='igrid_search',chi2_type='red',upper_limit=False): """ Use reddening maps to put additional constraints on the parameters. - + This constraint assumes that, if C{upper_limit=False}, given the distance to target, the reddening from the reddening maps is similar to the one derived from the SED. All models where this is not the case will be deemed improbable. - + If C{upper_limit=True}, all models with an E(B-V) value above C{ebv} will be considered improbable. - + When you don't set C{ebv}, the value will, by default, be derived from Drimmel maps if C{upper_limit=False} and from Schlegel if C{upper_limit=True}. You can change this behaviour by setting C{model} manually. - + @param distance: distance and uncertainty in parsec @type distance: tuple (float,float) @param ebv: E(B-V) reddening in magnitude @@ -2360,22 +2363,22 @@ def add_constraint_reddening(self,distance=None,ebv=None,e_ebv=0.1,Rv=3.1,\ grid['ci_'+chi2_type] = 1-combined_pvals else: grid['ci_'+chi2_type] = np.where(ebvs<=ebv,grid['ci_'+chi2_type],1.) - + #-- and replace the original confidence intervals, and re-order sa = np.argsort(grid['ci_'+chi2_type])[::-1] self.results[mtype]['grid'] = grid[sa] - + #-- update the confidence intervals self.calculate_confidence_intervals(mtype=mtype) logger.info('Added constraint: E(B-V)={0}+/-{1}'.format(ebv,e_ebv)) - + def add_constraint_angular_diameter(self,angdiam): raise NotImplementedError - + def add_constraint_mass(self,mass,mtype='igrid_search',chi2_type='red'): """ Add constraints on the mass. - + C{mass} must be a tuple, if the second element is smaller than the first, it is assumed to be (mu,sigma) from a normal distribution. If the second element is larger than the first, it is assumed to be (lower, upper) @@ -2391,15 +2394,15 @@ def add_constraint_mass(self,mass,mtype='igrid_search',chi2_type='red'): grid['ci_'+chi2_type] = 1-combined_pvals else: grid['ci_'+chi2_type] = np.where((masses<=mass[0]) | (mass[1]<=masses),1.,grid['ci_'+chi2_type]) - + #-- and replace the original confidence intervals, and re-order sa = np.argsort(grid['ci_'+chi2_type])[::-1] self.results[mtype]['grid'] = grid[sa] - + #-- update the confidence intervals self.calculate_confidence_intervals(mtype=mtype) logger.info('Added constraint: mass={0}+/-{1}'.format(*mass)) - + def add_constraint_evolution_models(self,models='siess2000',\ ylabels=['age','labs','radius'],e_y=None, function='linear',mtype='igrid_search',chi2_type='red'): @@ -2407,10 +2410,10 @@ def add_constraint_evolution_models(self,models='siess2000',\ Use stellar evolutionary models to put additional constraints on the parameters. """ grid = self.results[mtype]['grid'] - + #-- make sure y's and ylabels are iterable if isinstance(ylabels,str): - ylabels = [ylabels] + ylabels = [ylabels] #-- cycle over all yvalues and compute the pvalues pvalues = [1-grid['ci_'+chi2_type]] add_info = [] @@ -2452,7 +2455,7 @@ def add_constraint_evolution_models(self,models='siess2000',\ pl.scatter(grid['teff'],grid['logg'],c=y_computed,edgecolors='none',cmap=pl.cm.spectral) pl.colorbar() pvalues.append(y_prob) - + pl.figure() pl.subplot(221) pl.title('p1') @@ -2475,7 +2478,7 @@ def add_constraint_evolution_models(self,models='siess2000',\ pl.colorbar() pl.xlabel('labs') pl.ylabel('radius') - + #-- combine p values using Fisher's method combined_pvals = evaluate.fishers_method(pvalues) pl.subplot(224) @@ -2485,8 +2488,8 @@ def add_constraint_evolution_models(self,models='siess2000',\ pl.colorbar() pl.xlabel('labs') pl.ylabel('radius') - pl.show() - + pl.show() + #-- add the information to the grid (assume that if there one label # not in there already, none of them are if not 'c_'+ylabels[0] in grid.dtype.names: @@ -2496,31 +2499,31 @@ def add_constraint_evolution_models(self,models='siess2000',\ else: for info,ylabel in zip(add_info,ylabels): self.results[mtype]['grid']['c_'+ylabel] = add_info - + #-- and replace the original confidence intervals, and re-order self.results[mtype]['grid']['ci_'+chi2_type] = 1-combined_pvals - + sa = np.argsort(self.results[mtype]['grid']['ci_'+chi2_type])[::-1] self.results[mtype]['grid'] = self.results[mtype]['grid'][sa] - + #-- update the confidence intervals self.calculate_confidence_intervals(mtype=mtype) - logger.info('Added constraint: {0:s} via stellar models and replaced ci_{1:s} with combined CI'.format(', '.join(ylabels),chi2_type)) - - + logger.info('Added constraint: {0:s} via stellar models and replaced ci_{1:s} with combined CI'.format(', '.join(ylabels),chi2_type)) + + #} - + #{ Plotting routines - + def _label_dict(self, param): """ Returns the label belonging to a certain parameter - + If the label is not present, the function will just return param. """ #split parameter in param name and componentent number param, component = re.findall('(.*?)(\d?$)', param)[0] - + ldict = dict(teff='Effective temperature [K]',\ z='log (Metallicity Z [$Z_\odot$]) [dex]',\ logg=r'log (surface gravity [cm s$^{-2}$]) [dex]',\ @@ -2539,31 +2542,31 @@ def _label_dict(self, param): return param + " - " + component else: return param - + @standalone_figure def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0.95,d=None,**kwargs): """ Plot grid as scatter plot - + PrameterC{ptype} sets the colors of the scattered points (e.g., 'ci_red','z','ebv'). - + Example usage: - + First set the SED: - + >>> mysed = SED('HD180642') >>> mysed.load_fits() - + Then make the plots: - + >>> p = pl.figure() >>> p = pl.subplot(221);mysed.plot_grid(limit=None) >>> p = pl.subplot(222);mysed.plot_grid(x='ebv',y='z',limit=None) >>> p = pl.subplot(223);mysed.plot_grid(x='teff',y='ebv',limit=None) >>> p = pl.subplot(224);mysed.plot_grid(x='logg',y='z',limit=None) - + ]]include figure]]ivs_sed_builder_plot_grid_01.png] - + """ #-- if no distance is set, derive the most likely distance from the plx: #if d is None: @@ -2584,7 +2587,7 @@ def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0 #(radius,'Rsol'),unit='Msol') #-- compute angular diameter theta = 2*conversions.convert('sr','mas',self.results[mtype]['grid']['scale']) - + if limit is not None: region = self.results[mtype]['grid']['ci_red']<=limit else: @@ -2599,7 +2602,7 @@ def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0 else: colors = ptype[region] ptype = 'custom variable' - + if 'ci_' in ptype.lower(): colors *= 100. vmin = colors.min() @@ -2607,23 +2610,23 @@ def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0 else: vmin = kwargs.pop('vmin',colors.min()) vmax = kwargs.pop('vmax',colors.max()) - + #-- define abbrevation for plotting if x in self.results[mtype]['grid'].dtype.names: X = self.results[mtype]['grid'][x] else: X = locals()[x] - + if y in self.results[mtype]['grid'].dtype.names: Y = self.results[mtype]['grid'][y] else: Y = locals()[y] - + #-- make the plot if mtype == 'imc': pl.hexbin(X,Y,mincnt=1,cmap=pl.cm.spectral) #bins='log' ptype = 'mc' - + #-- set the limits pl.xlim(X.max(),X.min()) pl.ylim(Y.max(),Y.min()) @@ -2640,7 +2643,7 @@ def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0 #colors = self.results[mtype]['grid'][ptype][region] #else: #colors = locals()[ptype][region] - + #if 'ci' in ptype: #colors *= 100. #vmin = colors.min() @@ -2648,7 +2651,7 @@ def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0 #else: #vmin = kwargs.pop('vmin',colors.min()) #vmax = kwargs.pop('vmax',colors.max()) - + #-- grid scatter plot pl.scatter(X[region],Y[region], c=colors,edgecolors='none',cmap=pl.cm.spectral,vmin=vmin,vmax=vmax) @@ -2656,28 +2659,28 @@ def plot_grid(self,x='teff',y='logg',ptype='ci_red',mtype='igrid_search',limit=0 pl.xlim(X[region].max(),X[region].min()) pl.ylim(Y[region].max(),Y[region].min()) cbar = pl.colorbar() - + #-- mark best value pl.plot(X[-1],Y[-1],'r+',ms=40,mew=3) - + pl.xlabel(self._label_dict(x)) pl.ylabel(self._label_dict(y)) cbar.set_label(self._label_dict(ptype)) - + logger.info('Plotted %s-%s diagram of %s'%(x,y,ptype)) - + @standalone_figure def plot_CI2D(self,xpar='teff',ypar='logg',mtype='iminimize', ptype='ci_red',**kwargs): """ - Plot a 2D confidence intervall calculated using the CI2D computation + Plot a 2D confidence intervall calculated using the CI2D computation from calculate_iminimize_CI2D. Make sure you first calculate the grid you want using the calculate_iminimize_CI2D function. """ - + grid = self.results[mtype]['CI2D'][xpar+"-"+ypar][ptype] x = self.results[mtype]['CI2D'][xpar+"-"+ypar]['x'] y = self.results[mtype]['CI2D'][xpar+"-"+ypar]['y'] - + if ptype == 'ci_red' or ptype == 'ci_raw': grid = grid*100.0 levels = np.linspace(0,100,25) @@ -2686,32 +2689,32 @@ def plot_CI2D(self,xpar='teff',ypar='logg',mtype='iminimize', ptype='ci_red',**k grid = np.log10(grid) levels = np.linspace(np.min(grid), np.max(grid), 25) ticks = np.round(np.linspace(np.min(grid), np.max(grid), 11), 2) - + pl.contourf(x,y,grid,levels,**kwargs) pl.xlabel(self._label_dict(xpar)) pl.ylabel(self._label_dict(ypar)) cbar = pl.colorbar(fraction=0.08,ticks=ticks) cbar.set_label(ptype!='ci_chi2' and r'Probability' or r'$^{10}$log($\chi^2$)') - - - - + + + + @standalone_figure def plot_data(self,colors=False, plot_unselected=True, unit_wavelength='angstrom',unit_flux=None,**kwargs): """ Plot only the SED data. - + Extra kwargs are passed to plotting functions. - + The decorator provides an keyword C{savefig}. When set to C{True}, an image name is generated automatically (very long!), and the figure is closed. When set to a string, the image is saved with that name, and the figure is closed. When C{savefig} is not given, the image will stay open so that the user can still access all plot elements and do further enhancements. - + @param colors: if False, plot absolute values, otherwise plot colors (flux ratios) @type colors: boolean @@ -2729,11 +2732,11 @@ def plot_data(self,colors=False, plot_unselected=True, iscolor = np.array(master['color'],bool) photbands = master['photband'] indices = np.arange(len(master)) - + allsystems = np.array([i.split('.')[0] for i in photbands]) systems = sorted(set(allsystems)) color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(systems))] - + if not colors: color_cycle = itertools.cycle(color_cycle) pl.gca().set_xscale('log',nonposx='clip') @@ -2783,32 +2786,32 @@ def plot_data(self,colors=False, plot_unselected=True, leg.get_frame().set_alpha(0.5) pl.grid() pl.title(self.ID) - - + + @standalone_figure def plot_sed(self,colors=False,mtype='igrid_search',plot_deredded=False, plot_unselected=True,wave_units='AA',flux_units='erg/s/cm2',**kwargs): """ Plot a fitted SED together with the data. - + Example usage: - + First set the SED: - + >>> mysed = SED('HD180642') >>> mysed.load_fits() - + Then make the plots: - + >>> p = pl.figure() >>> p = pl.subplot(121) >>> mysed.plot_sed(colors=False) >>> p = pl.subplot(122) >>> mysed.plot_sed(colors=True) - + ]]include figure]]ivs_sed_builder_plot_sed_01.png] """ - + def plot_sed_getcolors(master,color_dict=None): myphotbands = [iphotb.split('.')[1] for iphotb in master['photband'][master['color']]] if not myphotbands: #-- If there are no colours none can be returned (added by joris 30-01-2012) @@ -2823,15 +2826,15 @@ def plot_sed_getcolors(master,color_dict=None): y = master['cmeas'] e_y = master['e_cmeas'] return x,y,e_y,color_dict - - + + x__,y__,e_y__,color_dict = plot_sed_getcolors(self.master) - + #-- get the color cycle systems = np.array([system.split('.')[0] for system in self.master['photband']],str) set_systems = sorted(list(set(systems))) color_cycle = itertools.cycle([pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(set_systems))]) - + #-- for plotting reasons, we translate every color to an integer for system in set_systems: color = color_cycle.next() @@ -2864,12 +2867,12 @@ def plot_sed_getcolors(master,color_dict=None): y = self.master['cmeas'][keep] e_y = self.master['e_cmeas'][keep] y,e_y = conversions.convert('erg/s/cm2/AA',flux_units,y,e_y,wave=(x,'AA')) - x = conversions.convert('AA',wave_units,x) - - + x = conversions.convert('AA',wave_units,x) + + p = pl.errorbar(x,y,yerr=e_y,fmt='o',label=system, capsize=10,ms=7,color=color,**kwargs) - + #-- exclude: label = np.any(keep) and '_nolegend_' or system keep = (systems==system) & (self.master['color']==colors) & -self.master['include'] @@ -2883,25 +2886,25 @@ def plot_sed_getcolors(master,color_dict=None): y = self.master['cmeas'][keep] e_y = self.master['e_cmeas'][keep] y,e_y = conversions.convert('erg/s/cm2/AA',flux_units,y,e_y,wave=(x,'AA')) - x = conversions.convert('AA',wave_units,x) + x = conversions.convert('AA',wave_units,x) pl.errorbar(x,y,yerr=e_y,fmt='o',label=label, capsize=10,ms=7,mew=2,color=color,mfc='1',mec=color,**kwargs) - + #-- only set logarithmic scale if absolute fluxes are plotted # and only plot the real model then if not colors: pl.gca().set_xscale('log',nonposx='clip') pl.gca().set_yscale('log',nonposy='clip') pl.gca().set_autoscale_on(False) - + #-- the model if mtype in self.results and 'model' in self.results[mtype]: wave,flux,flux_ur = self.results[mtype]['model'] - + flux = conversions.convert('erg/s/cm2/AA',flux_units,flux,wave=(wave,'AA')) flux_ur = conversions.convert('erg/s/cm2/AA',flux_units,flux_ur,wave=(wave,'AA')) - wave = conversions.convert('AA',wave_units,wave) - + wave = conversions.convert('AA',wave_units,wave) + pl.plot(wave,flux,'r-',**kwargs) if plot_deredded: pl.plot(wave,flux_ur,'k-',**kwargs) @@ -2914,7 +2917,7 @@ def plot_sed_getcolors(master,color_dict=None): pl.ylabel(r'Flux ratio') pl.xlabel('Color name') pl.xlim(min(xticks)-0.5,max(xticks)+0.5) - + pl.grid() if colors: leg = pl.legend(loc='best',prop=dict(size='x-small')) @@ -2938,38 +2941,38 @@ def plot_sed_getcolors(master,color_dict=None): pl.annotate('Teff=%d K\nlogg=%.2f cgs\nE(B-V)=%.3f mag\n$\Theta$=%.3g mas'%(teff,logg,ebv,angdiam), loc,xycoords='axes fraction') logger.info('Plotted SED as %s'%(colors and 'colors' or 'absolute fluxes')) - - @standalone_figure + + @standalone_figure def plot_chi2(self,colors=False,mtype='igrid_search'): """ Plot chi2 statistic for every datapoint included in the fit. - + To plot the statistic from the included absolute values, set C{colors=False}. - + To plot the statistic from the included colors, set C{colors=True}. - + Example usage: - + First set the SED: - + >>> mysed = SED('HD180642') >>> mysed.load_fits() - + Then make the plots: - + >>> p = pl.figure() >>> p = pl.subplot(121) >>> mysed.plot_chi2(colors=False) >>> p = pl.subplot(122) >>> mysed.plot_chi2(colors=True) - + ]]include figure]]ivs_sed_builder_plot_chi2_01.png] - + @param colors: flag to distinguish between colors and absolute values @type colors: boolean """ - + include_grid = self.master['include'] systems = np.array([system.split('.')[0] for system in self.master['photband'][include_grid]],str) set_systems = sorted(list(set(systems))) @@ -3001,17 +3004,17 @@ def plot_chi2(self,colors=False,mtype='igrid_search'): logger.info('Plotted CHI2 of %s'%(colors and 'colors' or 'absolute fluxes')) else: logger.info('%s not in results, no plot could be made.'%(mtype)) - - - @standalone_figure + + + @standalone_figure def plot_distance(self,mtype='igrid_search'): try: #-- necessary information (d_models,d_prob_models,radii) = self.results['igrid_search']['d_mod'] (d,dprob) = self.results['igrid_search']['d'] - + ax_d = pl.gca() - + gal = self.info['galpos'] #-- the plot dzoom = dprob>1e-4 @@ -3041,7 +3044,7 @@ def plot_distance(self,mtype='igrid_search'): logger.info('Plotted distance/reddening') except KeyError: logger.info('No distance/reddening plotted due to KeyError.') - + @standalone_figure def plot_grid_model(self,ptype='prob'): """ @@ -3053,7 +3056,7 @@ def plot_grid_model(self,ptype='prob'): (d_models,d_prob_models,radii) = self.results['igrid_search']['d_mod'] (d,dprob) = self.results['igrid_search']['d'] gal = self.info['galpos'] - + n = 75000 region = self.results['igrid_search']['grid']['ci_red']<0.95 total_prob = 100-(1-self.results['igrid_search']['grid']['ci_red'][cutlogg][-n:])*d_prob_models*100 @@ -3071,17 +3074,17 @@ def plot_grid_model(self,ptype='prob'): cbar = pl.colorbar() pl.xlabel('log (effective temperature [K]) [dex]') pl.ylabel(r'log (surface gravity [cm s$^{-2}$]) [dex]') - + if ptype=='prob': cbar.set_label('Probability (incl. plx) [%]') elif ptype=='radii': cbar.set_label('Model radii [$R_\odot$]') logger.info('Plotted teff-logg diagram of models (%s)'%(ptype)) - - + + @standalone_figure def plot_MW_side(self): - im = Image.open(config.get_datafile('images','NRmilkyway.tif')) + im = image.open(config.get_datafile('images','NRmilkyway.tif')) left,bottom,width,height = 0.0,0.0,1.0,1.0 startm,endm = 183,-177 startv,endv = -89,91 @@ -3100,10 +3103,10 @@ def plot_MW_side(self): pl.ylim(startv,endv) pl.xticks([]) pl.yticks([]) - + @standalone_figure def plot_MW_top(self): - im = Image.open(config.get_datafile('images','topmilkyway.jpg')) + im = image.open(config.get_datafile('images','topmilkyway.jpg')) pl.imshow(im,origin='lower') pl.box(on=False) pl.xticks([]) @@ -3113,7 +3116,7 @@ def plot_MW_top(self): gal = self.info['galpos'] pl.plot(2800,1720,'ro',ms=10) pl.plot([2800,-5000*np.sin(gal[0]/180.*np.pi)+2800],[1720,5000*np.cos(gal[0]/180.*np.pi)+1720],'r-',lw=2) - + #-- necessary information if 'igrid_search' in self.results and 'd_mod' in self.results['igrid_search']: (d_models,d_prob_models,radii) = self.results['igrid_search']['d_mod'] @@ -3121,20 +3124,20 @@ def plot_MW_top(self): else: d = np.linspace(0,1000,2) dprob = np.zeros(len(d)) - + x = d/10.*1.3 y = dprob*1000. theta = gal[0]/180.*np.pi + np.pi/2. x_ = np.cos(theta)*x - np.sin(theta)*y + 2800 y_ = np.sin(theta)*x + np.cos(theta)*y + 1720 - + pl.plot(x_,y_,'r-',lw=2) index = np.argmax(y) pl.plot(np.cos(theta)*x[index] + 2800,np.sin(theta)*x[index] + 1720,'rx',ms=15,mew=2) - + pl.xlim(xlims) pl.ylim(ylims) - + @standalone_figure def plot_finderchart(self,cmap_photometry=pl.cm.spectral,window_size=5.): """ @@ -3153,7 +3156,7 @@ def plot_finderchart(self,cmap_photometry=pl.cm.spectral,window_size=5.): toplot = self.master[keep_this] systems = np.array([system.split('.')[0] for system in toplot['photband']],str) set_systems = sorted(list(set(systems))) - color_cycle = itertools.cycle([cmap_photometry(j) for j in np.linspace(0, 1.0, len(set_systems))]) + color_cycle = itertools.cycle([cmap_photometry(j) for j in np.linspace(0, 1.0, len(set_systems))]) for system in set_systems: color = color_cycle.next() keep = systems==system @@ -3166,11 +3169,11 @@ def plot_finderchart(self,cmap_photometry=pl.cm.spectral,window_size=5.): pl.xlim(xlims) pl.ylim(ylims) pl.xlabel(r'Right ascension $\alpha$ [arcmin]') - pl.ylabel(r'Declination $\delta$ [arcmin]') + pl.ylabel(r'Declination $\delta$ [arcmin]') except: logger.warning('No image found of %s'%(self.ID)) pass - + if 'pm' in self.info: logger.info("Found proper motion info") ppm_ra,ppm_de = (self.info['pm']['pmRA'],self.info['pm']['epmRA']),(self.info['pm']['pmDE'],self.info['pm']['epmDE']) @@ -3191,9 +3194,9 @@ def plot_finderchart(self,cmap_photometry=pl.cm.spectral,window_size=5.): else: max_distance = 1000. e_max_distance = 100. - + tang_velo = 'Tan. vel. at %.0f+/-%.0f pc: '%(max_distance,e_max_distance) - + max_distance = conversions.convert('pc','km',max_distance,e_max_distance) ppm_ra = conversions.convert('mas/yr','rad/s',*ppm_ra) ppm_de = conversions.convert('mas/yr','rad/s',*ppm_de) @@ -3201,16 +3204,16 @@ def plot_finderchart(self,cmap_photometry=pl.cm.spectral,window_size=5.): x = unumpy.uarray([ppm_ra[0],ppm_ra[1]]) y = unumpy.uarray([ppm_de[0],ppm_de[1]]) velocity = max_distance*utan( usqrt(x**2+y**2)) - - + + pl.annotate(tang_velo + '%s km/s'%(velocity),xy=(0.05,0.2),xycoords='axes fraction',color='red') if 'Vel' in self.info and 'v' in self.info['Vel']: rad_velo = 'Rad. vel.: %.1f'%(self.info['Vel']['v']) if 'e' in self.info['Vel']: rad_velo += '+/-%.1f'%(self.info['Vel']['e']) pl.annotate(rad_velo+' km/s',xy=(0.05,0.15),xycoords='axes fraction',color='red') - - + + def make_plots(self): """ Make all available plots @@ -3223,29 +3226,29 @@ def make_plots(self): pl.subplot(rows,cols,2);self.plot_grid(ptype='ebv') pl.subplot(rows,cols,3);self.plot_grid(ptype='z') pl.subplot(rows,cols,4);self.plot_distance() - + pl.subplot(3,2,3);self.plot_sed(colors=False) pl.subplot(3,2,5);self.plot_sed(colors=True) - + pl.subplot(rows,cols,7);self.plot_chi2(colors=False) pl.subplot(rows,cols,11);self.plot_chi2(colors=True) - + #pl.subplot(rows,cols,8);self.plot_grid_model(ptype='prob') #pl.subplot(rows,cols,12);self.plot_grid_model(ptype='radii') - + pl.figure(figsize=(12,12)) pl.axes([0,0.0,1.0,0.5]);self.plot_MW_side() pl.axes([0,0.5,0.5,0.5]);self.plot_MW_top() pl.axes([0.5,0.5,0.5,0.5]);self.plot_finderchart() - + #} - + #{Input and output - + def save_photometry(self,photfile=None): """ Save master photometry to a file. - + @param photfile: name of the photfile. Defaults to C{starname.phot}. @type photfile: str """ @@ -3260,11 +3263,11 @@ def save_photometry(self,photfile=None): if not 'comments' in self.master.dtype.names: self.master = vizier.quality_check(self.master,self.ID) ascii.write_array(self.master,self.photfile,header=True,auto_width=True,use_float='%g',comments=['#'+json.dumps(self.info)]) - + def load_photometry(self,photfile=None): """ Load the contents of the photometry file to the master record. - + @param photfile: name of the photfile. Defaults to the value of C{self.photfile}. @type photfile: str """ @@ -3273,13 +3276,13 @@ def load_photometry(self,photfile=None): logger.info('Load photometry from file %s'%(self.photfile)) self.master,comments = ascii.read2recarray(self.photfile,return_comments=True) self.info = json.loads(comments[-3]) - - - + + + def save_fits(self,filename=None,overwrite=True): """ - Save content of SED object to a FITS file. - + Save content of SED object to a FITS file. + The .fits file will contain the following extensions if they are present in the object: 1) data (all previously collected photometric data = content of .phot file) 2) model_igrid_search (full best model SED table from grid_search) @@ -3288,14 +3291,14 @@ def save_fits(self,filename=None,overwrite=True): 5) model_imc (full best model SED table from monte carlo) 6) imc (CI from monte carlo = actual monte carlo samples) 7) synflux_imc (integrated synthetic fluxes for best model from monte carlo) - + @param filename: name of SED FITS file @type filename: string @param overwrite: overwrite old FITS file if true @type overwrite: boolean - + Example usage: - + >>> #mysed.save_fits() >>> #mysed.save_fits(filename='myname.fits') """ @@ -3305,7 +3308,7 @@ def save_fits(self,filename=None,overwrite=True): if os.path.isfile(filename): os.remove(filename) logger.info('Old FITS file removed') - + #-- write primary header #prim_header = {} #for key in self.info: @@ -3313,16 +3316,16 @@ def save_fits(self,filename=None,overwrite=True): #continue #prim_header[key] = self.info[key] #fits.write_recarray(np.array([[0]]),filename,header_dict=prim_header,ext=0) - + #-- write master data master = self.master.copy() fits.write_recarray(master,filename,header_dict=dict(extname='data')) - + #-- write the rest for mtype in self.results:#['igrid_search','imc']: eff_waves,synflux,photbands = self.results[mtype]['synflux'] chi2 = self.results[mtype]['chi2'] - + results_modeldict = dict(extname='model_'+mtype) results_griddict = dict(extname=mtype) keys = sorted(self.results[mtype]) @@ -3331,29 +3334,29 @@ def save_fits(self,filename=None,overwrite=True): for ikey in self.results[mtype][key]: if '_l' not in ikey and '_u' not in ikey and ikey != 'chisq': results_modeldict[ikey] = self.results[mtype][key][ikey] - results_griddict[ikey] = self.results[mtype][key][ikey] + results_griddict[ikey] = self.results[mtype][key][ikey] if key=='factor': results_griddict[key] = self.results[mtype][key] - + fits.write_array(list(self.results[mtype]['model']),filename, names=('wave','flux','dered_flux'), units=('AA','erg/s/cm2/AA','erg/s/cm2/AA'), header_dict=results_modeldict) if 'grid' in self.results[mtype]: fits.write_recarray(self.results[mtype]['grid'],filename,header_dict=results_griddict) - + results = np.rec.fromarrays([synflux,eff_waves,chi2],dtype=[('synflux','f8'),('mod_eff_wave','f8'),('chi2','f8')]) - + fits.write_recarray(results,filename,header_dict=dict(extname='synflux_'+mtype)) - + logger.info('Results saved to FITS file: %s'%(filename)) - - + + def load_fits(self,filename=None): """ - Load a previously made SED FITS file. Only works for SEDs saved with + Load a previously made SED FITS file. Only works for SEDs saved with the save_fits function after 14.06.2012. - + The .fits file can contain the following extensions: 1) data (all previously collected photometric data = content of .phot file) 2) model_igrid_search (full best model SED table from grid_search) @@ -3362,7 +3365,7 @@ def load_fits(self,filename=None): 5) model_imc (full best model SED table from monte carlo) 6) imc (CI from monte carlo = actual monte carlo samples) 7) synflux_imc (integrated synthetic fluxes for best model from monte carlo) - + @param filename: name of SED FITS file @type filename: string @rtype: bool @@ -3374,7 +3377,7 @@ def load_fits(self,filename=None): logger.warning('No previous results saved to FITS file {:s}'.format(filename)) return False ff = pyfits.open(filename) - + #-- observed photometry fields = ff['data'].columns.names master = np.rec.fromarrays([ff['data'].data.field(field) for field in fields],names=','.join(fields)) @@ -3384,11 +3387,11 @@ def load_fits(self,filename=None): for j,element in enumerate(master[name]): master[name][j] = element.strip() self.master = master - + #-- add dictionary that will contain the results if not hasattr(self,'results'): self.results = {} - + #-- grid search and MC results mtypes = [ext.header['extname'] for ext in ff[1:]] mtypes = list(set(mtypes) - set(['DATA'])) @@ -3399,8 +3402,8 @@ def load_fits(self,filename=None): self.results[mtype]['model'] = np.array(ff['model_'+mtype].data.field('wave'),dtype='float64'),np.array(ff['model_'+mtype].data.field('flux'),dtype='float64'),np.array(ff['model_'+mtype].data.field('dered_flux'),dtype='float64') self.results[mtype]['chi2'] = np.array(ff['synflux_'+mtype].data.field('chi2'),dtype='float64') self.results[mtype]['synflux'] = np.array(ff['synflux_'+mtype].data.field('mod_eff_wave'),dtype='float64'),np.array(ff['synflux_'+mtype].data.field('synflux'),dtype='float64'),self.master['photband'] - - + + for mtype in ['igrid_search','iminimize','imc']: try: fields = ff[mtype].columns.names @@ -3410,7 +3413,7 @@ def load_fits(self,filename=None): self.results[mtype]['grid'] = master if 'factor' in ff[mtype].header: self.results[mtype]['factor'] = np.array([ff[mtype].header['factor']])[0] - + headerkeys = ff[mtype].header.ascardlist().keys() for key in headerkeys[::-1]: for badkey in ['xtension','bitpix','naxis','pcount','gcount','tfields','ttype','tform','tunit','factor','extname']: @@ -3423,25 +3426,25 @@ def load_fits(self,filename=None): self.results[mtype]['CI'][key.lower()] = np.array([ff[mtype].header[key]])[0] except KeyError: continue - + #self.results['igrid_search'] = {} #fields = ff['igrid_search'].columns.names #master = np.rec.fromarrays([ff['igrid_search'].data.field(field) for field in fields],names=','.join(fields)) #self.results['igrid_search']['grid'] = master #self.results['igrid_search']['factor'] = ff['igrid_search'].header['factor'] - + #self.results['model'] = ff[2].data.field('wave'),ff[2].data.field('flux'),ff[2].data.field('dered_flux') #self.results['chi2'] = ff[4].data.field('chi2') #self.results['synflux'] = ff[4].data.field('mod_eff_wave'),ff[4].data.field('synflux'),ff[1].data.field('photband') - + ff.close() - + logger.info('Loaded previous results from FITS file: %s'%(filename)) return filename - + def save_hdf5(self, filename=None, update=True): """ - Save content of SED object to a HDF5 file. (HDF5 is the successor of FITS files, + Save content of SED object to a HDF5 file. (HDF5 is the successor of FITS files, providing a clearer structure of the saved content.) This way of saving is more thorough that save_fits(), fx. the CI2D confidence intervals are not save to a fits file, but are saved to a hdf5 file. @@ -3449,7 +3452,7 @@ def save_hdf5(self, filename=None, update=True): - sed.master (photometry) - sed.results (results from all fitting methods) - sed.constraints (extra constraints on the fits) - + @param filename: name of SED FITS file @type filename: string @param update: if True, an existing file will be updated with the current information, if @@ -3458,24 +3461,24 @@ def save_hdf5(self, filename=None, update=True): @return: the name of the output HDF5 file. @rtype: string """ - + if filename is None: filename = str(os.path.splitext(self.photfile)[0]+'.hdf5') - + data = dict() data['master'] = self.master data['results'] = self.results data['constraints'] = self.constraints - + hdf5.write_dict(data, filename, update=update) - + logger.info('Results saved to HDF5 file: %s'%(filename)) - return filename - + return filename + def load_hdf5(self,filename=None): """ Load a previously made SED from HDF5 file. - + @param filename: name of SED FITS file @type filename: string @return: True if HDF5 file could be loaded @@ -3486,27 +3489,27 @@ def load_hdf5(self,filename=None): if not os.path.isfile(filename): logger.warning('No previous results saved to HFD5 file {:s}'.format(filename)) return False - + data = hdf5.read2dict(filename) - + self.master = data.get('master', {}) self.results = data.get('results', {}) self.constraints = data.get('constraints', {}) - + logger.info('Loaded previous results from HDF5 file: %s'%(filename)) logger.debug('Loaded following datasets from HDF5 file:\n %s'%(data.keys())) return True - + def save_bibtex(self): """ Convert the bibcodes in a phot file to a bibtex file. - + The first line in the bibtex file contains a \citet command citing all photometry. """ filename = os.path.splitext(self.photfile)[0]+'.bib' crossmatch.make_bibtex(self.master,filename=filename) - + def save_summary(self,filename=None,CI_limit=None,method='igrid_search',chi2type='ci_red'): """ Save a summary of the results to an ASCII file. @@ -3514,10 +3517,10 @@ def save_summary(self,filename=None,CI_limit=None,method='igrid_search',chi2type #-- open the summary file to write the results if filename is None: filename = os.path.splitext(self.photfile)[0]+'.sum' - + if CI_limit is None: CI_limit = self.CI_limit - + #-- gather the results: grid_results = self.results[method]['grid'] start_CI = np.argmin(np.abs(grid_results[chi2type]-self.CI_limit)) @@ -3539,19 +3542,19 @@ def save_summary(self,filename=None,CI_limit=None,method='igrid_search',chi2type used_photbands = '#'+photbands+'\n' used_references = '#'+references comments = used_photometry+used_atmosphere+used_photbands+used_references - + contents = np.array([results]).T contents = np.rec.fromarrays(contents,names=names) ascii.write_array(contents,filename,auto_width=True,header=True, comments=comments.split('\n'),mode='a',use_float='%g') - + logger.info('Saved summary to {0}'.format(filename)) - + def ci2str(self, mtype='igrid_search'): """ Prints the stored confidence intervals of the given mtype, the confidence intervals are sorted alfabetically. - + @param mtype: the search type for which to show the ci @type mtype: str @return: The confidence intervals in a string @@ -3562,7 +3565,7 @@ def ci2str(self, mtype='igrid_search'): return res else: res = "Stored confidence intervals for type: %s \n"%(mtype) - + ci = self.results[mtype]['CI'] keys = [ key for key in ci.keys() if not '_' in key ] keys.sort() @@ -3571,13 +3574,13 @@ def ci2str(self, mtype='igrid_search'): res += "%15s : %g <= %g <= %g\n"%(key, ci[key+'_l'], ci[key], ci[key+'_u']) except Exception: res += "%15s : nan <= %g <= nan\n"%(key, ci[key]) - - return res.rstrip() - + + return res.rstrip() + #} class BinarySED(SED): - + def __init__(self,ID=None,photfile=None,plx=None,load_fits=True,load_hdf5=True,label='', **kwargs): """ Setup the Binary sed in the same way as a normal SED. @@ -3586,10 +3589,10 @@ def __init__(self,ID=None,photfile=None,plx=None,load_fits=True,load_hdf5=True,l """ super(BinarySED, self).__init__(ID=ID,photfile=photfile,plx=plx,label=label,\ load_fits=load_fits,load_hdf5=load_hdf5) - + self.set_constraints(**kwargs) - - + + def set_constraints(self, **kwargs): """ Add constraints that are used when fitting the Binary SED. Up till now the following @@ -3601,7 +3604,7 @@ def set_constraints(self, **kwargs): for key in kwargs.keys(): if kwargs[key] != None: self.constraints[key] = kwargs[key] - + def constraints2str(self): """ Summarizes all constraints in a string. @@ -3611,24 +3614,24 @@ def constraints2str(self): res += "Using constraint: %s = %s\n"%(key, self.constraints[key]) res = res[:-1] return res - + def igrid_search(self,points=100000,teffrange=None,loggrange=None,ebvrange=None,\ zrange=None,rvrange=((3.1,3.1),(3.1,3.1)),vradrange=((0,0),(0,0)),\ radrange=(None,None),compare=True,df=None,CI_limit=None,\ set_model=True, distance=None,**kwargs): """ Fit fundamental parameters using a (pre-integrated) grid search. - + If called consecutively, the ranges will be set to the CI_limit of previous estimations, unless set explicitly. - + If called for the first time, the ranges will be +/- np.inf by defaults, unless set explicitly. """ - + if CI_limit is None or CI_limit > 1.0: CI_limit = self.CI_limit - + #-- set defaults limits ranges = self.generate_ranges(teffrange=teffrange[0],\ loggrange=loggrange[0],ebvrange=ebvrange[0],\ @@ -3637,51 +3640,51 @@ def igrid_search(self,points=100000,teffrange=None,loggrange=None,ebvrange=None, logg2range=loggrange[1],ebv2range=ebvrange[1],\ z2range=zrange[1],rv2range=rvrange[1],vrad2range=vradrange[1], rad2range=radrange[1]) - + #-- grid search on all include data: extract the best CHI2 include_grid = self.master['include'] logger.info('The following measurements are included in the fitting process:\n%s'%\ (photometry2str(self.master[include_grid]))) - + logger.info('The following constraints are included in the fitting process:\n%s'%\ (self.constraints2str())) #-- build the grid, run over the grid and calculate the CHI2 masses = self.constraints.get('masses', None) - pars = fit.generate_grid_pix(self.master['photband'][include_grid], masses=masses, points=points, **ranges) - + pars = fit.generate_grid_pix(self.master['photband'][include_grid], masses=masses, points=points, **ranges) + chisqs,scales,escales,lumis = fit.igrid_search_pix(self.master['cmeas'][include_grid], self.master['e_cmeas'][include_grid], - self.master['photband'][include_grid], + self.master['photband'][include_grid], model_func=model.get_itable_pix,constraints=self.constraints, **pars) fitres = dict(chisq=chisqs, scale=scales, escale=escales, labs=lumis) - + #-- collect all the results in a record array self.collect_results(grid=pars, fitresults=fitres, mtype='igrid_search') - + #-- do the statistics self.calculate_statistics(df=df, ranges=ranges, mtype='igrid_search') - + #-- compute the confidence intervals ci = self.calculate_confidence_intervals(mtype='igrid_search',chi2_type='red',\ CI_limit=CI_limit) self.store_confidence_intervals(mtype='igrid_search', **ci) - + #-- remember the best model if set_model: self.set_best_model() - + def generate_fit_param(self, start_from='igrid_search', **pars): - """ - generates a dictionary with parameter information that can be handled by fit.iminimize + """ + generates a dictionary with parameter information that can be handled by fit.iminimize """ masses = self.constraints.get('masses', None) result = super(BinarySED, self).generate_fit_param(start_from=start_from, **pars) - + #-- Couple ebv and rv of both components result['ebv2_expr'] = 'ebv' result['rv2_expr'] = 'rv' - + if masses != None: #-- Couple the radii to the masses G, Msol, Rsol = constants.GG_cgs, constants.Msol_cgs, constants.Rsol_cgs @@ -3695,17 +3698,17 @@ def generate_fit_param(self, start_from='igrid_search', **pars): else: result['rad_value'] = self.results[start_from]['CI']['rad'] result['rad2_value'] = self.results[start_from]['CI']['rad2'] - + return result - + def iminimize(self, teff=(None,None), logg=(None,None), ebv=(None,None), z=(None,None), rv=(None,None), vrad=(0,0), teffrange=(None,None), loggrange=(None,None), ebvrange=(None,None), zrange=(None,None), rvrange=(None,None), - vradrange=(None,None), radrange=(None,None), compare=True, df=None, + vradrange=(None,None), radrange=(None,None), compare=True, df=None, distance=None, start_from='igrid_search', points=None, CI_limit=None, calc_ci=True, set_model=True, **kwargs): - """ Binary minimizer """ - + """ Binary minimizer """ + ranges = self.generate_ranges(teffrange=teffrange[0],\ loggrange=loggrange[0],ebvrange=ebvrange[0],\ zrange=zrange[0],rvrange=rvrange[0],vradrange=vradrange[0], @@ -3717,47 +3720,47 @@ def iminimize(self, teff=(None,None), logg=(None,None), ebv=(None,None), z=(None rv=rv[0], vrad=vrad[0],teff2=teff[1],logg2=logg[1],\ ebv2=ebv[1],z2=z[1],rv2=rv[1],vrad2=vrad[1], \ start_from=start_from, **ranges) - + #-- Print the used data and constraints include_grid = self.master['include'] logger.info('The following measurements are included in the fitting process:\n%s'%\ (photometry2str(self.master[include_grid]))) logger.info('The following constraints are included in the fitting process:\n%s'%\ (self.constraints2str())) - + #-- pass all ranges and starting values to the fitter kick_list = ['teff', 'logg', 'teff2', 'logg2', 'ebv'] grid, chisq, nfev, scale, lumis = fit.iminimize(self.master['cmeas'][include_grid], self.master['e_cmeas'][include_grid], self.master['photband'][include_grid], points=points, kick_list=kick_list, constraints=self.constraints, **pars) - + logger.info('Minimizer Succes with startpoints=%s, chi2=%s, nfev=%s'%(len(chisq), chisq[0], nfev[0])) #-- handle the results fitres = dict(chisq=chisq, nfev=nfev, scale=scale, labs=lumis) self.collect_results(grid=grid, fitresults=fitres, mtype='iminimize', selfact='chisq') - + #-- Do the statistics self.calculate_statistics(df=5, ranges=ranges, mtype='iminimize', selfact='chisq') - + #-- Store the confidence intervals ci = self._get_imin_ci(mtype='iminimize',**ranges) self.store_confidence_intervals(mtype='iminimize', **ci) if calc_ci: self.calculate_iminimize_CI(mtype='iminimize', CI_limit=CI_limit) - + #-- remember the best model - if set_model: self.set_best_model(mtype='iminimize') - - + if set_model: self.set_best_model(mtype='iminimize') + + def calculate_iminimize_CI(self, mtype='iminimize', CI_limit=0.66, **kwargs): """ Calculate the confidence intervals for each parameter using the lmfit - calculate confidence interval method. - - The calculated confidence intervals are stored in the results['CI'] - dictionary. If the method fails, or if the asked CI is outside the + calculate confidence interval method. + + The calculated confidence intervals are stored in the results['CI'] + dictionary. If the method fails, or if the asked CI is outside the provided ranges, those ranges will be set as CI. - + This method works in the same way as for a single SED, but it adds the mass as an extra constraint if the mass of both components is stored. """ @@ -3765,37 +3768,37 @@ def calculate_iminimize_CI(self, mtype='iminimize', CI_limit=0.66, **kwargs): super(BinarySED, self).calculate_iminimize_CI(mtype=mtype, CI_limit=CI_limit,\ masses=masses, constraints=self.constraints,\ **kwargs) - + def calculate_iminimize_CI2D(self,xpar, ypar, mtype='iminimize', limits=None, res=10, **kwargs): """ Calculated 2 dimentional confidence intervals for the given parameters, - using lmfit methods. - - The calculated confidence intervals are stored in the results['CI2D'] + using lmfit methods. + + The calculated confidence intervals are stored in the results['CI2D'] dictionary. - + This method works in the same way as for a single SED, but it adds the mass as an extra constraint if the mass of both components is stored. """ masses = self.constraints.get('masses',None) super(BinarySED, self).calculate_iminimize_CI2D(xpar, ypar, mtype=mtype,\ limits=limits, res=res, **kwargs) #,masses=masses) - + def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004', **kwargs): """ Get reddenend and unreddened model """ logger.info('Interpolating approximate full SED of best model') - + #-- synthetic flux include = self.master['include'] synflux = np.zeros(len(self.master['photband'])) keep = (self.master['cwave']<1.6e6) | np.isnan(self.master['cwave']) keep = keep & include - + if mtype in ['igrid_search', 'iminimize']: scale = self.results[mtype]['CI']['scale'] - + #-- get parameters pars, pars_ur = {}, {} for key in self.results[mtype]['CI'].keys(): @@ -3805,7 +3808,7 @@ def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004', **kwargs): pars_ur[key] = self.results[mtype]['CI'][key] else: pars_ur[key] = 0.0 - + #-- get (approximated) reddened and unreddened model wave,flux = model.get_table(**pars) wave_ur,flux_ur = model.get_table(**pars_ur) @@ -3813,9 +3816,9 @@ def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004', **kwargs): #-- get synthetic photometry synflux_,pars = model.get_itable(photbands=self.master['photband'][keep], **pars) flux,flux_ur = flux*scale,flux_ur*scale - + synflux[keep] = synflux_ - + synflux[-self.master['color']] *= scale chi2 = (self.master['cmeas']-synflux)**2/self.master['e_cmeas']**2 #-- calculate effective wavelengths of the photometric bands via the model @@ -3824,7 +3827,7 @@ def set_best_model(self,mtype='igrid_search',law='fitzpatrick2004', **kwargs): self.results[mtype]['model'] = wave,flux,flux_ur self.results[mtype]['synflux'] = eff_waves,synflux,self.master['photband'] self.results[mtype]['chi2'] = chi2 - + class Calibrator(SED): """ Convenience class for a photometric standard star or calibrator. @@ -3844,7 +3847,7 @@ def __init__(self,ID=None,photfile=None,plx=None,load_fits=False,label='',librar photfile = phot_files[index] super(Calibrator,self).__init__(photfile=photfile,plx=plx,load_fits=load_fits,label=label) self.set_model(wave,flux) - + class SampleSEDs(object): """ @@ -3853,7 +3856,7 @@ class SampleSEDs(object): def __init__(self,targets,**kwargs): """ Initialize a sample. - + This can be done either with a list of IDs or with a list of SED instances. The SED must exist! That is, for each ID, there must be a phot or FITS file. @@ -3876,24 +3879,24 @@ def __init__(self,targets,**kwargs): self.seds.append(mysed) self.targets.append(mysed.ID) logger.info("Loaded {} SEDs".format(len(self.seds))) - + def __iter__(self): """ Allow iteration over the SED instances. """ for sed in self.seds: yield sed - + def __len__(self): """ The length of a SampleSEDs instance is the number SED instances. """ return len(self.seds) - + def __getitem__(self,key): """ Implements various ways to get individual seds. - + Allows integer indexing, slicing, indexing with integer and boolean arrays. """ #-- via slicing @@ -3917,7 +3920,7 @@ def __getitem__(self,key): #-- that's all I can come up with else: raise TypeError("Cannot use arrays of type {} for indexing".format(key.dtype)) - + def summarize(self): #-- collect the names of all the different sources sources = {} @@ -3963,11 +3966,11 @@ def summarize(self): for i,name in enumerate(output.dtype.names): output[name] = summary[i] return output - + def get_data(self,source,photband,label=None): """ Get all data on a particular passband from a particular source. - + If label is not None, synthetic flux from a model will be added to the columns. """ @@ -3997,7 +4000,7 @@ def get_data(self,source,photband,label=None): output = np.rec.fromrecords(records,names=dtypes.names) output = np.array(output,dtype=dtypes) return output - + def get_confidence_interval(self,parameter='teff',mtype='igrid_search'): values = np.zeros((len(self),3)) for i,sed in enumerate(self): @@ -4007,8 +4010,8 @@ def get_confidence_interval(self,parameter='teff',mtype='igrid_search'): values[i,2] = sed.results[mtype]['CI'][parameter+'_u'] sed.clear() return values - - + + if __name__ == "__main__": @@ -4016,7 +4019,7 @@ def get_confidence_interval(self,parameter='teff',mtype='igrid_search'): import doctest import pprint from ivs.aux import loggers - + if not sys.argv[1:]: doctest.testmod() pl.show() @@ -4038,7 +4041,7 @@ def get_confidence_interval(self,parameter='teff',mtype='igrid_search'): os.unlink(mysed.photfile) logger.info('Removed %s'%(mysed.photfile)) raise SystemExit - + #-- clean up if os.path.isfile('HD180642.fits'): os.remove('HD180642.fits') @@ -4115,7 +4118,7 @@ def get_confidence_interval(self,parameter='teff',mtype='igrid_search'): pl.ylim(pars[:,2].min(),pars[:,2].max()) pl.colorbar() pl.show() - + diff --git a/sed/extinctionmodels.py b/sed/extinctionmodels.py index 4f3cd9b8f..0d88f9db0 100644 --- a/sed/extinctionmodels.py +++ b/sed/extinctionmodels.py @@ -131,7 +131,7 @@ def findext(lng, lat, model='drimmel', distance=None, **kwargs): # Astrophysics (ISSN 0004-6361), vol. 258, no. 1, p. 104-111.) # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -def findext_arenou(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='Av'): +def findext_arenou(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs): """ Find the predicted V-band extinction (Av) according to the 3D model for galactic extinction of Arenou et al, "Modelling @@ -196,7 +196,7 @@ def findext_arenou(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='A av = alpha*rr0 + beta*rr0**2. + (distance-rr0)*gamma #-- Marshall is standard in Ak, but you can change this: - redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V']) + redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'],**kwargs) return av/redflux[0] @@ -719,7 +719,7 @@ def findext_marshall(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm= @type redlaw: str @param Rv: Av/E(B-V) (standard: 3.1) @type Rv: float - @return: The extinction in K-band + @return: The extinction in V-band (or norm) @rtype: float """ @@ -775,7 +775,7 @@ def findext_marshall(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm= #-- Marshall is standard in Ak, but you can change this: #redwave, redflux = get_law(redlaw,Rv=Rv,wave_units='micron',norm='Av', wave=array([0.54,2.22])) - redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.K']) + redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.K'],**kwargs) return ak/redflux[0] #} @@ -1038,8 +1038,9 @@ def findext_drimmel(lng, lat, distance=None, rescaling=True, else: out = (absdisk + abspir + avloc).flatten() - #-- Marshall is standard in Ak, but you can change this: - redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V']) + #-- Drimmel is standard in Av, but you can change this: + # In case of norm=Av, the conversion factor is just 1. + redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'],**kwargs) return(out/redflux[0]) @@ -1462,7 +1463,7 @@ def findext_schlegel(ll, bb, distance = None, redlaw='cardelli1989', Rv=3.1, nor av = ebv*Rv #-- Marshall is standard in Ak, but you can change this: - redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V']) + redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'],**kwargs) return av/redflux[0] diff --git a/sed/reddening.py b/sed/reddening.py index 4d16369fc..3f048ca0b 100644 --- a/sed/reddening.py +++ b/sed/reddening.py @@ -184,7 +184,7 @@ def get_law(name,norm='E(B-V)',wave_units='AA',photbands=None,**kwargs): #-- set the units of the wavelengths - if wave_units != 'AA' and photbands is not None: + if wave_units != 'AA': wave = conversions.convert('AA',wave_units,wave) return wave,mag @@ -263,7 +263,8 @@ def chiar2006(Rv=3.1,curve='ism',**kwargs): """ Extinction curve at infrared wavelengths from Chiar and Tielens (2006) - We return A(lambda)/E(B-V), by multiplying A(lambda)/Av with Rv. + We return A(lambda)/Av, by multiplying A(lambda)/Ak Ak/Av=0.09 (see Chiar + and Tielens (2006). This is only defined for Rv=3.1. If it is different, this will raise an AssertionError @@ -297,6 +298,49 @@ def chiar2006(Rv=3.1,curve='ism',**kwargs): return wavelengths*1e4,alam_aV +@memoized +def fitz2004chiar2006(Rv=3.1,curve='ism',**kwargs): + """ + Combined and extrapolated extinction curve Fitzpatrick 2004 and + from Chiar and Tielens (2006). + + We return A(lambda)/E(B-V), by multiplying A(lambda)/Av with Rv. + + This is only defined for Rv=3.1. If it is different, this will raise an + AssertionError + + Extra kwags are to catch unwanted keyword arguments. + + @param Rv: Rv + @type Rv: float + @param curve: extinction curve + @type curve: string (one of 'gc' or 'ism', galactic centre or local ISM) + @return: wavelengths (A), A(lambda)/Av + @rtype: (ndarray,ndarray) + """ + + if curve.lower() not in ['ism','gc']: + raise ValueError,'No Fitzpatrick2004/Chiar2006 curve available for %s.'\ + %(curve) + + fn = 'fitzpatrick2004_chiar2006%s_extrapol.dat'%curve + source = os.path.join(basename,fn) + + #-- check Rv + assert(Rv==3.1) + wavelengths,alam_ak = ascii.read2array(source).T + + #-- Convert to AA + wavelengths *= 1e4 + + #-- Convert from Ak normalization to Av normalization. + norm_reddening = model.synthetic_flux(wavelengths,alam_ak,\ + ['JOHNSON.V','JOHNSON.K']) + ak_to_av = norm_reddening[1]/norm_reddening[0] + alam_aV = alam_ak * ak_to_av + + return wavelengths,alam_aV + @memoized def fitzpatrick1999(Rv=3.1,**kwargs): diff --git a/sed/redlaws/fitzpatrick2004_chiar2006gc_extrapol.dat b/sed/redlaws/fitzpatrick2004_chiar2006gc_extrapol.dat new file mode 100644 index 000000000..e036afd4c --- /dev/null +++ b/sed/redlaws/fitzpatrick2004_chiar2006gc_extrapol.dat @@ -0,0 +1,1430 @@ +#-- wavelength (micron) A_lambda/A_k +1.149e-01 4.301e+01 +1.161e-01 4.225e+01 +1.173e-01 4.152e+01 +1.185e-01 4.080e+01 +1.197e-01 4.009e+01 +1.210e-01 3.942e+01 +1.223e-01 3.876e+01 +1.236e-01 3.811e+01 +1.249e-01 3.749e+01 +1.263e-01 3.689e+01 +1.277e-01 3.630e+01 +1.291e-01 3.574e+01 +1.306e-01 3.519e+01 +1.321e-01 3.467e+01 +1.337e-01 3.417e+01 +1.352e-01 3.368e+01 +1.368e-01 3.322e+01 +1.385e-01 3.278e+01 +1.402e-01 3.236e+01 +1.419e-01 3.196e+01 +1.437e-01 3.159e+01 +1.455e-01 3.124e+01 +1.474e-01 3.091e+01 +1.493e-01 3.060e+01 +1.512e-01 3.032e+01 +1.533e-01 3.006e+01 +1.553e-01 2.984e+01 +1.575e-01 2.964e+01 +1.596e-01 2.947e+01 +1.619e-01 2.933e+01 +1.642e-01 2.924e+01 +1.666e-01 2.917e+01 +1.690e-01 2.912e+01 +1.716e-01 2.911e+01 +1.742e-01 2.914e+01 +1.768e-01 2.922e+01 +1.796e-01 2.935e+01 +1.824e-01 2.957e+01 +1.854e-01 2.988e+01 +1.884e-01 3.031e+01 +1.916e-01 3.090e+01 +1.948e-01 3.169e+01 +1.982e-01 3.271e+01 +2.017e-01 3.397e+01 +2.053e-01 3.544e+01 +2.090e-01 3.694e+01 +2.129e-01 3.813e+01 +2.169e-01 3.855e+01 +2.210e-01 3.792e+01 +2.254e-01 3.634e+01 +2.299e-01 3.426e+01 +2.346e-01 3.211e+01 +2.395e-01 3.016e+01 +2.446e-01 2.851e+01 +2.499e-01 2.715e+01 +2.554e-01 2.603e+01 +2.612e-01 2.510e+01 +2.673e-01 2.432e+01 +2.737e-01 2.362e+01 +2.803e-01 2.291e+01 +2.874e-01 2.225e+01 +2.947e-01 2.160e+01 +3.025e-01 2.100e+01 +3.107e-01 2.042e+01 +3.193e-01 1.988e+01 +3.284e-01 1.938e+01 +3.381e-01 1.892e+01 +3.483e-01 1.850e+01 +3.592e-01 1.809e+01 +3.708e-01 1.768e+01 +3.831e-01 1.725e+01 +3.964e-01 1.678e+01 +4.105e-01 1.625e+01 +4.257e-01 1.567e+01 +4.421e-01 1.505e+01 +4.598e-01 1.439e+01 +4.789e-01 1.369e+01 +4.998e-01 1.297e+01 +5.225e-01 1.223e+01 +5.473e-01 1.146e+01 +5.747e-01 1.068e+01 +6.050e-01 9.889e+00 +6.386e-01 9.094e+00 +6.761e-01 8.296e+00 +7.184e-01 7.505e+00 +7.663e-01 6.725e+00 +8.210e-01 5.960e+00 +8.842e-01 5.221e+00 +9.579e-01 4.516e+00 +1.045e+00 3.848e+00 +1.149e+00 3.232e+00 +1.240e+00 2.799e+00 +1.340e+00 2.378e+00 +1.440e+00 2.057e+00 +1.540e+00 1.807e+00 +1.640e+00 1.608e+00 +1.740e+00 1.447e+00 +1.840e+00 1.314e+00 +1.940e+00 1.203e+00 +2.040e+00 1.111e+00 +2.140e+00 1.031e+00 +2.240e+00 9.637e-01 +2.340e+00 9.051e-01 +2.440e+00 8.548e-01 +2.540e+00 8.096e-01 +2.640e+00 7.705e-01 +2.740e+00 7.356e-01 +2.840e+00 7.048e-01 +2.940e+00 6.770e-01 +3.040e+00 6.514e-01 +3.140e+00 6.298e-01 +3.240e+00 6.092e-01 +3.340e+00 5.907e-01 +3.440e+00 5.733e-01 +3.540e+00 5.589e-01 +3.640e+00 5.445e-01 +3.740e+00 5.312e-01 +3.840e+00 5.199e-01 +3.940e+00 5.096e-01 +4.040e+00 4.993e-01 +4.140e+00 4.901e-01 +4.240e+00 4.818e-01 +4.340e+00 4.736e-01 +4.440e+00 4.664e-01 +4.540e+00 4.592e-01 +4.640e+00 4.531e-01 +4.740e+00 4.479e-01 +4.840e+00 4.428e-01 +4.940e+00 4.377e-01 +5.040e+00 4.325e-01 +5.140e+00 4.284e-01 +5.240e+00 4.243e-01 +5.340e+00 4.212e-01 +5.440e+00 4.171e-01 +5.540e+00 4.140e-01 +5.640e+00 4.110e-01 +5.740e+00 4.089e-01 +5.840e+00 4.058e-01 +5.940e+00 4.038e-01 +6.040e+00 4.017e-01 +6.140e+00 3.997e-01 +6.240e+00 3.976e-01 +6.340e+00 3.966e-01 +6.440e+00 3.945e-01 +6.540e+00 3.935e-01 +6.640e+00 3.925e-01 +6.740e+00 3.914e-01 +6.840e+00 3.904e-01 +6.940e+00 3.894e-01 +7.040e+00 3.884e-01 +7.140e+00 3.884e-01 +7.240e+00 3.873e-01 +7.340e+00 3.873e-01 +7.440e+00 3.863e-01 +7.540e+00 3.863e-01 +7.640e+00 3.863e-01 +7.740e+00 3.863e-01 +7.840e+00 3.863e-01 +7.940e+00 3.863e-01 +8.100e+00 4.243e-01 +8.200e+00 4.592e-01 +8.300e+00 5.075e-01 +8.400e+00 5.671e-01 +8.500e+00 6.390e-01 +8.600e+00 7.181e-01 +8.700e+00 7.972e-01 +8.800e+00 8.712e-01 +8.900e+00 9.339e-01 +9.000e+00 9.986e-01 +9.100e+00 1.087e+00 +9.200e+00 1.198e+00 +9.300e+00 1.287e+00 +9.400e+00 1.353e+00 +9.500e+00 1.403e+00 +9.600e+00 1.429e+00 +9.700e+00 1.423e+00 +9.800e+00 1.391e+00 +9.900e+00 1.363e+00 +1.000e+01 1.341e+00 +1.010e+01 1.309e+00 +1.020e+01 1.262e+00 +1.030e+01 1.208e+00 +1.040e+01 1.157e+00 +1.050e+01 1.114e+00 +1.060e+01 1.078e+00 +1.070e+01 1.045e+00 +1.080e+01 1.014e+00 +1.090e+01 9.801e-01 +1.100e+01 9.421e-01 +1.110e+01 9.010e-01 +1.120e+01 8.579e-01 +1.130e+01 8.147e-01 +1.140e+01 7.746e-01 +1.150e+01 7.366e-01 +1.160e+01 7.007e-01 +1.170e+01 6.668e-01 +1.180e+01 6.329e-01 +1.190e+01 6.000e-01 +1.200e+01 5.671e-01 +1.210e+01 5.363e-01 +1.220e+01 5.096e-01 +1.230e+01 4.870e-01 +1.240e+01 4.705e-01 +1.250e+01 4.592e-01 +1.260e+01 4.531e-01 +1.270e+01 4.500e-01 +1.280e+01 4.469e-01 +1.290e+01 4.449e-01 +1.300e+01 4.428e-01 +1.310e+01 4.407e-01 +1.320e+01 4.387e-01 +1.330e+01 4.366e-01 +1.340e+01 4.356e-01 +1.350e+01 4.356e-01 +1.360e+01 4.356e-01 +1.370e+01 4.356e-01 +1.380e+01 4.356e-01 +1.390e+01 4.366e-01 +1.400e+01 4.377e-01 +1.410e+01 4.387e-01 +1.420e+01 4.397e-01 +1.430e+01 4.418e-01 +1.440e+01 4.438e-01 +1.450e+01 4.459e-01 +1.460e+01 4.500e-01 +1.470e+01 4.551e-01 +1.480e+01 4.623e-01 +1.490e+01 4.705e-01 +1.500e+01 4.798e-01 +1.510e+01 4.901e-01 +1.520e+01 5.014e-01 +1.530e+01 5.137e-01 +1.540e+01 5.260e-01 +1.550e+01 5.394e-01 +1.560e+01 5.527e-01 +1.570e+01 5.661e-01 +1.580e+01 5.805e-01 +1.590e+01 5.938e-01 +1.600e+01 6.062e-01 +1.610e+01 6.195e-01 +1.620e+01 6.318e-01 +1.630e+01 6.431e-01 +1.640e+01 6.555e-01 +1.650e+01 6.657e-01 +1.660e+01 6.770e-01 +1.670e+01 6.873e-01 +1.680e+01 6.966e-01 +1.690e+01 7.058e-01 +1.700e+01 7.140e-01 +1.710e+01 7.222e-01 +1.720e+01 7.294e-01 +1.730e+01 7.366e-01 +1.740e+01 7.438e-01 +1.750e+01 7.500e-01 +1.760e+01 7.562e-01 +1.770e+01 7.613e-01 +1.780e+01 7.664e-01 +1.790e+01 7.705e-01 +1.800e+01 7.746e-01 +1.810e+01 7.777e-01 +1.820e+01 7.808e-01 +1.830e+01 7.829e-01 +1.840e+01 7.839e-01 +1.850e+01 7.839e-01 +1.860e+01 7.839e-01 +1.870e+01 7.829e-01 +1.880e+01 7.808e-01 +1.890e+01 7.788e-01 +1.900e+01 7.767e-01 +1.910e+01 7.736e-01 +1.920e+01 7.695e-01 +1.930e+01 7.654e-01 +1.940e+01 7.613e-01 +1.950e+01 7.562e-01 +1.960e+01 7.520e-01 +1.970e+01 7.469e-01 +1.980e+01 7.418e-01 +1.990e+01 7.366e-01 +2.000e+01 7.315e-01 +2.010e+01 7.264e-01 +2.020e+01 7.212e-01 +2.030e+01 7.171e-01 +2.040e+01 7.120e-01 +2.050e+01 7.079e-01 +2.060e+01 7.048e-01 +2.070e+01 7.007e-01 +2.080e+01 6.976e-01 +2.090e+01 6.945e-01 +2.100e+01 6.914e-01 +2.110e+01 6.883e-01 +2.120e+01 6.853e-01 +2.130e+01 6.822e-01 +2.140e+01 6.791e-01 +2.150e+01 6.760e-01 +2.160e+01 6.740e-01 +2.170e+01 6.709e-01 +2.180e+01 6.678e-01 +2.190e+01 6.657e-01 +2.200e+01 6.627e-01 +2.210e+01 6.596e-01 +2.220e+01 6.575e-01 +2.230e+01 6.544e-01 +2.240e+01 6.524e-01 +2.250e+01 6.493e-01 +2.260e+01 6.473e-01 +2.270e+01 6.442e-01 +2.280e+01 6.421e-01 +2.290e+01 6.401e-01 +2.300e+01 6.380e-01 +2.310e+01 6.359e-01 +2.320e+01 6.339e-01 +2.330e+01 6.318e-01 +2.340e+01 6.298e-01 +2.350e+01 6.277e-01 +2.360e+01 6.257e-01 +2.370e+01 6.236e-01 +2.380e+01 6.216e-01 +2.390e+01 6.195e-01 +2.400e+01 6.175e-01 +2.410e+01 6.154e-01 +2.420e+01 6.133e-01 +2.430e+01 6.113e-01 +2.440e+01 6.092e-01 +2.450e+01 6.072e-01 +2.460e+01 6.051e-01 +2.470e+01 6.031e-01 +2.480e+01 6.010e-01 +2.490e+01 5.990e-01 +2.500e+01 5.969e-01 +2.510e+01 5.949e-01 +2.520e+01 5.928e-01 +2.530e+01 5.907e-01 +2.540e+01 5.887e-01 +2.550e+01 5.866e-01 +2.560e+01 5.836e-01 +2.570e+01 5.815e-01 +2.580e+01 5.794e-01 +2.590e+01 5.774e-01 +2.600e+01 5.753e-01 +2.610e+01 5.733e-01 +2.620e+01 5.712e-01 +2.630e+01 5.681e-01 +2.640e+01 5.661e-01 +2.650e+01 5.640e-01 +2.660e+01 5.620e-01 +2.670e+01 5.599e-01 +2.680e+01 5.579e-01 +2.690e+01 5.548e-01 +2.700e+01 5.527e-01 +2.710e+01 5.507e-01 +2.720e+01 5.486e-01 +2.730e+01 5.466e-01 +2.740e+01 5.445e-01 +2.750e+01 5.414e-01 +2.760e+01 5.394e-01 +2.770e+01 5.373e-01 +2.780e+01 5.353e-01 +2.790e+01 5.332e-01 +2.800e+01 5.312e-01 +2.810e+01 5.291e-01 +2.820e+01 5.270e-01 +2.830e+01 5.250e-01 +2.840e+01 5.219e-01 +2.850e+01 5.199e-01 +2.860e+01 5.178e-01 +2.870e+01 5.157e-01 +2.880e+01 5.137e-01 +2.890e+01 5.116e-01 +2.900e+01 5.096e-01 +2.910e+01 5.075e-01 +2.920e+01 5.055e-01 +2.930e+01 5.034e-01 +2.940e+01 5.014e-01 +2.950e+01 5.003e-01 +2.960e+01 4.983e-01 +2.970e+01 4.962e-01 +2.980e+01 4.942e-01 +2.990e+01 4.921e-01 +3.000e+01 4.901e-01 +3.010e+01 4.890e-01 +3.020e+01 4.870e-01 +3.030e+01 4.849e-01 +3.040e+01 4.829e-01 +3.050e+01 4.818e-01 +3.060e+01 4.798e-01 +3.070e+01 4.777e-01 +3.080e+01 4.767e-01 +3.090e+01 4.747e-01 +3.100e+01 4.736e-01 +3.110e+01 4.716e-01 +3.120e+01 4.705e-01 +3.130e+01 4.695e-01 +3.140e+01 4.675e-01 +3.150e+01 4.664e-01 +3.160e+01 4.644e-01 +3.170e+01 4.633e-01 +3.180e+01 4.623e-01 +3.190e+01 4.613e-01 +3.200e+01 4.603e-01 +3.210e+01 4.582e-01 +3.220e+01 4.572e-01 +3.230e+01 4.562e-01 +3.240e+01 4.551e-01 +3.250e+01 4.541e-01 +3.260e+01 4.531e-01 +3.270e+01 4.520e-01 +3.280e+01 4.510e-01 +3.290e+01 4.500e-01 +3.300e+01 4.490e-01 +3.310e+01 4.479e-01 +3.320e+01 4.469e-01 +3.330e+01 4.459e-01 +3.340e+01 4.449e-01 +3.350e+01 4.449e-01 +3.360e+01 4.438e-01 +3.370e+01 4.428e-01 +3.380e+01 4.418e-01 +3.390e+01 4.418e-01 +3.400e+01 4.407e-01 +3.410e+01 4.397e-01 +3.420e+01 4.397e-01 +3.430e+01 4.387e-01 +3.440e+01 4.387e-01 +3.450e+01 4.387e-01 +3.460e+01 4.387e-01 +3.470e+01 4.377e-01 +3.480e+01 4.377e-01 +3.490e+01 4.377e-01 +3.500e+01 4.377e-01 +3.510e+01 4.354e-01 +3.607e+01 4.147e-01 +3.703e+01 3.954e-01 +3.800e+01 3.775e-01 +3.896e+01 3.608e-01 +3.993e+01 3.453e-01 +4.090e+01 3.307e-01 +4.186e+01 3.171e-01 +4.283e+01 3.044e-01 +4.379e+01 2.924e-01 +4.476e+01 2.811e-01 +4.572e+01 2.705e-01 +4.669e+01 2.605e-01 +4.766e+01 2.511e-01 +4.862e+01 2.422e-01 +4.959e+01 2.338e-01 +5.055e+01 2.258e-01 +5.152e+01 2.182e-01 +5.249e+01 2.111e-01 +5.345e+01 2.042e-01 +5.442e+01 1.978e-01 +5.538e+01 1.916e-01 +5.635e+01 1.857e-01 +5.731e+01 1.801e-01 +5.828e+01 1.748e-01 +5.925e+01 1.697e-01 +6.021e+01 1.648e-01 +6.118e+01 1.602e-01 +6.214e+01 1.557e-01 +6.311e+01 1.515e-01 +6.408e+01 1.474e-01 +6.504e+01 1.435e-01 +6.601e+01 1.397e-01 +6.697e+01 1.361e-01 +6.794e+01 1.326e-01 +6.891e+01 1.293e-01 +6.987e+01 1.261e-01 +7.084e+01 1.230e-01 +7.180e+01 1.201e-01 +7.277e+01 1.172e-01 +7.373e+01 1.145e-01 +7.470e+01 1.118e-01 +7.567e+01 1.093e-01 +7.663e+01 1.068e-01 +7.760e+01 1.044e-01 +7.856e+01 1.021e-01 +7.953e+01 9.989e-02 +8.050e+01 9.774e-02 +8.146e+01 9.566e-02 +8.243e+01 9.366e-02 +8.339e+01 9.171e-02 +8.436e+01 8.983e-02 +8.533e+01 8.801e-02 +8.629e+01 8.624e-02 +8.726e+01 8.453e-02 +8.822e+01 8.287e-02 +8.919e+01 8.127e-02 +9.015e+01 7.971e-02 +9.112e+01 7.819e-02 +9.209e+01 7.672e-02 +9.305e+01 7.529e-02 +9.402e+01 7.391e-02 +9.498e+01 7.256e-02 +9.595e+01 7.125e-02 +9.692e+01 6.998e-02 +9.788e+01 6.874e-02 +9.885e+01 6.753e-02 +9.981e+01 6.636e-02 +1.008e+02 6.522e-02 +1.017e+02 6.411e-02 +1.027e+02 6.303e-02 +1.037e+02 6.198e-02 +1.046e+02 6.095e-02 +1.056e+02 5.995e-02 +1.066e+02 5.898e-02 +1.075e+02 5.803e-02 +1.085e+02 5.710e-02 +1.095e+02 5.620e-02 +1.104e+02 5.532e-02 +1.114e+02 5.446e-02 +1.124e+02 5.362e-02 +1.133e+02 5.280e-02 +1.143e+02 5.200e-02 +1.153e+02 5.121e-02 +1.162e+02 5.045e-02 +1.172e+02 4.971e-02 +1.182e+02 4.898e-02 +1.191e+02 4.826e-02 +1.201e+02 4.757e-02 +1.211e+02 4.689e-02 +1.220e+02 4.622e-02 +1.230e+02 4.557e-02 +1.240e+02 4.493e-02 +1.249e+02 4.431e-02 +1.259e+02 4.370e-02 +1.269e+02 4.310e-02 +1.278e+02 4.252e-02 +1.288e+02 4.195e-02 +1.298e+02 4.138e-02 +1.307e+02 4.084e-02 +1.317e+02 4.030e-02 +1.327e+02 3.977e-02 +1.336e+02 3.926e-02 +1.346e+02 3.875e-02 +1.356e+02 3.825e-02 +1.365e+02 3.777e-02 +1.375e+02 3.729e-02 +1.384e+02 3.683e-02 +1.394e+02 3.637e-02 +1.404e+02 3.592e-02 +1.413e+02 3.548e-02 +1.423e+02 3.505e-02 +1.433e+02 3.462e-02 +1.442e+02 3.421e-02 +1.452e+02 3.380e-02 +1.462e+02 3.340e-02 +1.471e+02 3.300e-02 +1.481e+02 3.262e-02 +1.491e+02 3.224e-02 +1.500e+02 3.186e-02 +1.510e+02 3.150e-02 +1.520e+02 3.114e-02 +1.529e+02 3.079e-02 +1.539e+02 3.044e-02 +1.549e+02 3.010e-02 +1.558e+02 2.976e-02 +1.568e+02 2.943e-02 +1.578e+02 2.911e-02 +1.587e+02 2.879e-02 +1.597e+02 2.848e-02 +1.607e+02 2.817e-02 +1.616e+02 2.787e-02 +1.626e+02 2.757e-02 +1.636e+02 2.728e-02 +1.645e+02 2.699e-02 +1.655e+02 2.671e-02 +1.665e+02 2.643e-02 +1.674e+02 2.616e-02 +1.684e+02 2.589e-02 +1.694e+02 2.562e-02 +1.703e+02 2.536e-02 +1.713e+02 2.510e-02 +1.723e+02 2.485e-02 +1.732e+02 2.460e-02 +1.742e+02 2.436e-02 +1.752e+02 2.412e-02 +1.761e+02 2.388e-02 +1.771e+02 2.365e-02 +1.780e+02 2.342e-02 +1.790e+02 2.319e-02 +1.800e+02 2.296e-02 +1.809e+02 2.274e-02 +1.819e+02 2.253e-02 +1.829e+02 2.231e-02 +1.838e+02 2.210e-02 +1.848e+02 2.190e-02 +1.858e+02 2.169e-02 +1.867e+02 2.149e-02 +1.877e+02 2.129e-02 +1.887e+02 2.110e-02 +1.896e+02 2.090e-02 +1.906e+02 2.071e-02 +1.916e+02 2.052e-02 +1.925e+02 2.034e-02 +1.935e+02 2.016e-02 +1.945e+02 1.998e-02 +1.954e+02 1.980e-02 +1.964e+02 1.963e-02 +1.974e+02 1.945e-02 +1.983e+02 1.928e-02 +1.993e+02 1.911e-02 +2.003e+02 1.895e-02 +2.012e+02 1.879e-02 +2.022e+02 1.862e-02 +2.032e+02 1.847e-02 +2.041e+02 1.831e-02 +2.051e+02 1.815e-02 +2.061e+02 1.800e-02 +2.070e+02 1.785e-02 +2.080e+02 1.770e-02 +2.090e+02 1.755e-02 +2.099e+02 1.741e-02 +2.109e+02 1.727e-02 +2.119e+02 1.712e-02 +2.128e+02 1.698e-02 +2.138e+02 1.685e-02 +2.148e+02 1.671e-02 +2.157e+02 1.658e-02 +2.167e+02 1.644e-02 +2.176e+02 1.631e-02 +2.186e+02 1.618e-02 +2.196e+02 1.605e-02 +2.205e+02 1.593e-02 +2.215e+02 1.580e-02 +2.225e+02 1.568e-02 +2.234e+02 1.556e-02 +2.244e+02 1.544e-02 +2.254e+02 1.532e-02 +2.263e+02 1.520e-02 +2.273e+02 1.509e-02 +2.283e+02 1.497e-02 +2.292e+02 1.486e-02 +2.302e+02 1.475e-02 +2.312e+02 1.463e-02 +2.321e+02 1.453e-02 +2.331e+02 1.442e-02 +2.341e+02 1.431e-02 +2.350e+02 1.420e-02 +2.360e+02 1.410e-02 +2.370e+02 1.400e-02 +2.379e+02 1.389e-02 +2.389e+02 1.379e-02 +2.399e+02 1.369e-02 +2.408e+02 1.360e-02 +2.418e+02 1.350e-02 +2.428e+02 1.340e-02 +2.437e+02 1.331e-02 +2.447e+02 1.321e-02 +2.457e+02 1.312e-02 +2.466e+02 1.303e-02 +2.476e+02 1.293e-02 +2.486e+02 1.284e-02 +2.495e+02 1.275e-02 +2.505e+02 1.267e-02 +2.515e+02 1.258e-02 +2.524e+02 1.249e-02 +2.534e+02 1.241e-02 +2.544e+02 1.232e-02 +2.553e+02 1.224e-02 +2.563e+02 1.216e-02 +2.572e+02 1.207e-02 +2.582e+02 1.199e-02 +2.592e+02 1.191e-02 +2.601e+02 1.183e-02 +2.611e+02 1.175e-02 +2.621e+02 1.168e-02 +2.630e+02 1.160e-02 +2.640e+02 1.152e-02 +2.650e+02 1.145e-02 +2.659e+02 1.137e-02 +2.669e+02 1.130e-02 +2.679e+02 1.123e-02 +2.688e+02 1.115e-02 +2.698e+02 1.108e-02 +2.708e+02 1.101e-02 +2.717e+02 1.094e-02 +2.727e+02 1.087e-02 +2.737e+02 1.080e-02 +2.746e+02 1.073e-02 +2.756e+02 1.066e-02 +2.766e+02 1.060e-02 +2.775e+02 1.053e-02 +2.785e+02 1.047e-02 +2.795e+02 1.040e-02 +2.804e+02 1.034e-02 +2.814e+02 1.027e-02 +2.824e+02 1.021e-02 +2.833e+02 1.015e-02 +2.843e+02 1.009e-02 +2.853e+02 1.002e-02 +2.862e+02 9.963e-03 +2.872e+02 9.903e-03 +2.882e+02 9.843e-03 +2.891e+02 9.784e-03 +2.901e+02 9.725e-03 +2.911e+02 9.667e-03 +2.920e+02 9.610e-03 +2.930e+02 9.553e-03 +2.940e+02 9.497e-03 +2.949e+02 9.441e-03 +2.959e+02 9.385e-03 +2.968e+02 9.330e-03 +2.978e+02 9.276e-03 +2.988e+02 9.222e-03 +2.997e+02 9.169e-03 +3.007e+02 9.116e-03 +3.017e+02 9.063e-03 +3.026e+02 9.011e-03 +3.036e+02 8.960e-03 +3.046e+02 8.909e-03 +3.055e+02 8.858e-03 +3.065e+02 8.808e-03 +3.075e+02 8.758e-03 +3.084e+02 8.709e-03 +3.094e+02 8.660e-03 +3.104e+02 8.611e-03 +3.113e+02 8.563e-03 +3.123e+02 8.516e-03 +3.133e+02 8.469e-03 +3.142e+02 8.422e-03 +3.152e+02 8.375e-03 +3.162e+02 8.329e-03 +3.171e+02 8.284e-03 +3.181e+02 8.239e-03 +3.191e+02 8.194e-03 +3.200e+02 8.149e-03 +3.210e+02 8.105e-03 +3.220e+02 8.061e-03 +3.229e+02 8.018e-03 +3.239e+02 7.975e-03 +3.249e+02 7.933e-03 +3.258e+02 7.890e-03 +3.268e+02 7.848e-03 +3.278e+02 7.807e-03 +3.287e+02 7.765e-03 +3.297e+02 7.725e-03 +3.307e+02 7.684e-03 +3.316e+02 7.644e-03 +3.326e+02 7.604e-03 +3.336e+02 7.564e-03 +3.345e+02 7.525e-03 +3.355e+02 7.486e-03 +3.365e+02 7.447e-03 +3.374e+02 7.409e-03 +3.384e+02 7.371e-03 +3.393e+02 7.333e-03 +3.403e+02 7.296e-03 +3.413e+02 7.259e-03 +3.422e+02 7.222e-03 +3.432e+02 7.185e-03 +3.442e+02 7.149e-03 +3.451e+02 7.113e-03 +3.461e+02 7.078e-03 +3.471e+02 7.042e-03 +3.480e+02 7.007e-03 +3.490e+02 6.972e-03 +3.500e+02 6.938e-03 +3.509e+02 6.903e-03 +3.519e+02 6.869e-03 +3.529e+02 6.835e-03 +3.538e+02 6.802e-03 +3.548e+02 6.768e-03 +3.558e+02 6.735e-03 +3.567e+02 6.703e-03 +3.577e+02 6.670e-03 +3.587e+02 6.638e-03 +3.596e+02 6.606e-03 +3.606e+02 6.574e-03 +3.616e+02 6.542e-03 +3.625e+02 6.511e-03 +3.635e+02 6.480e-03 +3.645e+02 6.449e-03 +3.654e+02 6.418e-03 +3.664e+02 6.388e-03 +3.674e+02 6.358e-03 +3.683e+02 6.328e-03 +3.693e+02 6.298e-03 +3.703e+02 6.268e-03 +3.712e+02 6.239e-03 +3.722e+02 6.210e-03 +3.732e+02 6.181e-03 +3.741e+02 6.152e-03 +3.751e+02 6.124e-03 +3.761e+02 6.096e-03 +3.770e+02 6.068e-03 +3.780e+02 6.040e-03 +3.789e+02 6.012e-03 +3.799e+02 5.985e-03 +3.809e+02 5.957e-03 +3.818e+02 5.930e-03 +3.828e+02 5.903e-03 +3.838e+02 5.877e-03 +3.847e+02 5.850e-03 +3.857e+02 5.824e-03 +3.867e+02 5.797e-03 +3.876e+02 5.772e-03 +3.886e+02 5.746e-03 +3.896e+02 5.720e-03 +3.905e+02 5.695e-03 +3.915e+02 5.669e-03 +3.925e+02 5.644e-03 +3.934e+02 5.619e-03 +3.944e+02 5.595e-03 +3.954e+02 5.570e-03 +3.963e+02 5.546e-03 +3.973e+02 5.521e-03 +3.983e+02 5.497e-03 +3.992e+02 5.473e-03 +4.002e+02 5.450e-03 +4.012e+02 5.426e-03 +4.021e+02 5.403e-03 +4.031e+02 5.379e-03 +4.041e+02 5.356e-03 +4.050e+02 5.333e-03 +4.060e+02 5.310e-03 +4.070e+02 5.288e-03 +4.079e+02 5.265e-03 +4.089e+02 5.243e-03 +4.099e+02 5.221e-03 +4.108e+02 5.199e-03 +4.118e+02 5.177e-03 +4.128e+02 5.155e-03 +4.137e+02 5.133e-03 +4.147e+02 5.112e-03 +4.157e+02 5.090e-03 +4.166e+02 5.069e-03 +4.176e+02 5.048e-03 +4.185e+02 5.027e-03 +4.195e+02 5.006e-03 +4.205e+02 4.986e-03 +4.214e+02 4.965e-03 +4.224e+02 4.945e-03 +4.234e+02 4.924e-03 +4.243e+02 4.904e-03 +4.253e+02 4.884e-03 +4.263e+02 4.864e-03 +4.272e+02 4.845e-03 +4.282e+02 4.825e-03 +4.292e+02 4.805e-03 +4.301e+02 4.786e-03 +4.311e+02 4.767e-03 +4.321e+02 4.748e-03 +4.330e+02 4.728e-03 +4.340e+02 4.710e-03 +4.350e+02 4.691e-03 +4.359e+02 4.672e-03 +4.369e+02 4.653e-03 +4.379e+02 4.635e-03 +4.388e+02 4.617e-03 +4.398e+02 4.598e-03 +4.408e+02 4.580e-03 +4.417e+02 4.562e-03 +4.427e+02 4.544e-03 +4.437e+02 4.527e-03 +4.446e+02 4.509e-03 +4.456e+02 4.491e-03 +4.466e+02 4.474e-03 +4.475e+02 4.456e-03 +4.485e+02 4.439e-03 +4.495e+02 4.422e-03 +4.504e+02 4.405e-03 +4.514e+02 4.388e-03 +4.524e+02 4.371e-03 +4.533e+02 4.354e-03 +4.543e+02 4.338e-03 +4.553e+02 4.321e-03 +4.562e+02 4.305e-03 +4.572e+02 4.288e-03 +4.581e+02 4.272e-03 +4.591e+02 4.256e-03 +4.601e+02 4.240e-03 +4.610e+02 4.224e-03 +4.620e+02 4.208e-03 +4.630e+02 4.192e-03 +4.639e+02 4.177e-03 +4.649e+02 4.161e-03 +4.659e+02 4.146e-03 +4.668e+02 4.130e-03 +4.678e+02 4.115e-03 +4.688e+02 4.099e-03 +4.697e+02 4.084e-03 +4.707e+02 4.069e-03 +4.717e+02 4.054e-03 +4.726e+02 4.039e-03 +4.736e+02 4.025e-03 +4.746e+02 4.010e-03 +4.755e+02 3.995e-03 +4.765e+02 3.981e-03 +4.775e+02 3.966e-03 +4.784e+02 3.952e-03 +4.794e+02 3.937e-03 +4.804e+02 3.923e-03 +4.813e+02 3.909e-03 +4.823e+02 3.895e-03 +4.833e+02 3.881e-03 +4.842e+02 3.867e-03 +4.852e+02 3.853e-03 +4.862e+02 3.839e-03 +4.871e+02 3.826e-03 +4.881e+02 3.812e-03 +4.891e+02 3.799e-03 +4.900e+02 3.785e-03 +4.910e+02 3.772e-03 +4.920e+02 3.758e-03 +4.929e+02 3.745e-03 +4.939e+02 3.732e-03 +4.949e+02 3.719e-03 +4.958e+02 3.706e-03 +4.968e+02 3.693e-03 +4.977e+02 3.680e-03 +4.987e+02 3.667e-03 +4.997e+02 3.654e-03 +5.006e+02 3.642e-03 +5.016e+02 3.629e-03 +5.026e+02 3.617e-03 +5.035e+02 3.604e-03 +5.045e+02 3.592e-03 +5.055e+02 3.579e-03 +5.064e+02 3.567e-03 +5.074e+02 3.555e-03 +5.084e+02 3.543e-03 +5.093e+02 3.531e-03 +5.103e+02 3.519e-03 +5.113e+02 3.507e-03 +5.122e+02 3.495e-03 +5.132e+02 3.483e-03 +5.142e+02 3.471e-03 +5.151e+02 3.459e-03 +5.161e+02 3.448e-03 +5.171e+02 3.436e-03 +5.180e+02 3.425e-03 +5.190e+02 3.413e-03 +5.200e+02 3.402e-03 +5.209e+02 3.390e-03 +5.219e+02 3.379e-03 +5.229e+02 3.368e-03 +5.238e+02 3.357e-03 +5.248e+02 3.346e-03 +5.258e+02 3.335e-03 +5.267e+02 3.324e-03 +5.277e+02 3.313e-03 +5.287e+02 3.302e-03 +5.296e+02 3.291e-03 +5.306e+02 3.280e-03 +5.316e+02 3.269e-03 +5.325e+02 3.259e-03 +5.335e+02 3.248e-03 +5.345e+02 3.238e-03 +5.354e+02 3.227e-03 +5.364e+02 3.217e-03 +5.374e+02 3.206e-03 +5.383e+02 3.196e-03 +5.393e+02 3.186e-03 +5.402e+02 3.175e-03 +5.412e+02 3.165e-03 +5.422e+02 3.155e-03 +5.431e+02 3.145e-03 +5.441e+02 3.135e-03 +5.451e+02 3.125e-03 +5.460e+02 3.115e-03 +5.470e+02 3.105e-03 +5.480e+02 3.095e-03 +5.489e+02 3.085e-03 +5.499e+02 3.076e-03 +5.509e+02 3.066e-03 +5.518e+02 3.056e-03 +5.528e+02 3.047e-03 +5.538e+02 3.037e-03 +5.547e+02 3.028e-03 +5.557e+02 3.018e-03 +5.567e+02 3.009e-03 +5.576e+02 2.999e-03 +5.586e+02 2.990e-03 +5.596e+02 2.981e-03 +5.605e+02 2.972e-03 +5.615e+02 2.962e-03 +5.625e+02 2.953e-03 +5.634e+02 2.944e-03 +5.644e+02 2.935e-03 +5.654e+02 2.926e-03 +5.663e+02 2.917e-03 +5.673e+02 2.908e-03 +5.683e+02 2.899e-03 +5.692e+02 2.890e-03 +5.702e+02 2.882e-03 +5.712e+02 2.873e-03 +5.721e+02 2.864e-03 +5.731e+02 2.855e-03 +5.741e+02 2.847e-03 +5.750e+02 2.838e-03 +5.760e+02 2.830e-03 +5.770e+02 2.821e-03 +5.779e+02 2.813e-03 +5.789e+02 2.804e-03 +5.798e+02 2.796e-03 +5.808e+02 2.787e-03 +5.818e+02 2.779e-03 +5.827e+02 2.771e-03 +5.837e+02 2.763e-03 +5.847e+02 2.754e-03 +5.856e+02 2.746e-03 +5.866e+02 2.738e-03 +5.876e+02 2.730e-03 +5.885e+02 2.722e-03 +5.895e+02 2.714e-03 +5.905e+02 2.706e-03 +5.914e+02 2.698e-03 +5.924e+02 2.690e-03 +5.934e+02 2.682e-03 +5.943e+02 2.674e-03 +5.953e+02 2.666e-03 +5.963e+02 2.659e-03 +5.972e+02 2.651e-03 +5.982e+02 2.643e-03 +5.992e+02 2.636e-03 +6.001e+02 2.628e-03 +6.011e+02 2.620e-03 +6.021e+02 2.613e-03 +6.030e+02 2.605e-03 +6.040e+02 2.598e-03 +6.050e+02 2.590e-03 +6.059e+02 2.583e-03 +6.069e+02 2.576e-03 +6.079e+02 2.568e-03 +6.088e+02 2.561e-03 +6.098e+02 2.554e-03 +6.108e+02 2.546e-03 +6.117e+02 2.539e-03 +6.127e+02 2.532e-03 +6.137e+02 2.525e-03 +6.146e+02 2.518e-03 +6.156e+02 2.510e-03 +6.166e+02 2.503e-03 +6.175e+02 2.496e-03 +6.185e+02 2.489e-03 +6.194e+02 2.482e-03 +6.204e+02 2.475e-03 +6.214e+02 2.468e-03 +6.223e+02 2.462e-03 +6.233e+02 2.455e-03 +6.243e+02 2.448e-03 +6.252e+02 2.441e-03 +6.262e+02 2.434e-03 +6.272e+02 2.428e-03 +6.281e+02 2.421e-03 +6.291e+02 2.414e-03 +6.301e+02 2.407e-03 +6.310e+02 2.401e-03 +6.320e+02 2.394e-03 +6.330e+02 2.388e-03 +6.339e+02 2.381e-03 +6.349e+02 2.375e-03 +6.359e+02 2.368e-03 +6.368e+02 2.362e-03 +6.378e+02 2.355e-03 +6.388e+02 2.349e-03 +6.397e+02 2.342e-03 +6.407e+02 2.336e-03 +6.417e+02 2.330e-03 +6.426e+02 2.323e-03 +6.436e+02 2.317e-03 +6.446e+02 2.311e-03 +6.455e+02 2.305e-03 +6.465e+02 2.299e-03 +6.475e+02 2.292e-03 +6.484e+02 2.286e-03 +6.494e+02 2.280e-03 +6.504e+02 2.274e-03 +6.513e+02 2.268e-03 +6.523e+02 2.262e-03 +6.533e+02 2.256e-03 +6.542e+02 2.250e-03 +6.552e+02 2.244e-03 +6.562e+02 2.238e-03 +6.571e+02 2.232e-03 +6.581e+02 2.226e-03 +6.590e+02 2.220e-03 +6.600e+02 2.214e-03 +6.610e+02 2.209e-03 +6.619e+02 2.203e-03 +6.629e+02 2.197e-03 +6.639e+02 2.191e-03 +6.648e+02 2.186e-03 +6.658e+02 2.180e-03 +6.668e+02 2.174e-03 +6.677e+02 2.169e-03 +6.687e+02 2.163e-03 +6.697e+02 2.157e-03 +6.706e+02 2.152e-03 +6.716e+02 2.146e-03 +6.726e+02 2.141e-03 +6.735e+02 2.135e-03 +6.745e+02 2.130e-03 +6.755e+02 2.124e-03 +6.764e+02 2.119e-03 +6.774e+02 2.113e-03 +6.784e+02 2.108e-03 +6.793e+02 2.102e-03 +6.803e+02 2.097e-03 +6.813e+02 2.092e-03 +6.822e+02 2.086e-03 +6.832e+02 2.081e-03 +6.842e+02 2.076e-03 +6.851e+02 2.070e-03 +6.861e+02 2.065e-03 +6.871e+02 2.060e-03 +6.880e+02 2.055e-03 +6.890e+02 2.050e-03 +6.900e+02 2.044e-03 +6.909e+02 2.039e-03 +6.919e+02 2.034e-03 +6.929e+02 2.029e-03 +6.938e+02 2.024e-03 +6.948e+02 2.019e-03 +6.958e+02 2.014e-03 +6.967e+02 2.009e-03 +6.977e+02 2.004e-03 +6.986e+02 1.999e-03 +6.996e+02 1.994e-03 +7.006e+02 1.989e-03 +7.015e+02 1.984e-03 +7.025e+02 1.979e-03 +7.035e+02 1.974e-03 +7.044e+02 1.969e-03 +7.054e+02 1.965e-03 +7.064e+02 1.960e-03 +7.073e+02 1.955e-03 +7.083e+02 1.950e-03 +7.093e+02 1.945e-03 +7.102e+02 1.941e-03 +7.112e+02 1.936e-03 +7.122e+02 1.931e-03 +7.131e+02 1.926e-03 +7.141e+02 1.922e-03 +7.151e+02 1.917e-03 +7.160e+02 1.912e-03 +7.170e+02 1.908e-03 +7.180e+02 1.903e-03 +7.189e+02 1.899e-03 +7.199e+02 1.894e-03 +7.209e+02 1.889e-03 +7.218e+02 1.885e-03 +7.228e+02 1.880e-03 +7.238e+02 1.876e-03 +7.247e+02 1.871e-03 +7.257e+02 1.867e-03 +7.267e+02 1.862e-03 +7.276e+02 1.858e-03 +7.286e+02 1.853e-03 +7.296e+02 1.849e-03 +7.305e+02 1.845e-03 +7.315e+02 1.840e-03 +7.325e+02 1.836e-03 +7.334e+02 1.832e-03 +7.344e+02 1.827e-03 +7.354e+02 1.823e-03 +7.363e+02 1.819e-03 +7.373e+02 1.814e-03 +7.383e+02 1.810e-03 +7.392e+02 1.806e-03 +7.402e+02 1.802e-03 +7.411e+02 1.797e-03 +7.421e+02 1.793e-03 +7.431e+02 1.789e-03 +7.440e+02 1.785e-03 +7.450e+02 1.781e-03 +7.460e+02 1.776e-03 +7.469e+02 1.772e-03 +7.479e+02 1.768e-03 +7.489e+02 1.764e-03 +7.498e+02 1.760e-03 +7.508e+02 1.756e-03 +7.518e+02 1.752e-03 +7.527e+02 1.748e-03 +7.537e+02 1.744e-03 +7.547e+02 1.740e-03 +7.556e+02 1.736e-03 +7.566e+02 1.732e-03 +7.576e+02 1.728e-03 +7.585e+02 1.724e-03 +7.595e+02 1.720e-03 +7.605e+02 1.716e-03 +7.614e+02 1.712e-03 +7.624e+02 1.708e-03 +7.634e+02 1.704e-03 +7.643e+02 1.700e-03 +7.653e+02 1.697e-03 +7.663e+02 1.693e-03 +7.672e+02 1.689e-03 +7.682e+02 1.685e-03 +7.692e+02 1.681e-03 +7.701e+02 1.677e-03 +7.711e+02 1.674e-03 +7.721e+02 1.670e-03 +7.730e+02 1.666e-03 +7.740e+02 1.662e-03 +7.750e+02 1.659e-03 +7.759e+02 1.655e-03 +7.769e+02 1.651e-03 +7.779e+02 1.648e-03 +7.788e+02 1.644e-03 +7.798e+02 1.640e-03 +7.807e+02 1.637e-03 +7.817e+02 1.633e-03 +7.827e+02 1.629e-03 +7.836e+02 1.626e-03 +7.846e+02 1.622e-03 +7.856e+02 1.619e-03 +7.865e+02 1.615e-03 +7.875e+02 1.611e-03 +7.885e+02 1.608e-03 +7.894e+02 1.604e-03 +7.904e+02 1.601e-03 +7.914e+02 1.597e-03 +7.923e+02 1.594e-03 +7.933e+02 1.590e-03 +7.943e+02 1.587e-03 +7.952e+02 1.583e-03 +7.962e+02 1.580e-03 +7.972e+02 1.576e-03 +7.981e+02 1.573e-03 +7.991e+02 1.570e-03 +8.001e+02 1.566e-03 +8.010e+02 1.563e-03 +8.020e+02 1.559e-03 +8.030e+02 1.556e-03 +8.039e+02 1.553e-03 +8.049e+02 1.549e-03 +8.059e+02 1.546e-03 +8.068e+02 1.543e-03 +8.078e+02 1.539e-03 +8.088e+02 1.536e-03 +8.097e+02 1.533e-03 +8.107e+02 1.529e-03 +8.117e+02 1.526e-03 +8.126e+02 1.523e-03 +8.136e+02 1.520e-03 +8.146e+02 1.516e-03 +8.155e+02 1.513e-03 +8.165e+02 1.510e-03 +8.175e+02 1.507e-03 +8.184e+02 1.504e-03 +8.194e+02 1.500e-03 +8.203e+02 1.497e-03 +8.213e+02 1.494e-03 +8.223e+02 1.491e-03 +8.232e+02 1.488e-03 +8.242e+02 1.485e-03 +8.252e+02 1.481e-03 +8.261e+02 1.478e-03 +8.271e+02 1.475e-03 +8.281e+02 1.472e-03 +8.290e+02 1.469e-03 +8.300e+02 1.466e-03 +8.310e+02 1.463e-03 +8.319e+02 1.460e-03 +8.329e+02 1.457e-03 +8.339e+02 1.454e-03 +8.348e+02 1.451e-03 +8.358e+02 1.448e-03 +8.368e+02 1.445e-03 +8.377e+02 1.442e-03 +8.387e+02 1.439e-03 +8.397e+02 1.436e-03 +8.406e+02 1.433e-03 +8.416e+02 1.430e-03 +8.426e+02 1.427e-03 +8.435e+02 1.424e-03 +8.445e+02 1.421e-03 +8.455e+02 1.418e-03 +8.464e+02 1.415e-03 +8.474e+02 1.412e-03 +8.484e+02 1.409e-03 +8.493e+02 1.406e-03 +8.503e+02 1.404e-03 +8.513e+02 1.401e-03 +8.522e+02 1.398e-03 +8.532e+02 1.395e-03 +8.542e+02 1.392e-03 +8.551e+02 1.389e-03 +8.561e+02 1.387e-03 +8.571e+02 1.384e-03 +8.580e+02 1.381e-03 +8.590e+02 1.378e-03 +8.599e+02 1.375e-03 +8.609e+02 1.373e-03 +8.619e+02 1.370e-03 +8.628e+02 1.367e-03 +8.638e+02 1.364e-03 +8.648e+02 1.362e-03 +8.657e+02 1.359e-03 +8.667e+02 1.356e-03 +8.677e+02 1.353e-03 +8.686e+02 1.351e-03 +8.696e+02 1.348e-03 +8.706e+02 1.345e-03 +8.715e+02 1.343e-03 +8.725e+02 1.340e-03 +8.735e+02 1.337e-03 +8.744e+02 1.335e-03 +8.754e+02 1.332e-03 +8.764e+02 1.329e-03 +8.773e+02 1.327e-03 +8.783e+02 1.324e-03 +8.793e+02 1.321e-03 +8.802e+02 1.319e-03 +8.812e+02 1.316e-03 +8.822e+02 1.314e-03 +8.831e+02 1.311e-03 +8.841e+02 1.308e-03 +8.851e+02 1.306e-03 +8.860e+02 1.303e-03 +8.870e+02 1.301e-03 +8.880e+02 1.298e-03 +8.889e+02 1.296e-03 +8.899e+02 1.293e-03 +8.909e+02 1.291e-03 +8.918e+02 1.288e-03 +8.928e+02 1.286e-03 +8.938e+02 1.283e-03 +8.947e+02 1.281e-03 +8.957e+02 1.278e-03 +8.967e+02 1.276e-03 +8.976e+02 1.273e-03 +8.986e+02 1.271e-03 +8.995e+02 1.268e-03 +9.005e+02 1.266e-03 +9.015e+02 1.263e-03 +9.024e+02 1.261e-03 +9.034e+02 1.259e-03 +9.044e+02 1.256e-03 +9.053e+02 1.254e-03 +9.063e+02 1.251e-03 +9.073e+02 1.249e-03 +9.082e+02 1.247e-03 +9.092e+02 1.244e-03 +9.102e+02 1.242e-03 +9.111e+02 1.239e-03 +9.121e+02 1.237e-03 +9.131e+02 1.235e-03 +9.140e+02 1.232e-03 +9.150e+02 1.230e-03 +9.160e+02 1.228e-03 +9.169e+02 1.225e-03 +9.179e+02 1.223e-03 +9.189e+02 1.221e-03 +9.198e+02 1.218e-03 +9.208e+02 1.216e-03 +9.218e+02 1.214e-03 +9.227e+02 1.212e-03 +9.237e+02 1.209e-03 +9.247e+02 1.207e-03 +9.256e+02 1.205e-03 +9.266e+02 1.202e-03 +9.276e+02 1.200e-03 +9.285e+02 1.198e-03 +9.295e+02 1.196e-03 +9.305e+02 1.193e-03 +9.314e+02 1.191e-03 +9.324e+02 1.189e-03 +9.334e+02 1.187e-03 +9.343e+02 1.185e-03 +9.353e+02 1.182e-03 +9.363e+02 1.180e-03 +9.372e+02 1.178e-03 +9.382e+02 1.176e-03 +9.392e+02 1.174e-03 +9.401e+02 1.171e-03 +9.411e+02 1.169e-03 +9.420e+02 1.167e-03 +9.430e+02 1.165e-03 +9.440e+02 1.163e-03 +9.449e+02 1.161e-03 +9.459e+02 1.159e-03 +9.469e+02 1.156e-03 +9.478e+02 1.154e-03 +9.488e+02 1.152e-03 +9.498e+02 1.150e-03 +9.507e+02 1.148e-03 +9.517e+02 1.146e-03 +9.527e+02 1.144e-03 +9.536e+02 1.142e-03 +9.546e+02 1.140e-03 +9.556e+02 1.138e-03 +9.565e+02 1.136e-03 +9.575e+02 1.133e-03 +9.585e+02 1.131e-03 +9.594e+02 1.129e-03 +9.604e+02 1.127e-03 +9.614e+02 1.125e-03 +9.623e+02 1.123e-03 +9.633e+02 1.121e-03 +9.643e+02 1.119e-03 +9.652e+02 1.117e-03 +9.662e+02 1.115e-03 +9.672e+02 1.113e-03 +9.681e+02 1.111e-03 +9.691e+02 1.109e-03 +9.701e+02 1.107e-03 +9.710e+02 1.105e-03 +9.720e+02 1.103e-03 +9.730e+02 1.101e-03 +9.739e+02 1.099e-03 +9.749e+02 1.097e-03 +9.759e+02 1.095e-03 +9.768e+02 1.093e-03 +9.778e+02 1.092e-03 +9.788e+02 1.090e-03 +9.797e+02 1.088e-03 +9.807e+02 1.086e-03 +9.816e+02 1.084e-03 +9.826e+02 1.082e-03 +9.836e+02 1.080e-03 +9.845e+02 1.078e-03 +9.855e+02 1.076e-03 +9.865e+02 1.074e-03 +9.874e+02 1.072e-03 +9.884e+02 1.070e-03 +9.894e+02 1.069e-03 +9.903e+02 1.067e-03 +9.913e+02 1.065e-03 +9.923e+02 1.063e-03 +9.932e+02 1.061e-03 +9.942e+02 1.059e-03 +9.952e+02 1.057e-03 +9.961e+02 1.056e-03 +9.971e+02 1.054e-03 +9.981e+02 1.052e-03 +9.990e+02 1.050e-03 +1.000e+03 1.048e-03 \ No newline at end of file diff --git a/sed/redlaws/fitzpatrick2004_chiar2006ism_extrapol.dat b/sed/redlaws/fitzpatrick2004_chiar2006ism_extrapol.dat new file mode 100644 index 000000000..915d4bf3c --- /dev/null +++ b/sed/redlaws/fitzpatrick2004_chiar2006ism_extrapol.dat @@ -0,0 +1,1350 @@ +#-- wavelength (micron) A_lambda/A_k +1.149e-01 4.301e+01 +1.161e-01 4.225e+01 +1.173e-01 4.152e+01 +1.185e-01 4.080e+01 +1.197e-01 4.009e+01 +1.210e-01 3.942e+01 +1.223e-01 3.876e+01 +1.236e-01 3.811e+01 +1.249e-01 3.749e+01 +1.263e-01 3.689e+01 +1.277e-01 3.630e+01 +1.291e-01 3.574e+01 +1.306e-01 3.519e+01 +1.321e-01 3.467e+01 +1.337e-01 3.417e+01 +1.352e-01 3.368e+01 +1.368e-01 3.322e+01 +1.385e-01 3.278e+01 +1.402e-01 3.236e+01 +1.419e-01 3.196e+01 +1.437e-01 3.159e+01 +1.455e-01 3.124e+01 +1.474e-01 3.091e+01 +1.493e-01 3.060e+01 +1.512e-01 3.032e+01 +1.533e-01 3.006e+01 +1.553e-01 2.984e+01 +1.575e-01 2.964e+01 +1.596e-01 2.947e+01 +1.619e-01 2.933e+01 +1.642e-01 2.924e+01 +1.666e-01 2.917e+01 +1.690e-01 2.912e+01 +1.716e-01 2.911e+01 +1.742e-01 2.914e+01 +1.768e-01 2.922e+01 +1.796e-01 2.935e+01 +1.824e-01 2.957e+01 +1.854e-01 2.988e+01 +1.884e-01 3.031e+01 +1.916e-01 3.090e+01 +1.948e-01 3.169e+01 +1.982e-01 3.271e+01 +2.017e-01 3.397e+01 +2.053e-01 3.544e+01 +2.090e-01 3.694e+01 +2.129e-01 3.813e+01 +2.169e-01 3.855e+01 +2.210e-01 3.792e+01 +2.254e-01 3.634e+01 +2.299e-01 3.426e+01 +2.346e-01 3.211e+01 +2.395e-01 3.016e+01 +2.446e-01 2.851e+01 +2.499e-01 2.715e+01 +2.554e-01 2.603e+01 +2.612e-01 2.510e+01 +2.673e-01 2.432e+01 +2.737e-01 2.362e+01 +2.803e-01 2.291e+01 +2.874e-01 2.225e+01 +2.947e-01 2.160e+01 +3.025e-01 2.100e+01 +3.107e-01 2.042e+01 +3.193e-01 1.988e+01 +3.284e-01 1.938e+01 +3.381e-01 1.892e+01 +3.483e-01 1.850e+01 +3.592e-01 1.809e+01 +3.708e-01 1.768e+01 +3.831e-01 1.725e+01 +3.964e-01 1.678e+01 +4.105e-01 1.625e+01 +4.257e-01 1.567e+01 +4.421e-01 1.505e+01 +4.598e-01 1.439e+01 +4.789e-01 1.369e+01 +4.998e-01 1.297e+01 +5.225e-01 1.223e+01 +5.473e-01 1.146e+01 +5.747e-01 1.068e+01 +6.050e-01 9.889e+00 +6.386e-01 9.094e+00 +6.761e-01 8.296e+00 +7.184e-01 7.505e+00 +7.663e-01 6.725e+00 +8.210e-01 5.960e+00 +8.842e-01 5.221e+00 +9.579e-01 4.516e+00 +1.045e+00 3.848e+00 +1.149e+00 3.232e+00 +1.240e+00 2.799e+00 +1.340e+00 2.378e+00 +1.440e+00 2.057e+00 +1.540e+00 1.807e+00 +1.640e+00 1.608e+00 +1.740e+00 1.447e+00 +1.840e+00 1.314e+00 +1.940e+00 1.203e+00 +2.040e+00 1.111e+00 +2.140e+00 1.031e+00 +2.240e+00 9.637e-01 +2.340e+00 9.051e-01 +2.440e+00 8.548e-01 +2.540e+00 8.096e-01 +2.640e+00 7.705e-01 +2.740e+00 7.356e-01 +2.840e+00 7.048e-01 +2.940e+00 6.770e-01 +3.040e+00 6.514e-01 +3.140e+00 6.298e-01 +3.240e+00 6.092e-01 +3.340e+00 5.907e-01 +3.440e+00 5.733e-01 +3.540e+00 5.589e-01 +3.640e+00 5.445e-01 +3.740e+00 5.312e-01 +3.840e+00 5.199e-01 +3.940e+00 5.096e-01 +4.040e+00 4.993e-01 +4.140e+00 4.901e-01 +4.240e+00 4.818e-01 +4.340e+00 4.736e-01 +4.440e+00 4.664e-01 +4.540e+00 4.592e-01 +4.640e+00 4.531e-01 +4.740e+00 4.479e-01 +4.840e+00 4.428e-01 +4.940e+00 4.377e-01 +5.040e+00 4.325e-01 +5.140e+00 4.284e-01 +5.240e+00 4.243e-01 +5.340e+00 4.212e-01 +5.440e+00 4.171e-01 +5.540e+00 4.140e-01 +5.640e+00 4.110e-01 +5.740e+00 4.089e-01 +5.840e+00 4.058e-01 +5.940e+00 4.038e-01 +6.040e+00 4.017e-01 +6.140e+00 3.997e-01 +6.240e+00 3.976e-01 +6.340e+00 3.966e-01 +6.440e+00 3.945e-01 +6.540e+00 3.935e-01 +6.640e+00 3.925e-01 +6.740e+00 3.914e-01 +6.840e+00 3.904e-01 +6.940e+00 3.894e-01 +7.040e+00 3.884e-01 +7.140e+00 3.884e-01 +7.240e+00 3.873e-01 +7.340e+00 3.873e-01 +7.440e+00 3.863e-01 +7.540e+00 3.863e-01 +7.640e+00 3.863e-01 +7.740e+00 3.863e-01 +7.840e+00 3.863e-01 +7.940e+00 3.863e-01 +8.100e+00 4.181e-01 +8.200e+00 4.397e-01 +8.300e+00 4.623e-01 +8.400e+00 4.880e-01 +8.500e+00 5.199e-01 +8.600e+00 5.568e-01 +8.700e+00 5.979e-01 +8.800e+00 6.390e-01 +8.900e+00 6.812e-01 +9.000e+00 7.264e-01 +9.100e+00 7.777e-01 +9.200e+00 8.363e-01 +9.300e+00 8.959e-01 +9.400e+00 9.421e-01 +9.500e+00 9.719e-01 +9.600e+00 9.986e-01 +9.700e+00 1.023e+00 +9.800e+00 1.034e+00 +9.900e+00 1.027e+00 +1.000e+01 1.010e+00 +1.010e+01 9.873e-01 +1.020e+01 9.637e-01 +1.030e+01 9.390e-01 +1.040e+01 9.133e-01 +1.050e+01 8.887e-01 +1.060e+01 8.661e-01 +1.070e+01 8.455e-01 +1.080e+01 8.270e-01 +1.090e+01 8.096e-01 +1.100e+01 7.921e-01 +1.110e+01 7.746e-01 +1.120e+01 7.572e-01 +1.130e+01 7.377e-01 +1.140e+01 7.181e-01 +1.150e+01 6.986e-01 +1.160e+01 6.812e-01 +1.170e+01 6.647e-01 +1.180e+01 6.503e-01 +1.190e+01 6.370e-01 +1.200e+01 6.257e-01 +1.210e+01 6.144e-01 +1.220e+01 6.020e-01 +1.230e+01 5.907e-01 +1.240e+01 5.794e-01 +1.250e+01 5.702e-01 +1.260e+01 5.620e-01 +1.270e+01 5.548e-01 +1.280e+01 5.507e-01 +1.290e+01 5.476e-01 +1.300e+01 5.445e-01 +1.310e+01 5.435e-01 +1.320e+01 5.425e-01 +1.330e+01 5.414e-01 +1.340e+01 5.425e-01 +1.350e+01 5.425e-01 +1.360e+01 5.435e-01 +1.370e+01 5.445e-01 +1.380e+01 5.466e-01 +1.390e+01 5.486e-01 +1.400e+01 5.507e-01 +1.410e+01 5.538e-01 +1.420e+01 5.558e-01 +1.430e+01 5.589e-01 +1.440e+01 5.610e-01 +1.450e+01 5.640e-01 +1.460e+01 5.671e-01 +1.470e+01 5.702e-01 +1.480e+01 5.733e-01 +1.490e+01 5.764e-01 +1.500e+01 5.794e-01 +1.510e+01 5.836e-01 +1.520e+01 5.866e-01 +1.530e+01 5.907e-01 +1.540e+01 5.938e-01 +1.550e+01 5.979e-01 +1.560e+01 6.010e-01 +1.570e+01 6.051e-01 +1.580e+01 6.092e-01 +1.590e+01 6.123e-01 +1.600e+01 6.164e-01 +1.610e+01 6.205e-01 +1.620e+01 6.246e-01 +1.630e+01 6.288e-01 +1.640e+01 6.329e-01 +1.650e+01 6.370e-01 +1.660e+01 6.421e-01 +1.670e+01 6.462e-01 +1.680e+01 6.503e-01 +1.690e+01 6.544e-01 +1.700e+01 6.596e-01 +1.710e+01 6.637e-01 +1.720e+01 6.678e-01 +1.730e+01 6.729e-01 +1.740e+01 6.770e-01 +1.750e+01 6.812e-01 +1.760e+01 6.863e-01 +1.770e+01 6.904e-01 +1.780e+01 6.945e-01 +1.790e+01 6.986e-01 +1.800e+01 7.017e-01 +1.810e+01 7.058e-01 +1.820e+01 7.089e-01 +1.830e+01 7.120e-01 +1.840e+01 7.140e-01 +1.850e+01 7.171e-01 +1.860e+01 7.181e-01 +1.870e+01 7.202e-01 +1.880e+01 7.212e-01 +1.890e+01 7.222e-01 +1.900e+01 7.222e-01 +1.910e+01 7.222e-01 +1.920e+01 7.212e-01 +1.930e+01 7.202e-01 +1.940e+01 7.192e-01 +1.950e+01 7.181e-01 +1.960e+01 7.161e-01 +1.970e+01 7.140e-01 +1.980e+01 7.120e-01 +1.990e+01 7.089e-01 +2.000e+01 7.068e-01 +2.010e+01 7.038e-01 +2.020e+01 6.996e-01 +2.030e+01 6.966e-01 +2.040e+01 6.925e-01 +2.050e+01 6.894e-01 +2.060e+01 6.853e-01 +2.070e+01 6.812e-01 +2.080e+01 6.770e-01 +2.090e+01 6.729e-01 +2.100e+01 6.678e-01 +2.110e+01 6.637e-01 +2.120e+01 6.596e-01 +2.130e+01 6.555e-01 +2.140e+01 6.503e-01 +2.150e+01 6.462e-01 +2.160e+01 6.421e-01 +2.170e+01 6.380e-01 +2.180e+01 6.339e-01 +2.190e+01 6.298e-01 +2.200e+01 6.257e-01 +2.210e+01 6.216e-01 +2.220e+01 6.185e-01 +2.230e+01 6.144e-01 +2.240e+01 6.103e-01 +2.250e+01 6.062e-01 +2.260e+01 6.031e-01 +2.270e+01 5.990e-01 +2.280e+01 5.959e-01 +2.290e+01 5.918e-01 +2.300e+01 5.887e-01 +2.310e+01 5.846e-01 +2.320e+01 5.815e-01 +2.330e+01 5.774e-01 +2.340e+01 5.743e-01 +2.350e+01 5.702e-01 +2.360e+01 5.671e-01 +2.370e+01 5.630e-01 +2.380e+01 5.599e-01 +2.390e+01 5.568e-01 +2.400e+01 5.527e-01 +2.410e+01 5.496e-01 +2.420e+01 5.466e-01 +2.430e+01 5.435e-01 +2.440e+01 5.394e-01 +2.450e+01 5.363e-01 +2.460e+01 5.332e-01 +2.470e+01 5.301e-01 +2.480e+01 5.270e-01 +2.490e+01 5.240e-01 +2.500e+01 5.209e-01 +2.510e+01 5.178e-01 +2.520e+01 5.147e-01 +2.530e+01 5.116e-01 +2.540e+01 5.086e-01 +2.550e+01 5.055e-01 +2.560e+01 5.024e-01 +2.570e+01 5.003e-01 +2.580e+01 4.973e-01 +2.590e+01 4.942e-01 +2.600e+01 4.911e-01 +2.610e+01 4.890e-01 +2.620e+01 4.860e-01 +2.630e+01 4.829e-01 +2.640e+01 4.808e-01 +2.650e+01 4.777e-01 +2.660e+01 4.747e-01 +2.670e+01 4.726e-01 +2.680e+01 4.695e-01 +2.690e+01 4.675e-01 +2.700e+01 4.644e-01 +2.710e+01 4.613e-01 +2.807e+01 4.329e-01 +2.905e+01 4.071e-01 +3.002e+01 3.837e-01 +3.100e+01 3.622e-01 +3.197e+01 3.426e-01 +3.294e+01 3.246e-01 +3.392e+01 3.080e-01 +3.489e+01 2.927e-01 +3.586e+01 2.786e-01 +3.684e+01 2.654e-01 +3.781e+01 2.533e-01 +3.879e+01 2.419e-01 +3.976e+01 2.314e-01 +4.073e+01 2.215e-01 +4.171e+01 2.123e-01 +4.268e+01 2.037e-01 +4.366e+01 1.955e-01 +4.463e+01 1.879e-01 +4.560e+01 1.808e-01 +4.658e+01 1.740e-01 +4.755e+01 1.677e-01 +4.853e+01 1.617e-01 +4.950e+01 1.560e-01 +5.047e+01 1.506e-01 +5.145e+01 1.455e-01 +5.242e+01 1.407e-01 +5.339e+01 1.361e-01 +5.437e+01 1.317e-01 +5.534e+01 1.276e-01 +5.632e+01 1.236e-01 +5.729e+01 1.199e-01 +5.826e+01 1.163e-01 +5.924e+01 1.129e-01 +6.021e+01 1.096e-01 +6.119e+01 1.065e-01 +6.216e+01 1.035e-01 +6.313e+01 1.007e-01 +6.411e+01 9.792e-02 +6.508e+01 9.530e-02 +6.605e+01 9.279e-02 +6.703e+01 9.038e-02 +6.800e+01 8.806e-02 +6.898e+01 8.583e-02 +6.995e+01 8.370e-02 +7.092e+01 8.164e-02 +7.190e+01 7.966e-02 +7.287e+01 7.775e-02 +7.385e+01 7.592e-02 +7.482e+01 7.415e-02 +7.579e+01 7.244e-02 +7.677e+01 7.080e-02 +7.774e+01 6.921e-02 +7.872e+01 6.767e-02 +7.969e+01 6.619e-02 +8.066e+01 6.476e-02 +8.164e+01 6.338e-02 +8.261e+01 6.204e-02 +8.358e+01 6.074e-02 +8.456e+01 5.949e-02 +8.553e+01 5.828e-02 +8.651e+01 5.710e-02 +8.748e+01 5.596e-02 +8.845e+01 5.486e-02 +8.943e+01 5.379e-02 +9.040e+01 5.275e-02 +9.138e+01 5.174e-02 +9.235e+01 5.076e-02 +9.332e+01 4.981e-02 +9.430e+01 4.889e-02 +9.527e+01 4.799e-02 +9.625e+01 4.712e-02 +9.722e+01 4.628e-02 +9.819e+01 4.546e-02 +9.917e+01 4.465e-02 +1.001e+02 4.388e-02 +1.011e+02 4.312e-02 +1.021e+02 4.238e-02 +1.031e+02 4.166e-02 +1.040e+02 4.096e-02 +1.050e+02 4.028e-02 +1.060e+02 3.962e-02 +1.070e+02 3.897e-02 +1.079e+02 3.834e-02 +1.089e+02 3.773e-02 +1.099e+02 3.713e-02 +1.109e+02 3.654e-02 +1.118e+02 3.597e-02 +1.128e+02 3.541e-02 +1.138e+02 3.487e-02 +1.147e+02 3.434e-02 +1.157e+02 3.382e-02 +1.167e+02 3.331e-02 +1.177e+02 3.282e-02 +1.186e+02 3.234e-02 +1.196e+02 3.186e-02 +1.206e+02 3.140e-02 +1.216e+02 3.095e-02 +1.225e+02 3.051e-02 +1.235e+02 3.008e-02 +1.245e+02 2.966e-02 +1.255e+02 2.924e-02 +1.264e+02 2.884e-02 +1.274e+02 2.844e-02 +1.284e+02 2.806e-02 +1.294e+02 2.768e-02 +1.303e+02 2.730e-02 +1.313e+02 2.694e-02 +1.323e+02 2.659e-02 +1.333e+02 2.624e-02 +1.342e+02 2.590e-02 +1.352e+02 2.556e-02 +1.362e+02 2.523e-02 +1.371e+02 2.491e-02 +1.381e+02 2.460e-02 +1.391e+02 2.429e-02 +1.401e+02 2.398e-02 +1.410e+02 2.369e-02 +1.420e+02 2.339e-02 +1.430e+02 2.311e-02 +1.440e+02 2.283e-02 +1.449e+02 2.255e-02 +1.459e+02 2.228e-02 +1.469e+02 2.202e-02 +1.479e+02 2.176e-02 +1.488e+02 2.150e-02 +1.498e+02 2.125e-02 +1.508e+02 2.100e-02 +1.518e+02 2.076e-02 +1.527e+02 2.052e-02 +1.537e+02 2.029e-02 +1.547e+02 2.006e-02 +1.557e+02 1.984e-02 +1.566e+02 1.961e-02 +1.576e+02 1.940e-02 +1.586e+02 1.918e-02 +1.595e+02 1.897e-02 +1.605e+02 1.877e-02 +1.615e+02 1.856e-02 +1.625e+02 1.836e-02 +1.634e+02 1.817e-02 +1.644e+02 1.797e-02 +1.654e+02 1.778e-02 +1.664e+02 1.760e-02 +1.673e+02 1.741e-02 +1.683e+02 1.723e-02 +1.693e+02 1.705e-02 +1.703e+02 1.688e-02 +1.712e+02 1.671e-02 +1.722e+02 1.654e-02 +1.732e+02 1.637e-02 +1.742e+02 1.620e-02 +1.751e+02 1.604e-02 +1.761e+02 1.588e-02 +1.771e+02 1.573e-02 +1.781e+02 1.557e-02 +1.790e+02 1.542e-02 +1.800e+02 1.527e-02 +1.810e+02 1.512e-02 +1.819e+02 1.498e-02 +1.829e+02 1.483e-02 +1.839e+02 1.469e-02 +1.849e+02 1.455e-02 +1.858e+02 1.442e-02 +1.868e+02 1.428e-02 +1.878e+02 1.415e-02 +1.888e+02 1.402e-02 +1.897e+02 1.389e-02 +1.907e+02 1.376e-02 +1.917e+02 1.364e-02 +1.927e+02 1.351e-02 +1.936e+02 1.339e-02 +1.946e+02 1.327e-02 +1.956e+02 1.315e-02 +1.966e+02 1.303e-02 +1.975e+02 1.292e-02 +1.985e+02 1.280e-02 +1.995e+02 1.269e-02 +2.004e+02 1.258e-02 +2.014e+02 1.247e-02 +2.024e+02 1.236e-02 +2.034e+02 1.226e-02 +2.043e+02 1.215e-02 +2.053e+02 1.205e-02 +2.063e+02 1.195e-02 +2.073e+02 1.185e-02 +2.082e+02 1.175e-02 +2.092e+02 1.165e-02 +2.102e+02 1.155e-02 +2.112e+02 1.146e-02 +2.121e+02 1.136e-02 +2.131e+02 1.127e-02 +2.141e+02 1.118e-02 +2.151e+02 1.108e-02 +2.160e+02 1.099e-02 +2.170e+02 1.091e-02 +2.180e+02 1.082e-02 +2.190e+02 1.073e-02 +2.199e+02 1.065e-02 +2.209e+02 1.056e-02 +2.219e+02 1.048e-02 +2.228e+02 1.040e-02 +2.238e+02 1.032e-02 +2.248e+02 1.024e-02 +2.258e+02 1.016e-02 +2.267e+02 1.008e-02 +2.277e+02 1.000e-02 +2.287e+02 9.924e-03 +2.297e+02 9.848e-03 +2.306e+02 9.773e-03 +2.316e+02 9.700e-03 +2.326e+02 9.627e-03 +2.336e+02 9.554e-03 +2.345e+02 9.483e-03 +2.355e+02 9.413e-03 +2.365e+02 9.343e-03 +2.375e+02 9.274e-03 +2.384e+02 9.206e-03 +2.394e+02 9.139e-03 +2.404e+02 9.072e-03 +2.414e+02 9.006e-03 +2.423e+02 8.941e-03 +2.433e+02 8.877e-03 +2.443e+02 8.813e-03 +2.452e+02 8.751e-03 +2.462e+02 8.688e-03 +2.472e+02 8.627e-03 +2.482e+02 8.566e-03 +2.491e+02 8.506e-03 +2.501e+02 8.446e-03 +2.511e+02 8.387e-03 +2.521e+02 8.329e-03 +2.530e+02 8.272e-03 +2.540e+02 8.215e-03 +2.550e+02 8.158e-03 +2.560e+02 8.102e-03 +2.569e+02 8.047e-03 +2.579e+02 7.993e-03 +2.589e+02 7.939e-03 +2.599e+02 7.885e-03 +2.608e+02 7.832e-03 +2.618e+02 7.780e-03 +2.628e+02 7.728e-03 +2.638e+02 7.677e-03 +2.647e+02 7.626e-03 +2.657e+02 7.576e-03 +2.667e+02 7.526e-03 +2.676e+02 7.477e-03 +2.686e+02 7.428e-03 +2.696e+02 7.380e-03 +2.706e+02 7.332e-03 +2.715e+02 7.285e-03 +2.725e+02 7.238e-03 +2.735e+02 7.192e-03 +2.745e+02 7.146e-03 +2.754e+02 7.100e-03 +2.764e+02 7.055e-03 +2.774e+02 7.011e-03 +2.784e+02 6.967e-03 +2.793e+02 6.923e-03 +2.803e+02 6.880e-03 +2.813e+02 6.837e-03 +2.823e+02 6.795e-03 +2.832e+02 6.753e-03 +2.842e+02 6.711e-03 +2.852e+02 6.670e-03 +2.862e+02 6.629e-03 +2.871e+02 6.589e-03 +2.881e+02 6.549e-03 +2.891e+02 6.509e-03 +2.900e+02 6.470e-03 +2.910e+02 6.431e-03 +2.920e+02 6.392e-03 +2.930e+02 6.354e-03 +2.939e+02 6.316e-03 +2.949e+02 6.279e-03 +2.959e+02 6.242e-03 +2.969e+02 6.205e-03 +2.978e+02 6.168e-03 +2.988e+02 6.132e-03 +2.998e+02 6.096e-03 +3.008e+02 6.061e-03 +3.017e+02 6.026e-03 +3.027e+02 5.991e-03 +3.037e+02 5.956e-03 +3.047e+02 5.922e-03 +3.056e+02 5.888e-03 +3.066e+02 5.855e-03 +3.076e+02 5.821e-03 +3.085e+02 5.788e-03 +3.095e+02 5.755e-03 +3.105e+02 5.723e-03 +3.115e+02 5.691e-03 +3.124e+02 5.659e-03 +3.134e+02 5.627e-03 +3.144e+02 5.596e-03 +3.154e+02 5.565e-03 +3.163e+02 5.534e-03 +3.173e+02 5.504e-03 +3.183e+02 5.473e-03 +3.193e+02 5.443e-03 +3.202e+02 5.414e-03 +3.212e+02 5.384e-03 +3.222e+02 5.355e-03 +3.232e+02 5.326e-03 +3.241e+02 5.297e-03 +3.251e+02 5.268e-03 +3.261e+02 5.240e-03 +3.271e+02 5.212e-03 +3.280e+02 5.184e-03 +3.290e+02 5.157e-03 +3.300e+02 5.129e-03 +3.309e+02 5.102e-03 +3.319e+02 5.075e-03 +3.329e+02 5.049e-03 +3.339e+02 5.022e-03 +3.348e+02 4.996e-03 +3.358e+02 4.970e-03 +3.368e+02 4.944e-03 +3.378e+02 4.918e-03 +3.387e+02 4.893e-03 +3.397e+02 4.868e-03 +3.407e+02 4.843e-03 +3.417e+02 4.818e-03 +3.426e+02 4.793e-03 +3.436e+02 4.769e-03 +3.446e+02 4.745e-03 +3.456e+02 4.721e-03 +3.465e+02 4.697e-03 +3.475e+02 4.673e-03 +3.485e+02 4.650e-03 +3.495e+02 4.626e-03 +3.504e+02 4.603e-03 +3.514e+02 4.580e-03 +3.524e+02 4.557e-03 +3.533e+02 4.535e-03 +3.543e+02 4.512e-03 +3.553e+02 4.490e-03 +3.563e+02 4.468e-03 +3.572e+02 4.446e-03 +3.582e+02 4.424e-03 +3.592e+02 4.403e-03 +3.602e+02 4.381e-03 +3.611e+02 4.360e-03 +3.621e+02 4.339e-03 +3.631e+02 4.318e-03 +3.641e+02 4.297e-03 +3.650e+02 4.277e-03 +3.660e+02 4.256e-03 +3.670e+02 4.236e-03 +3.680e+02 4.216e-03 +3.689e+02 4.196e-03 +3.699e+02 4.176e-03 +3.709e+02 4.156e-03 +3.719e+02 4.137e-03 +3.728e+02 4.117e-03 +3.738e+02 4.098e-03 +3.748e+02 4.079e-03 +3.757e+02 4.060e-03 +3.767e+02 4.041e-03 +3.777e+02 4.022e-03 +3.787e+02 4.004e-03 +3.796e+02 3.985e-03 +3.806e+02 3.967e-03 +3.816e+02 3.949e-03 +3.826e+02 3.931e-03 +3.835e+02 3.913e-03 +3.845e+02 3.895e-03 +3.855e+02 3.877e-03 +3.865e+02 3.860e-03 +3.874e+02 3.842e-03 +3.884e+02 3.825e-03 +3.894e+02 3.808e-03 +3.904e+02 3.791e-03 +3.913e+02 3.774e-03 +3.923e+02 3.757e-03 +3.933e+02 3.740e-03 +3.943e+02 3.723e-03 +3.952e+02 3.707e-03 +3.962e+02 3.690e-03 +3.972e+02 3.674e-03 +3.981e+02 3.658e-03 +3.991e+02 3.642e-03 +4.001e+02 3.626e-03 +4.011e+02 3.610e-03 +4.020e+02 3.594e-03 +4.030e+02 3.579e-03 +4.040e+02 3.563e-03 +4.050e+02 3.548e-03 +4.059e+02 3.533e-03 +4.069e+02 3.517e-03 +4.079e+02 3.502e-03 +4.089e+02 3.487e-03 +4.098e+02 3.472e-03 +4.108e+02 3.458e-03 +4.118e+02 3.443e-03 +4.128e+02 3.428e-03 +4.137e+02 3.414e-03 +4.147e+02 3.399e-03 +4.157e+02 3.385e-03 +4.166e+02 3.371e-03 +4.176e+02 3.357e-03 +4.186e+02 3.343e-03 +4.196e+02 3.329e-03 +4.205e+02 3.315e-03 +4.215e+02 3.301e-03 +4.225e+02 3.287e-03 +4.235e+02 3.274e-03 +4.244e+02 3.260e-03 +4.254e+02 3.247e-03 +4.264e+02 3.234e-03 +4.274e+02 3.220e-03 +4.283e+02 3.207e-03 +4.293e+02 3.194e-03 +4.303e+02 3.181e-03 +4.313e+02 3.168e-03 +4.322e+02 3.155e-03 +4.332e+02 3.143e-03 +4.342e+02 3.130e-03 +4.352e+02 3.117e-03 +4.361e+02 3.105e-03 +4.371e+02 3.092e-03 +4.381e+02 3.080e-03 +4.390e+02 3.068e-03 +4.400e+02 3.055e-03 +4.410e+02 3.043e-03 +4.420e+02 3.031e-03 +4.429e+02 3.019e-03 +4.439e+02 3.007e-03 +4.449e+02 2.995e-03 +4.459e+02 2.984e-03 +4.468e+02 2.972e-03 +4.478e+02 2.960e-03 +4.488e+02 2.949e-03 +4.498e+02 2.937e-03 +4.507e+02 2.926e-03 +4.517e+02 2.915e-03 +4.527e+02 2.903e-03 +4.537e+02 2.892e-03 +4.546e+02 2.881e-03 +4.556e+02 2.870e-03 +4.566e+02 2.859e-03 +4.576e+02 2.848e-03 +4.585e+02 2.837e-03 +4.595e+02 2.826e-03 +4.605e+02 2.816e-03 +4.614e+02 2.805e-03 +4.624e+02 2.794e-03 +4.634e+02 2.784e-03 +4.644e+02 2.773e-03 +4.653e+02 2.763e-03 +4.663e+02 2.752e-03 +4.673e+02 2.742e-03 +4.683e+02 2.732e-03 +4.692e+02 2.722e-03 +4.702e+02 2.711e-03 +4.712e+02 2.701e-03 +4.722e+02 2.691e-03 +4.731e+02 2.681e-03 +4.741e+02 2.671e-03 +4.751e+02 2.662e-03 +4.761e+02 2.652e-03 +4.770e+02 2.642e-03 +4.780e+02 2.632e-03 +4.790e+02 2.623e-03 +4.800e+02 2.613e-03 +4.809e+02 2.604e-03 +4.819e+02 2.594e-03 +4.829e+02 2.585e-03 +4.838e+02 2.575e-03 +4.848e+02 2.566e-03 +4.858e+02 2.557e-03 +4.868e+02 2.548e-03 +4.877e+02 2.539e-03 +4.887e+02 2.529e-03 +4.897e+02 2.520e-03 +4.907e+02 2.511e-03 +4.916e+02 2.502e-03 +4.926e+02 2.494e-03 +4.936e+02 2.485e-03 +4.946e+02 2.476e-03 +4.955e+02 2.467e-03 +4.965e+02 2.458e-03 +4.975e+02 2.450e-03 +4.985e+02 2.441e-03 +4.994e+02 2.433e-03 +5.004e+02 2.424e-03 +5.014e+02 2.416e-03 +5.024e+02 2.407e-03 +5.033e+02 2.399e-03 +5.043e+02 2.390e-03 +5.053e+02 2.382e-03 +5.062e+02 2.374e-03 +5.072e+02 2.366e-03 +5.082e+02 2.358e-03 +5.092e+02 2.350e-03 +5.101e+02 2.341e-03 +5.111e+02 2.333e-03 +5.121e+02 2.325e-03 +5.131e+02 2.317e-03 +5.140e+02 2.310e-03 +5.150e+02 2.302e-03 +5.160e+02 2.294e-03 +5.170e+02 2.286e-03 +5.179e+02 2.278e-03 +5.189e+02 2.271e-03 +5.199e+02 2.263e-03 +5.209e+02 2.255e-03 +5.218e+02 2.248e-03 +5.228e+02 2.240e-03 +5.238e+02 2.233e-03 +5.247e+02 2.225e-03 +5.257e+02 2.218e-03 +5.267e+02 2.211e-03 +5.277e+02 2.203e-03 +5.286e+02 2.196e-03 +5.296e+02 2.189e-03 +5.306e+02 2.181e-03 +5.316e+02 2.174e-03 +5.325e+02 2.167e-03 +5.335e+02 2.160e-03 +5.345e+02 2.153e-03 +5.355e+02 2.146e-03 +5.364e+02 2.139e-03 +5.374e+02 2.132e-03 +5.384e+02 2.125e-03 +5.394e+02 2.118e-03 +5.403e+02 2.111e-03 +5.413e+02 2.104e-03 +5.423e+02 2.098e-03 +5.433e+02 2.091e-03 +5.442e+02 2.084e-03 +5.452e+02 2.077e-03 +5.462e+02 2.071e-03 +5.471e+02 2.064e-03 +5.481e+02 2.058e-03 +5.491e+02 2.051e-03 +5.501e+02 2.044e-03 +5.510e+02 2.038e-03 +5.520e+02 2.031e-03 +5.530e+02 2.025e-03 +5.540e+02 2.019e-03 +5.549e+02 2.012e-03 +5.559e+02 2.006e-03 +5.569e+02 2.000e-03 +5.579e+02 1.993e-03 +5.588e+02 1.987e-03 +5.598e+02 1.981e-03 +5.608e+02 1.975e-03 +5.618e+02 1.969e-03 +5.627e+02 1.962e-03 +5.637e+02 1.956e-03 +5.647e+02 1.950e-03 +5.657e+02 1.944e-03 +5.666e+02 1.938e-03 +5.676e+02 1.932e-03 +5.686e+02 1.926e-03 +5.695e+02 1.920e-03 +5.705e+02 1.914e-03 +5.715e+02 1.909e-03 +5.725e+02 1.903e-03 +5.734e+02 1.897e-03 +5.744e+02 1.891e-03 +5.754e+02 1.885e-03 +5.764e+02 1.880e-03 +5.773e+02 1.874e-03 +5.783e+02 1.868e-03 +5.793e+02 1.863e-03 +5.803e+02 1.857e-03 +5.812e+02 1.851e-03 +5.822e+02 1.846e-03 +5.832e+02 1.840e-03 +5.842e+02 1.835e-03 +5.851e+02 1.829e-03 +5.861e+02 1.824e-03 +5.871e+02 1.818e-03 +5.881e+02 1.813e-03 +5.890e+02 1.808e-03 +5.900e+02 1.802e-03 +5.910e+02 1.797e-03 +5.919e+02 1.791e-03 +5.929e+02 1.786e-03 +5.939e+02 1.781e-03 +5.949e+02 1.776e-03 +5.958e+02 1.770e-03 +5.968e+02 1.765e-03 +5.978e+02 1.760e-03 +5.988e+02 1.755e-03 +5.997e+02 1.750e-03 +6.007e+02 1.745e-03 +6.017e+02 1.740e-03 +6.027e+02 1.735e-03 +6.036e+02 1.730e-03 +6.046e+02 1.725e-03 +6.056e+02 1.720e-03 +6.066e+02 1.715e-03 +6.075e+02 1.710e-03 +6.085e+02 1.705e-03 +6.095e+02 1.700e-03 +6.105e+02 1.695e-03 +6.114e+02 1.690e-03 +6.124e+02 1.685e-03 +6.134e+02 1.680e-03 +6.143e+02 1.676e-03 +6.153e+02 1.671e-03 +6.163e+02 1.666e-03 +6.173e+02 1.661e-03 +6.182e+02 1.657e-03 +6.192e+02 1.652e-03 +6.202e+02 1.647e-03 +6.212e+02 1.643e-03 +6.221e+02 1.638e-03 +6.231e+02 1.633e-03 +6.241e+02 1.629e-03 +6.251e+02 1.624e-03 +6.260e+02 1.620e-03 +6.270e+02 1.615e-03 +6.280e+02 1.611e-03 +6.290e+02 1.606e-03 +6.299e+02 1.602e-03 +6.309e+02 1.597e-03 +6.319e+02 1.593e-03 +6.328e+02 1.588e-03 +6.338e+02 1.584e-03 +6.348e+02 1.580e-03 +6.358e+02 1.575e-03 +6.367e+02 1.571e-03 +6.377e+02 1.567e-03 +6.387e+02 1.562e-03 +6.397e+02 1.558e-03 +6.406e+02 1.554e-03 +6.416e+02 1.550e-03 +6.426e+02 1.545e-03 +6.436e+02 1.541e-03 +6.445e+02 1.537e-03 +6.455e+02 1.533e-03 +6.465e+02 1.529e-03 +6.475e+02 1.525e-03 +6.484e+02 1.520e-03 +6.494e+02 1.516e-03 +6.504e+02 1.512e-03 +6.514e+02 1.508e-03 +6.523e+02 1.504e-03 +6.533e+02 1.500e-03 +6.543e+02 1.496e-03 +6.552e+02 1.492e-03 +6.562e+02 1.488e-03 +6.572e+02 1.484e-03 +6.582e+02 1.480e-03 +6.591e+02 1.476e-03 +6.601e+02 1.472e-03 +6.611e+02 1.468e-03 +6.621e+02 1.465e-03 +6.630e+02 1.461e-03 +6.640e+02 1.457e-03 +6.650e+02 1.453e-03 +6.660e+02 1.449e-03 +6.669e+02 1.445e-03 +6.679e+02 1.442e-03 +6.689e+02 1.438e-03 +6.699e+02 1.434e-03 +6.708e+02 1.430e-03 +6.718e+02 1.427e-03 +6.728e+02 1.423e-03 +6.738e+02 1.419e-03 +6.747e+02 1.415e-03 +6.757e+02 1.412e-03 +6.767e+02 1.408e-03 +6.776e+02 1.404e-03 +6.786e+02 1.401e-03 +6.796e+02 1.397e-03 +6.806e+02 1.394e-03 +6.815e+02 1.390e-03 +6.825e+02 1.387e-03 +6.835e+02 1.383e-03 +6.845e+02 1.379e-03 +6.854e+02 1.376e-03 +6.864e+02 1.372e-03 +6.874e+02 1.369e-03 +6.884e+02 1.365e-03 +6.893e+02 1.362e-03 +6.903e+02 1.358e-03 +6.913e+02 1.355e-03 +6.923e+02 1.352e-03 +6.932e+02 1.348e-03 +6.942e+02 1.345e-03 +6.952e+02 1.341e-03 +6.962e+02 1.338e-03 +6.971e+02 1.335e-03 +6.981e+02 1.331e-03 +6.991e+02 1.328e-03 +7.000e+02 1.325e-03 +7.010e+02 1.321e-03 +7.020e+02 1.318e-03 +7.030e+02 1.315e-03 +7.039e+02 1.311e-03 +7.049e+02 1.308e-03 +7.059e+02 1.305e-03 +7.069e+02 1.302e-03 +7.078e+02 1.299e-03 +7.088e+02 1.295e-03 +7.098e+02 1.292e-03 +7.108e+02 1.289e-03 +7.117e+02 1.286e-03 +7.127e+02 1.283e-03 +7.137e+02 1.279e-03 +7.147e+02 1.276e-03 +7.156e+02 1.273e-03 +7.166e+02 1.270e-03 +7.176e+02 1.267e-03 +7.186e+02 1.264e-03 +7.195e+02 1.261e-03 +7.205e+02 1.258e-03 +7.215e+02 1.255e-03 +7.224e+02 1.252e-03 +7.234e+02 1.249e-03 +7.244e+02 1.246e-03 +7.254e+02 1.243e-03 +7.263e+02 1.240e-03 +7.273e+02 1.237e-03 +7.283e+02 1.234e-03 +7.293e+02 1.231e-03 +7.302e+02 1.228e-03 +7.312e+02 1.225e-03 +7.322e+02 1.222e-03 +7.332e+02 1.219e-03 +7.341e+02 1.216e-03 +7.351e+02 1.213e-03 +7.361e+02 1.210e-03 +7.371e+02 1.207e-03 +7.380e+02 1.204e-03 +7.390e+02 1.202e-03 +7.400e+02 1.199e-03 +7.409e+02 1.196e-03 +7.419e+02 1.193e-03 +7.429e+02 1.190e-03 +7.439e+02 1.187e-03 +7.448e+02 1.185e-03 +7.458e+02 1.182e-03 +7.468e+02 1.179e-03 +7.478e+02 1.176e-03 +7.487e+02 1.174e-03 +7.497e+02 1.171e-03 +7.507e+02 1.168e-03 +7.517e+02 1.165e-03 +7.526e+02 1.163e-03 +7.536e+02 1.160e-03 +7.546e+02 1.157e-03 +7.556e+02 1.155e-03 +7.565e+02 1.152e-03 +7.575e+02 1.149e-03 +7.585e+02 1.147e-03 +7.595e+02 1.144e-03 +7.604e+02 1.141e-03 +7.614e+02 1.139e-03 +7.624e+02 1.136e-03 +7.633e+02 1.134e-03 +7.643e+02 1.131e-03 +7.653e+02 1.128e-03 +7.663e+02 1.126e-03 +7.672e+02 1.123e-03 +7.682e+02 1.121e-03 +7.692e+02 1.118e-03 +7.702e+02 1.116e-03 +7.711e+02 1.113e-03 +7.721e+02 1.110e-03 +7.731e+02 1.108e-03 +7.741e+02 1.105e-03 +7.750e+02 1.103e-03 +7.760e+02 1.100e-03 +7.770e+02 1.098e-03 +7.780e+02 1.095e-03 +7.789e+02 1.093e-03 +7.799e+02 1.091e-03 +7.809e+02 1.088e-03 +7.819e+02 1.086e-03 +7.828e+02 1.083e-03 +7.838e+02 1.081e-03 +7.848e+02 1.078e-03 +7.857e+02 1.076e-03 +7.867e+02 1.074e-03 +7.877e+02 1.071e-03 +7.887e+02 1.069e-03 +7.896e+02 1.066e-03 +7.906e+02 1.064e-03 +7.916e+02 1.062e-03 +7.926e+02 1.059e-03 +7.935e+02 1.057e-03 +7.945e+02 1.055e-03 +7.955e+02 1.052e-03 +7.965e+02 1.050e-03 +7.974e+02 1.048e-03 +7.984e+02 1.045e-03 +7.994e+02 1.043e-03 +8.004e+02 1.041e-03 +8.013e+02 1.039e-03 +8.023e+02 1.036e-03 +8.033e+02 1.034e-03 +8.043e+02 1.032e-03 +8.052e+02 1.030e-03 +8.062e+02 1.027e-03 +8.072e+02 1.025e-03 +8.081e+02 1.023e-03 +8.091e+02 1.021e-03 +8.101e+02 1.018e-03 +8.111e+02 1.016e-03 +8.120e+02 1.014e-03 +8.130e+02 1.012e-03 +8.140e+02 1.010e-03 +8.150e+02 1.008e-03 +8.159e+02 1.005e-03 +8.169e+02 1.003e-03 +8.179e+02 1.001e-03 +8.189e+02 9.990e-04 +8.198e+02 9.968e-04 +8.208e+02 9.947e-04 +8.218e+02 9.926e-04 +8.228e+02 9.905e-04 +8.237e+02 9.884e-04 +8.247e+02 9.863e-04 +8.257e+02 9.842e-04 +8.267e+02 9.821e-04 +8.276e+02 9.800e-04 +8.286e+02 9.779e-04 +8.296e+02 9.759e-04 +8.305e+02 9.738e-04 +8.315e+02 9.717e-04 +8.325e+02 9.697e-04 +8.335e+02 9.677e-04 +8.344e+02 9.656e-04 +8.354e+02 9.636e-04 +8.364e+02 9.616e-04 +8.374e+02 9.596e-04 +8.383e+02 9.576e-04 +8.393e+02 9.556e-04 +8.403e+02 9.536e-04 +8.413e+02 9.516e-04 +8.422e+02 9.496e-04 +8.432e+02 9.476e-04 +8.442e+02 9.457e-04 +8.452e+02 9.437e-04 +8.461e+02 9.418e-04 +8.471e+02 9.398e-04 +8.481e+02 9.379e-04 +8.490e+02 9.359e-04 +8.500e+02 9.340e-04 +8.510e+02 9.321e-04 +8.520e+02 9.302e-04 +8.529e+02 9.283e-04 +8.539e+02 9.264e-04 +8.549e+02 9.245e-04 +8.559e+02 9.226e-04 +8.568e+02 9.207e-04 +8.578e+02 9.188e-04 +8.588e+02 9.169e-04 +8.598e+02 9.150e-04 +8.607e+02 9.132e-04 +8.617e+02 9.113e-04 +8.627e+02 9.095e-04 +8.637e+02 9.076e-04 +8.646e+02 9.058e-04 +8.656e+02 9.040e-04 +8.666e+02 9.021e-04 +8.676e+02 9.003e-04 +8.685e+02 8.985e-04 +8.695e+02 8.967e-04 +8.705e+02 8.949e-04 +8.714e+02 8.931e-04 +8.724e+02 8.913e-04 +8.734e+02 8.895e-04 +8.744e+02 8.877e-04 +8.753e+02 8.859e-04 +8.763e+02 8.842e-04 +8.773e+02 8.824e-04 +8.783e+02 8.806e-04 +8.792e+02 8.789e-04 +8.802e+02 8.771e-04 +8.812e+02 8.754e-04 +8.822e+02 8.737e-04 +8.831e+02 8.719e-04 +8.841e+02 8.702e-04 +8.851e+02 8.685e-04 +8.861e+02 8.668e-04 +8.870e+02 8.650e-04 +8.880e+02 8.633e-04 +8.890e+02 8.616e-04 +8.900e+02 8.599e-04 +8.909e+02 8.582e-04 +8.919e+02 8.566e-04 +8.929e+02 8.549e-04 +8.938e+02 8.532e-04 +8.948e+02 8.515e-04 +8.958e+02 8.499e-04 +8.968e+02 8.482e-04 +8.977e+02 8.465e-04 +8.987e+02 8.449e-04 +8.997e+02 8.433e-04 +9.007e+02 8.416e-04 +9.016e+02 8.400e-04 +9.026e+02 8.383e-04 +9.036e+02 8.367e-04 +9.046e+02 8.351e-04 +9.055e+02 8.335e-04 +9.065e+02 8.319e-04 +9.075e+02 8.303e-04 +9.085e+02 8.287e-04 +9.094e+02 8.271e-04 +9.104e+02 8.255e-04 +9.114e+02 8.239e-04 +9.124e+02 8.223e-04 +9.133e+02 8.207e-04 +9.143e+02 8.192e-04 +9.153e+02 8.176e-04 +9.162e+02 8.160e-04 +9.172e+02 8.145e-04 +9.182e+02 8.129e-04 +9.192e+02 8.114e-04 +9.201e+02 8.098e-04 +9.211e+02 8.083e-04 +9.221e+02 8.067e-04 +9.231e+02 8.052e-04 +9.240e+02 8.037e-04 +9.250e+02 8.022e-04 +9.260e+02 8.006e-04 +9.270e+02 7.991e-04 +9.279e+02 7.976e-04 +9.289e+02 7.961e-04 +9.299e+02 7.946e-04 +9.309e+02 7.931e-04 +9.318e+02 7.916e-04 +9.328e+02 7.901e-04 +9.338e+02 7.887e-04 +9.348e+02 7.872e-04 +9.357e+02 7.857e-04 +9.367e+02 7.842e-04 +9.377e+02 7.828e-04 +9.386e+02 7.813e-04 +9.396e+02 7.799e-04 +9.406e+02 7.784e-04 +9.416e+02 7.769e-04 +9.425e+02 7.755e-04 +9.435e+02 7.741e-04 +9.445e+02 7.726e-04 +9.455e+02 7.712e-04 +9.464e+02 7.698e-04 +9.474e+02 7.683e-04 +9.484e+02 7.669e-04 +9.494e+02 7.655e-04 +9.503e+02 7.641e-04 +9.513e+02 7.627e-04 +9.523e+02 7.613e-04 +9.533e+02 7.599e-04 +9.542e+02 7.585e-04 +9.552e+02 7.571e-04 +9.562e+02 7.557e-04 +9.571e+02 7.543e-04 +9.581e+02 7.530e-04 +9.591e+02 7.516e-04 +9.601e+02 7.502e-04 +9.610e+02 7.488e-04 +9.620e+02 7.475e-04 +9.630e+02 7.461e-04 +9.640e+02 7.448e-04 +9.649e+02 7.434e-04 +9.659e+02 7.421e-04 +9.669e+02 7.407e-04 +9.679e+02 7.394e-04 +9.688e+02 7.380e-04 +9.698e+02 7.367e-04 +9.708e+02 7.354e-04 +9.718e+02 7.340e-04 +9.727e+02 7.327e-04 +9.737e+02 7.314e-04 +9.747e+02 7.301e-04 +9.757e+02 7.288e-04 +9.766e+02 7.275e-04 +9.776e+02 7.262e-04 +9.786e+02 7.249e-04 +9.795e+02 7.236e-04 +9.805e+02 7.223e-04 +9.815e+02 7.210e-04 +9.825e+02 7.197e-04 +9.834e+02 7.184e-04 +9.844e+02 7.171e-04 +9.854e+02 7.159e-04 +9.864e+02 7.146e-04 +9.873e+02 7.133e-04 +9.883e+02 7.121e-04 +9.893e+02 7.108e-04 +9.903e+02 7.095e-04 +9.912e+02 7.083e-04 +9.922e+02 7.070e-04 +9.932e+02 7.058e-04 +9.942e+02 7.045e-04 +9.951e+02 7.033e-04 +9.961e+02 7.021e-04 +9.971e+02 7.008e-04 +9.981e+02 6.996e-04 +9.990e+02 6.984e-04 +1.000e+03 6.971e-04 \ No newline at end of file diff --git a/sigproc/lmfit/minimizer.py b/sigproc/lmfit/minimizer.py index 1713de241..f4ed3b36a 100644 --- a/sigproc/lmfit/minimizer.py +++ b/sigproc/lmfit/minimizer.py @@ -18,7 +18,7 @@ from scipy.optimize import leastsq as scipy_leastsq from scipy.optimize import fmin as scipy_fmin -from scipy.optimize import anneal as scipy_anneal +#from scipy.optimize import anneal as scipy_anneal from scipy.optimize.lbfgsb import fmin_l_bfgs_b as scipy_lbfgsb # check for scipy.optimize.minimize @@ -279,7 +279,7 @@ def anneal(self, schedule='cauchy', **kws): sakws.update(self.kws) sakws.update(kws) print("WARNING: scipy anneal appears unusable!") - saout = scipy_anneal(self.penalty, self.vars, **sakws) + #saout = scipy_anneal(self.penalty, self.vars, **sakws) self.sa_out = saout return diff --git a/sigproc/lmfit/minimizer.py~ b/sigproc/lmfit/minimizer.py~ deleted file mode 100644 index 62fc85660..000000000 --- a/sigproc/lmfit/minimizer.py~ +++ /dev/null @@ -1,429 +0,0 @@ -""" -Simple minimizer is a wrapper around scipy.leastsq, allowing a -user to build a fitting model as a function of general purpose -Fit Parameters that can be fixed or floated, bounded, and written -as a simple expression of other Fit Parameters. - -The user sets up a model in terms of instance of Parameters, writes a -function-to-be-minimized (residual function) in terms of these Parameters. - - Copyright (c) 2011 Matthew Newville, The University of Chicago - -""" - -from numpy import sqrt - -from scipy.optimize import leastsq as scipy_leastsq -from scipy.optimize import anneal as scipy_anneal -from scipy.optimize import fmin_l_bfgs_b as scipy_lbfgsb - -try: - from collections import OrderedDict -except ImportError: - from ordereddict import OrderedDict - -from .asteval import Interpreter -from .astutils import NameFinder, valid_symbol_name - -class Parameters(OrderedDict): - """a custom dictionary of Parameters. All keys must be - strings, and valid Python symbol names, and all values - must be Parameters. - - Custom methods: - --------------- - - add() - add_many() - """ - def __init__(self, *args, **kwds): - OrderedDict.__init__(self) - self.update(*args, **kwds) - - def __setitem__(self, key, value): - if key not in self: - if not valid_symbol_name(key): - raise KeyError("'%s' is not a valid Parameters name" % key) - if value is not None and not isinstance(value, Parameter): - raise ValueError("'%s' is not a Parameter" % value) - OrderedDict.__setitem__(self, key, value) - value.name = key - - def add(self, name, value=None, vary=True, expr=None, - min=None, max=None): - """convenience function for adding a Parameter: - with p = Parameters() - p.add(name, value=XX, ....) - - is equivalent to - p[name] = Parameter(name=name, value=XX, .... - """ - self.__setitem__(name, Parameter(value=value, name=name, vary=vary, - expr=expr, min=min, max=max)) - - def add_many(self, *parlist): - """convenience function for adding a list of Parameters: - Here, you must provide a sequence of tuples, each containing - at least the name. The order in each tuple is the following: - name, value, vary, min, max, expr - with p = Parameters() - p.add_many( (name1, val1, True, None, None, None), - (name2, val2, True, 0.0, None, None), - (name3, val3, False, None, None, None), - (name4, val4)) - - """ - for para in parlist: - self.add(*para) - -class Parameter(object): - """A Parameter is the basic Parameter going - into Fit Model. The Parameter holds many attributes: - value, vary, max_value, min_value, constraint expression. - The value and min/max values will be be set to floats. - """ - def __init__(self, name=None, value=None, vary=True, - min=None, max=None, expr=None, **kws): - self.name = name - self.value = value - self.init_value = value - self.min = min - self.max = max - self.vary = vary - self.expr = expr - self.stderr = None - self.correl = None - - def __repr__(self): - s = [] - if self.name is not None: - s.append("'%s'" % self.name) - val = repr(self.value) - if self.vary and self.stderr is not None: - val = "value=%s +/- %.3g" % (repr(self.value), self.stderr) - elif not self.vary: - val = "value=%s (fixed)" % (repr(self.value)) - s.append(val) - s.append("bounds=[%s:%s]" % (repr(self.min), repr(self.max))) - if self.expr is not None: - s.append("expr='%s'" % (self.expr)) - return "" % ', '.join(s) - -class MinimizerException(Exception): - """General Purpose Exception""" - def __init__(self, msg): - Exception.__init__(self) - self.msg = msg - - def __str__(self): - return "\n%s" % (self.msg) - -def check_ast_errors(error): - """check for errors derived from asteval, raise MinimizerException""" - if len(error) > 0: - msg = [] - for err in error: - msg = '\n'.join(err.get_error()) - raise MinimizerException(msg) - - -class Minimizer(object): - """general minimizer""" - err_nonparam = \ - "params must be a minimizer.Parameters() instance or list of Parameters()" - err_maxfev = """Too many function calls (max set to %i)! Use: - minimize(func, params, ...., maxfev=NNN) -or set leastsq_kws['maxfev'] to increase this maximum.""" - - def __init__(self, userfcn, params, fcn_args=None, fcn_kws=None, - iter_cb=None, scale_covar=True, **kws): - self.userfcn = userfcn - print type(params) - self.__set_params(params) - self.userargs = fcn_args - if self.userargs is None: - self.userargs = [] - - self.userkws = fcn_kws - if self.userkws is None: - self.userkws = {} - self.kws = kws - self.iter_cb = iter_cb - self.scale_covar = scale_covar - self.nfev = 0 - self.nfree = 0 - self.message = None - self.var_map = [] - self.jacfcn = None - self.asteval = Interpreter() - self.namefinder = NameFinder() - self.__prepared = False - - def __update_paramval(self, name): - """ - update parameter value, including setting bounds. - For a constrained parameter (one with an expr defined), - this first updates (recursively) all parameters on which - the parameter depends (using the 'deps' field). - """ - # Has this param already been updated? - # if this is called as an expression dependency, - # it may have been! - if self.updated[name]: - return - - par = self.params[name] - val = par.value - if par.expr is not None: - for dep in par.deps: - self.__update_paramval(dep) - val = self.asteval.run(par.ast) - check_ast_errors(self.asteval.error) - # apply min/max - if par.min is not None: - val = max(val, par.min) - if par.max is not None: - val = min(val, par.max) - - self.asteval.symtable[name] = par.value = float(val) - self.updated[name] = True - - def __residual(self, fvars): - """ - residual function used for least-squares fit. - With the new, candidate values of fvars (the fitting variables), - this evaluates all parameters, including setting bounds and - evaluating constraints, and then passes those to the - user-supplied function to calculate the residual. - """ - # set parameter values - for varname, val in zip(self.var_map, fvars): - self.params[varname].value = val - self.nfev = self.nfev + 1 - - self.updated = dict([(name, False) for name in self.params]) - for name in self.params: - self.__update_paramval(name) - - out = self.userfcn(self.params, *self.userargs, **self.userkws) - if hasattr(self.iter_cb, '__call__'): - self.iter_cb(self.params, self.nfev, out, - *self.userargs, **self.userkws) - return out - - def __jacobian(self, fvars): - """ - analytical jacobian to be used with the Levenberg-Marquardt - - modified 02-01-2012 by Glenn Jones, Aberystwyth University - """ - for varname, val in zip(self.var_map, fvars): - self.params[varname].value = val - self.nfev = self.nfev + 1 - - self.updated = dict([(name, False) for name in self.params]) - for name in self.params: - self.__update_paramval(name) - - # computing the jacobian - return self.jacfcn(self.params, *self.userargs, **self.userkws) - - def __set_params(self, params): - """ set internal self.params from a Parameters object or - a list/tuple of Parameters""" - if params is None or isinstance(params, Parameters): - self.params = params - elif isinstance(params, (list, tuple)): - _params = Parameters() - for _par in params: - if not isinstance(_par, Parameter): - raise MinimizerException(self.err_nonparam) - else: - _params[_par.name] = _par - self.params = _params - else: - raise MinimizerException(self.err_nonparam) - - def prepare_fit(self, params=None): - """prepare parameters for fit""" - # determine which parameters are actually variables - # and which are defined expressions. - if params is None and self.params is not None and self.__prepared: - return - if params is not None and self.params is None: - self.__set_params(params) - self.nfev = 0 - self.var_map = [] - self.vars = [] - self.vmin = [] - self.vmax = [] - for name, par in self.params.items(): - if par.expr is not None: - par.ast = self.asteval.parse(par.expr) - check_ast_errors(self.asteval.error) - par.vary = False - par.deps = [] - self.namefinder.names = [] - self.namefinder.generic_visit(par.ast) - for symname in self.namefinder.names: - if (symname in self.params and - symname not in par.deps): - par.deps.append(symname) - elif par.vary: - self.var_map.append(name) - self.vars.append(par.value) - self.vmin.append(par.min) - self.vmax.append(par.max) - - self.asteval.symtable[name] = par.value - par.init_value = par.value - if par.name is None: - par.name = name - - self.nvarys = len(self.vars) - - # now evaluate make sure initial values - # are used to set values of the defined expressions. - # this also acts as a check of expression syntax. - self.updated = dict([(name, False) for name in self.params]) - for name in self.params: - self.__update_paramval(name) - self.__prepared = True - - def anneal(self, schedule='cauchy', **kws): - """ - use simulated annealing - """ - sched = 'fast' - if schedule in ('cauchy', 'boltzmann'): - sched = schedule - - self.prepare_fit() - sakws = dict(full_output=1, schedule=sched, - maxiter = 2000 * (self.nvarys + 1)) - - sakws.update(self.kws) - sakws.update(kws) - - def penalty(params): - "local penalty function -- anneal wants sum-squares residual" - r = self.__residual(params) - return (r*r).sum() - - saout = scipy_anneal(penalty, self.vars, **sakws) - self.sa_out = saout - return - - def lbfgsb(self, **kws): - """ - use l-bfgs-b minimization - """ - self.prepare_fit() - lb_kws = dict(factr=1000.0, approx_grad=True, m=20, - maxfun = 2000 * (self.nvarys + 1), - bounds = zip(self.vmin, self.vmax)) - lb_kws.update(self.kws) - lb_kws.update(kws) - def penalty(params): - "local penalty function -- lbgfsb wants sum-squares residual" - r = self.__residual(params) - return (r*r).sum() - - xout, fout, info = scipy_lbfgsb(penalty, self.vars, **lb_kws) - - self.nfev = info['funcalls'] - self.message = info['task'] - - - def leastsq(self, scale_covar=True, **kws): - """ - use Levenberg-Marquardt minimization to perform fit. - This assumes that ModelParameters have been stored, - and a function to minimize has been properly set up. - - This wraps scipy.optimize.leastsq, and keyward arguments are passed - directly as options to scipy.optimize.leastsq - - When possible, this calculates the estimated uncertainties and - variable correlations from the covariance matrix. - - writes outputs to many internal attributes, and - returns True if fit was successful, False if not. - """ - self.prepare_fit() - lskws = dict(full_output=1, xtol=1.e-7, ftol=1.e-7, - gtol=1.e-7, maxfev=1000*(self.nvarys+1), Dfun=None) - - lskws.update(self.kws) - lskws.update(kws) - - if lskws['Dfun'] is not None: - self.jacfcn = lskws['Dfun'] - lskws['Dfun'] = self.__jacobian - - lsout = scipy_leastsq(self.__residual, self.vars, **lskws) - vbest, cov, infodict, errmsg, ier = lsout - - self.residual = resid = infodict['fvec'] - - self.ier = ier - self.lmdif_message = errmsg - self.message = 'Fit succeeded.' - self.success = ier in [1, 2, 3, 4] - - if ier == 0: - self.message = 'Invalid Input Parameters.' - elif ier == 5: - self.message = self.err_maxfev % lskws['maxfev'] - else: - self.message = 'Tolerance seems to be too small.' - - self.nfev = infodict['nfev'] - self.ndata = len(resid) - - sum_sqr = (resid**2).sum() - self.chisqr = sum_sqr - self.nfree = (self.ndata - self.nvarys) - self.redchi = sum_sqr / self.nfree - - for par in self.params.values(): - par.stderr = 0 - par.correl = None - if hasattr(par, 'ast'): - delattr(par, 'ast') - - if cov is None: - self.errorbars = False - self.message = '%s. Could not estimate error-bars' - else: - self.errorbars = True - self.covar = cov - if self.scale_covar: - cov = cov * sum_sqr / self.nfree - for ivar, varname in enumerate(self.var_map): - par = self.params[varname] - par.stderr = sqrt(cov[ivar, ivar]) - par.correl = {} - for jvar, varn2 in enumerate(self.var_map): - if jvar != ivar: - par.correl[varn2] = (cov[ivar, jvar]/ - (par.stderr * sqrt(cov[jvar, jvar]))) - - return self.success - -def minimize(fcn, params, engine='leastsq', args=None, kws=None, - scale_covar=True, - iter_cb=None, **fit_kws): - """simple minimization function, - finding the values for the params which give the - minimal sum-of-squares of the array return by fcn - """ - fitter = Minimizer(fcn, params, fcn_args=args, fcn_kws=kws, - iter_cb=iter_cb, scale_covar=scale_covar, **fit_kws) - if engine == 'anneal': - fitter.anneal() - elif engine == 'lbfgsb': - fitter.lbfgsb() - else: - fitter.leastsq() - return fitter