Python WMAP9.lookback_time方法代码示例

# 需要导入模块: from astropy.cosmology import WMAP9 [as 别名]
# 或者: from astropy.cosmology.WMAP9 import lookback_time [as 别名]
def cosmoLookBack(redshift,
                  WMAP9=False,
                  H0=70.0,
                  Om0=0.30,
                  Planck15=False,
                  Myr=False):
    """
    Get the Look-back Time at redshift=z.

    This is simply a wrapper of astropy.cosmology
    The input redsfhit can be an array
    """
    if WMAP9:
        from astropy.cosmology import WMAP9 as cosmo
    elif Planck15:
        from astropy.cosmology import Planck15 as cosmo
    else:
        from astropy.cosmology import FlatLambdaCDM
        cosmo = FlatLambdaCDM(H0=H0, Om0=Om0)

    lbt = cosmo.lookback_time(redshift)

    if not Myr:
        return lbt.value
    else:
        return lbt.to(u.Myr).value

# 需要导入模块: from astropy.cosmology import WMAP9 [as 别名]
# 或者: from astropy.cosmology.WMAP9 import lookback_time [as 别名]
                            chival = np.sum((obj_fluxes/obj_fluxerrs - const*th_fluxes/obj_fluxerrs)**2)
                            print chival
                            raw_input()
                            if chival < output[m,7]:
                                output[m,:] = np.array([m+1, all_obj_specz[m], z, ages[j], float(tauvals[l]), EBV, const, chival])
                print "Age: " + str(ages[j]) + ", tau: " + tauvals[l] + ", redshift: " + str(z)
                np.savetxt("photoz_grid84876.txt", output, header="obj_no spec_z phot_z age tau EBV norm chi")

