Python Planck13.age方法代碼示例

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import age [as 別名]
def z_to_time(z,cosmology=None):
    """
    gives the time, in years, since the big bang given an astropy cosmology
    and a redshift.  if no cosmology passed in, assumes Planck13
    """
    if cosmology is None:
        from astropy.cosmology import Planck13 as cosmology
    from yt import YTArray
    t = cosmology.age(z)
    return YTArray(t*t.unit.in_units('yr'),'yr')

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import age [as 別名]
t=5*u.s
v=s/t
print 'v= ',v,'\n'

cms=u.cm/u.s
from astropy.units import imperial
mph=imperial.mile/u.h

from astropy.cosmology import Planck13 #imports 2013 results of Planck
#from astropy import cosmology #imporst cosmology
#cosmology.core.set_current(cosmology.core.Planck13) #sets current cosmology model
from math import pi

print Planck13.H(0) #gives H at certain z
print Planck13.age(0)

from astropy.cosmology import FlatLambdaCDM
cosmo=FlatLambdaCDM(H0=70, Om0=.3)
print cosmo.H(0)


# read/write files with astropy
import numpy as np
from astropy.table import table

x=np.asarray([1,6,5,3,29])
y_measurements=np.asarray([10,20,30,40,50])

measurements=table.Table([x,y_measurements],names=('x','y'))
print measurements #saving this as a file will organize it as a table

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import age [as 別名]
def firefly_single(parameters):

	"""
	The routine for a single run of FIREFLY.
	It is called from firefly_job, test_firefly, or can be run 
	interactively for a custom single SED.
	
	In the interactive case, one sets the 'custom' and 'interactive' parameters
	in the parameter file, then one enters at an interative prompt:
	> from firefly_single import firefly_single
	> firefly_single('[locationpathofdata]','custom/[outputdirname]','./parameters.dat')
	One is then able to view the output plots at the native X window. 

	This routine retrieves the options from the parameters file,
	including the locations of the data file and models to be used.

	It then opens the data file, model files, matches their resolutions
	(downgrading to instrumental+velocity disp resolution if necessary)
	fits, then produces output files and plots. 

	INPUTS: 
	- options_file: location of the parameter file (default: ./parameters.dat) 

	No outputs.
	"""

    # data_file, output_dir
	

	data = get_data(parameters)

	# restrict_ages can be default (allow age<+1 Gyr age uni), off (allow all ages), strict (age<age uni only)
	if parameters['restrict_ages'] == 'default':
		age_universe = Planck13.age(data['redshift'])
		parameters['age_limits'][1] = np.log10(age_universe.value+1.0)+9.0 # log(yr units)
	elif parameters['restrict_ages'] == 'strict':
		age_universe = Planck13.age(data['redshift'])
		parameters['age_limits'][1] = np.log10(age_universe.value)+9.0

	# Get the models with observation information needed for downgrading.
	for mi,mm in enumerate(parameters['model_libs']):
		for ii in parameters['imfs']:
			deltal = parameters['deltal_libs'][mi]
			model_wave_int,model_flux_int, age,metal = \
				get_model(parameters,mm,ii,deltal,data['vdisp'],data['wavelength'],data['r_instrument'],data['ebv_mw'])

			print "Matching data to models..."
			wave,data_flux,error_flux,model_flux_raw = \
				match_data_models(data['wavelength'],data['flux'],data['flags'],data['error'],model_wave_int,model_flux_int,\
					parameters['wave_limits'][0],parameters['wave_limits'][1])
			print "Normalising all spectra."
			model_flux,mass_factors = normalise_spec(data_flux,model_flux_raw)

			# Get filtered values IF dust is on!
			if parameters['hpf_mode']=='on':

				print "Determining attenuation curve through HPF fitting"
				best_ebv,attenuation_curve = determine_attenuation(wave,data_flux,error_flux,model_flux,parameters,age,metal)

				if parameters['plot_diagnostics']:
					print "Best ebv is "+str(best_ebv)
					plt.plot(attenuation_curve)
					plt.title("Attenuation curve")
					plt.show()


				# Apply curve to models and renormalise:
				print "Curve found! Applying to models..."
				model_flux_atten = np.zeros(np.shape(model_flux_raw))
				for m in range(len(model_flux_raw)):
					model_flux_atten[m] 	= attenuation_curve * model_flux_raw[m]
				model_flux,mass_factors = normalise_spec(data_flux,model_flux_atten)
				print "Fitting with attenuated models..."
				light_weights,chis,branch = fitter(wave,data_flux,error_flux,model_flux,parameters)
			
			elif parameters['hpf_mode'] == 'hpf_only':

				print "Using filtered values to determing SP properties only."
				smoothing_length = parameters['dust_smoothing_length']
				hpf_data    = hpf(data_flux)
				hpf_models  = np.zeros(np.shape(model_flux))
				for m in range(len(model_flux)):
					hpf_models[m] = hpf(model_flux[m])


				zero_dat = np.where( (np.isnan(hpf_data)) & (np.isinf(hpf_data)) )
				hpf_data[zero_dat] = 0.0
				for m in range(len(model_flux)):
					hpf_models[m,zero_dat] = 0.0
				hpf_error    = np.zeros(len(error_flux))
				hpf_error[:] = np.median(error_flux)/np.median(data_flux) * np.median(hpf_data)
				hpf_error[zero_dat] = np.max(hpf_error)*999999.9

				best_ebv = -99
				hpf_models,mass_factors = normalise_spec(hpf_data,hpf_models)
				light_weights,chis,branch = fitter(wave,hpf_data,hpf_error,hpf_models,parameters)

			elif parameters['hpf_mode'] == 'off':
				raise NotImplementedError("Not using a HPF and fitting using model curves not implemented yet")
				# use loop over dust curve, but this will take a while!