"""
#the chi squared array is now 2D over EBV and the different objects

for j in range(128): 
    for l in range(3):
        th_mag_array = np.loadtxt("synmags_T" + tauvals[l] + "/synmags/synmags_age_" + str(ages[j]) + ".txt", usecols=(1,2,3,4,5,6,7,8,9,10,11,12))
        for i in range(1, 501): 
            z = 0.01*i
            th_mags = th_mag_array[i-1,:]
            if ages[j]*(10**-9) < 14.00 - WMAP9.lookback_time(z).value:
                for k in range(41):
                    EBV = 0.05*k
                    th_fluxes = 10**((23.9 - EBV*coef - th_mags)/2.5) #microjanskys
                    const =  np.sum(fulldet_obj_fluxes*th_fluxes/fulldet_obj_fluxerrs**2, axis=1)/np.sum(th_fluxes**2/fulldet_obj_fluxerrs**2, axis=1)
                    chival = np.sum((fulldet_obj_fluxes/fulldet_obj_fluxerrs - (const*(th_fluxes/fulldet_obj_fluxerrs).T).T)**2, axis=1)
                    for m in range(len(fulldet_obj_fluxes)):
                        if chival[m] < output[m,7]:
                            output[m,:] = np.array([m+1, all_obj_specz[m], z, ages[j], float(tauvals[l]), EBV, const[m], chival[m]])
                print "Age: " + str(ages[j]) + ", tau: " + tauvals[l] + ", redshift: " + str(z)
                np.savetxt("photoz_grid86548.txt", output, header="obj_no spec_z phot_z age tau EBV norm chi")

# 需要导入模块: from astropy.cosmology import WMAP9 [as 别名]
# 或者: from astropy.cosmology.WMAP9 import lookback_time [as 别名]
        offset_errs[0, i, 0] = 0.

all_obj_fluxes = all_obj_fluxes + offsets
all_obj_fluxerrs = np.sqrt(all_obj_fluxerrs**2 + offset_errs**2)


### Build array for photometric redshift and parameter outputs
output = np.zeros(len(all_obj_fluxes)*8, dtype="float")
output.shape = (len(all_obj_fluxes), 8)
output[:,7] = 9999999999.


### Perform model fitting using a 2D chi squared array over object and redshift with max age of stellar pop physically determined
tauvals = np.array([0.05, 1, 10], dtype="str")
zarr = np.arange(0.01, 5.001, 0.01)
lbtarr = WMAP9.lookback_time(zarr).value
EBVarr = np.expand_dims(np.expand_dims(np.expand_dims(np.arange(0, 1.50001, 0.025), axis=0), axis=0), axis=0)
    
for j in range(minage, maxage+1): 
    arg = 0
    while ages[j]*(10**-9) < 14.0 - lbtarr[arg] and arg < 499:
        arg = arg+1
    for l in range(3):
        print "Age: " + str(ages[j]) + ", fitting to redshift " + str((arg+1)*0.01) + ", with tau: " + tauvals[l]
        th_mag_array = np.expand_dims(np.expand_dims(np.loadtxt("models/" + tauvals[l] + "/synmags_age_" + str(ages[j]) + ".txt", usecols=(1,2,3,4,5,6,7,8,9,10,11,12))[0:arg, :].T, axis=0), axis=3)
        print "th_mag_array computed"
        print th_mag_array.nbytes
        th_flux_array = (10**((23.9 - EBVarr*coef - th_mag_array)/2.5)) #microjanskys
        print "th_flux_array computed"
        print th_flux_array.nbytes
        const =  np.expand_dims(np.sum(all_obj_fluxes*th_flux_array/all_obj_fluxerrs**2, axis=1)/np.sum(th_flux_array**2/all_obj_fluxerrs**2, axis=1), axis=1)

# 需要导入模块: from astropy.cosmology import WMAP9 [as 别名]
# 或者: from astropy.cosmology.WMAP9 import lookback_time [as 别名]
    filter_bounds, filter_widths = make_bins(filter[:,0], make_rhs="True")
    top = np.zeros(len(spec_widths)*len(filter_widths))
    bottom = np.copy(top)

    top = (np.expand_dims(filter_bounds[1:], axis=1) - np.expand_dims(spec_bounds[:-1], axis=0))/np.expand_dims(filter_widths, axis=1)
    bottom = (np.expand_dims(spec_bounds[1:], axis=0) - np.expand_dims(filter_bounds[:-1], axis=1))/np.expand_dims(filter_widths, axis=1)
    part1 = np.logical_and(top>=1, bottom>=1).astype("float")
    part2 = np.logical_and(np.logical_and(0<top, top<1), bottom>=1).astype("float")
    part3 = np.logical_and(np.logical_and(0<bottom, bottom<1), top>=1).astype("float")
    part4 = np.logical_and(np.logical_and(0<top, top<1), np.logical_and(0<bottom, bottom<1)).astype("float")
    crossover = part1 + top*part2 + bottom*part3 + (top+bottom-1)*part4 
    flux_total = np.sum(np.sum(crossover*np.expand_dims(filter[:,1], axis=1)*np.expand_dims(spec[:,1], axis=0)*np.expand_dims(filter_widths, axis=1), axis=1), axis=0)
    return flux_total
    
zarr = np.arange(0.01, 7.001, 0.01)
lbtarr = cosmo.lookback_time(zarr).value
D = np.loadtxt("../IGM_Da_Db.txt")

### Load up the data necessary to build the "filter" to rebin the model spectra onto

hdulist1 = fits.open("../../VANDELS_data/spectra/sc_206806_UDS_P1M1_MR_Q1_029_1.fits")
wavzpt = hdulist1[0].header["CRVAL1"]
dwav = hdulist1[0].header["CDELT1"]
fluxes1 = hdulist1[4].data

maxwav = wavzpt + dwav*(len(fluxes1))

objwavs_nobin = np.arange(wavzpt, maxwav, dwav)

objwavs = specbin(objwavs_nobin, 2)

# 需要导入模块: from astropy.cosmology import WMAP9 [as 别名]
# 或者: from astropy.cosmology.WMAP9 import lookback_time [as 别名]
    filter_bounds, filter_widths = make_bins(filter[:,0], make_rhs="True")
    top = np.zeros(len(spec_widths)*len(filter_widths))
    bottom = np.copy(top)

    top = (np.expand_dims(filter_bounds[1:], axis=1) - np.expand_dims(spec_bounds[:-1], axis=0))/np.expand_dims(filter_widths, axis=1)
    bottom = (np.expand_dims(spec_bounds[1:], axis=0) - np.expand_dims(filter_bounds[:-1], axis=1))/np.expand_dims(filter_widths, axis=1)
    part1 = np.logical_and(top>=1, bottom>=1).astype("float")
    part2 = np.logical_and(np.logical_and(0<top, top<1), bottom>=1).astype("float")
    part3 = np.logical_and(np.logical_and(0<bottom, bottom<1), top>=1).astype("float")
    part4 = np.logical_and(np.logical_and(0<top, top<1), np.logical_and(0<bottom, bottom<1)).astype("float")
    crossover = part1 + top*part2 + bottom*part3 + (top+bottom-1)*part4 
    flux_total = np.sum(np.sum(crossover*np.expand_dims(filter[:,1], axis=1)*np.expand_dims(spec[:,1], axis=0)*np.expand_dims(filter_widths, axis=1), axis=1), axis=0)
    return flux_total
    
zarr = np.arange(0.01, 5.001, 0.01)
lbtarr = cosmo.lookback_time(zarr).value


### Load up the data necessary to build the "filter" to rebin the model spectra onto

hdulist1 = fits.open("../VANDELS_data/spectra/sc_206806_UDS_P1M1_MR_Q1_029_1.fits")
wavzpt = hdulist1[0].header["CRVAL1"]
dwav = hdulist1[0].header["CDELT1"]
fluxes1 = hdulist1[4].data#*10**19
print wavzpt
print dwav
raw_input()

maxwav = wavzpt + dwav*(len(fluxes1))
objwavs_nobin = np.arange(wavzpt, maxwav, dwav)
tauvals = np.array([0.05, 1, 10], dtype="str")