#.........這裏部分代碼省略.........

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import age [as 別名]
def Semimajor_to_Period(A,M1,M2):
    Mtot = M1 + M2
    from numpy import sqrt,pi
    return sqrt((4*pi*pi) * A**3 / (phys_const.G*Mtot))

def Semimajor_to_Coalescence(a,M1,M2):
    mu = (M1*M2)/(M1+M2)
    Mtot = M1+M2
    return 5./256.*(phys_const.clight**5 * a**4)/(phys_const.G**3 * Mtot**2 * mu)

def Period_to_Coalescence(P,M1,M2):
    a = Period_to_Semimajor(P,M1,M2)
    return Semimajor_to_Coalescence(a,M1,M2)

th = cosmo.age(0)
thubble = YTQuantity(th*th.unit.in_units('yr'),'yr')

M1 = 30. * units.Msun
M2 = 30. * units.Msun
Mtot = M1 + M2
# mu = (M**2/Mtot)

Pmin = 2.*units.day
amin = Period_to_Semimajor(Pmin,M1,M2)
tmin = Semimajor_to_Coalescence(amin,M1,M2)
# amin = (((phys_const.G*Mtot)/(4.*pi*pi)) * Pmin**2)**(1./3.) #.in_units('kpc')
# tmin = 5./256.*(phys_const.clight**5 * amin**4)/(phys_const.G**3 * Mtot**2 * mu)

Pmax = 20.*units.day
amax = Period_to_Semimajor(Pmax,M1,M2)

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import age [as 別名]
def firefly_single(parameters):

    """
	The routine for a single run of FIREFLY.
	It is called from firefly_job, test_firefly, or can be run 
	interactively for a custom single SED.
	
	In the interactive case, one sets the 'custom' and 'interactive' parameters
	in the parameter file, then one enters at an interative prompt:
	> from firefly_single import firefly_single
	> firefly_single('[locationpathofdata]','custom/[outputdirname]','./parameters.dat')
	One is then able to view the output plots at the native X window. 

	This routine retrieves the options from the parameters file,
	including the locations of the data file and models to be used.

	It then opens the data file, model files, matches their resolutions
	(downgrading to instrumental+velocity disp resolution if necessary)
	fits, then produces output files and plots. 

	INPUTS: 
	- options_file: location of the parameter file (default: ./parameters.dat) 

	No outputs.
	"""

    # data_file, output_dir

    data = get_data(parameters)

    # restrict_ages can be default (allow age<+1 Gyr age uni), off (allow all ages), strict (age<age uni only)
    if parameters["restrict_ages"] == "default":
        age_universe = Planck13.age(data["redshift"])
        parameters["age_limits"][1] = np.log10(age_universe.value + 1.0) + 9.0  # log(yr units)
    elif parameters["restrict_ages"] == "strict":
        age_universe = Planck13.age(data["redshift"])
        parameters["age_limits"][1] = np.log10(age_universe.value) + 9.0

        # Get the models with observation information needed for downgrading.
    for mi, mm in enumerate(parameters["model_libs"]):
        for ii in parameters["imfs"]:
            deltal = parameters["deltal_libs"][mi]
            model_wave_int, model_flux_int, age, metal = get_model(
                parameters, mm, ii, deltal, data["vdisp"], data["wavelength"], data["r_instrument"], data["ebv_mw"]
            )

            print "Matching data to models..."
            wave, data_flux, error_flux, model_flux_raw = match_data_models(
                data["wavelength"],
                data["flux"],
                data["flags"],
                data["error"],
                model_wave_int,
                model_flux_int,
                parameters["wave_limits"][0],
                parameters["wave_limits"][1],
            )
            print "Normalising all spectra."
            model_flux, mass_factors = normalise_spec(data_flux, model_flux_raw)

            # Get filtered values IF dust is on!
            if parameters["hpf_mode"] == "on":

                print "Determining attenuation curve through HPF fitting"
                best_ebv, attenuation_curve = determine_attenuation(
                    wave, data_flux, error_flux, model_flux, parameters, age, metal
                )

                if parameters["plot_diagnostics"]:
                    print "Best ebv is " + str(best_ebv)
                    plt.plot(attenuation_curve)
                    plt.title("Attenuation curve")
                    plt.show()

                    # Apply curve to models and renormalise:
                print "Curve found! Applying to models..."
                model_flux_atten = np.zeros(np.shape(model_flux_raw))
                for m in range(len(model_flux_raw)):
                    model_flux_atten[m] = attenuation_curve * model_flux_raw[m]
                model_flux, mass_factors = normalise_spec(data_flux, model_flux_atten)
                print "Fitting with attenuated models..."
                light_weights_int, chis_int, branch = fitter(wave, data_flux, error_flux, model_flux, parameters)

            elif parameters["hpf_mode"] == "hpf_only":

                print "Using filtered values to determing SP properties only."
                smoothing_length = parameters["dust_smoothing_length"]
                hpf_data = hpf(data_flux)
                hpf_models = np.zeros(np.shape(model_flux))
                for m in range(len(model_flux)):
                    hpf_models[m] = hpf(model_flux[m])

                zero_dat = np.where((np.isnan(hpf_data)) & (np.isinf(hpf_data)))
                hpf_data[zero_dat] = 0.0
                for m in range(len(model_flux)):
                    hpf_models[m, zero_dat] = 0.0
                hpf_error = np.zeros(len(error_flux))
                hpf_error[:] = np.median(error_flux) / np.median(data_flux) * np.median(hpf_data)
                hpf_error[zero_dat] = np.max(hpf_error) * 999999.9

#.........這裏部分代碼省略.........