Python Class.get_perturbations方法代码示例

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import get_perturbations [as 别名]

#.........这里部分代码省略.........

        if texname:
            self.set_texnames({varied_name: texname})
        elif key not in self.texnames:  # texname will not be set at this stage. No check required
            self.set_texnames({varied_name: varied_name})

        if (not update) or (key not in self.computed.keys()):
            self.computed[key] = od()

        for val in values:
            # key = "{}={}".format(varied_name, val)
            params["parameters_smg"] = inip.vary_params(params["parameters_smg"], [[index_variable, val]])

            # It might be after the try to not store empty dictionaries.
            # Nevertheless, I find more useful having them to keep track of
            # those failed and, perhaps, to implement a method to obtain them
            # with Omega_smg_debug.
            d = self.computed[key][val] = {}

            self.cosmo.empty()
            self.cosmo.set(params)

            try:
                self.cosmo.compute()
            except Exception, e:
                print "Error: skipping {}={}".format(key, val)
                if cosmo_msg:
                    print e

                continue

            d['tunned'] = self.cosmo.get_current_derived_parameters(['tuning_parameter'])['tuning_parameter']

            for lst in [[back, 'back', self.cosmo.get_background],
                        [thermo, 'thermo', self.cosmo.get_thermodynamics],
                        [prim, 'prim', self.cosmo.get_thermodynamics]]:
                if lst[0]:
                    output = lst[2]()
                    if lst[0][0] == 'all':
                        d[lst[1]] = output
                    else:
                        d[lst[1]] = {}
                        for item in back:
                            if type(item) is list:
                                d[lst[1]].update({item[0]: output[item[0]][item[1]]})
                            else:
                                d[lst[1]].update({item: output[item]})

            if pert:
                # Perturbation is tricky because it can accept two optional
                # argument for get_perturbations and this method returns a
                # dictionary {'kind_of_pert': [{variable: list_values}]}, where
                # each item in the list is for a k (chosen in params).
                if type(pert[0]) is dict:
                    output = self.cosmo.get_perturbations(pert[0]['z'], pert[0]['output_format'])
                    if pert[1] == 'all':
                        d['pert'] = output
                else:
                    output = self.cosmo.get_perturbations()
                    if pert[0] == 'all':
                        d['pert'] = output

                if (type(pert[0]) is not dict) and (pert[0] != 'all'):
                    d['pert'] = {}
                    for subkey, lst in output.iteritems():
                        d['pert'].update({subkey: []})
                        for n, kdic in enumerate(lst):  # Each item is for a k
                            d['pert'][subkey].append({})
                            for item in pert:
                                if type(item) is list:
                                    d['pert'][subkey][n].update({item[0]: kdic[item[0]][item[1]]})
                                else:
                                    d['pert'][subkey][n].update({item: kdic[item]})

            for i in extra:
                exec('d[i] = self.cosmo.{}'.format(i))

            try:
                d['cl'] = self.__store_cl(self.cosmo.raw_cl())
            except CosmoSevereError:
                pass

            try:
                d['lcl'] = self.__store_cl(self.cosmo.lensed_cl())
            except CosmoSevereError:
                pass

            try:
                d['dcl'] = self.cosmo.density_cl()
            except CosmoSevereError:
                pass


            if ("output" in self.cosmo.pars) and ('mPk' in self.cosmo.pars['output']):
                k_array = np.linspace(*pk)
                pk_array = np.array([self.cosmo.pk(k, 0) for k in k_array])

                d['pk'] = {'k': k_array, 'pk': pk_array}

            self.cosmo.struct_cleanup()

# 需要导入模块: from classy import Class [as 别名]
# 或者: from classy.Class import get_perturbations [as 别名]
                   # Take fixed value for primordial Helium (instead of automatic BBN adjustment)
                   'YHe':0.246,
                   # other options and settings
                   'compute damping scale':'yes', # needed to output the time of damping scale crossing
                   'gauge':'newtonian'}
##############
#
# call CLASS
#
M = Class()
M.set(common_settings)
M.compute()
#
# load perturbations
#
all_k = M.get_perturbations()  # this potentially constains scalars/tensors and all k values
print all_k['scalar'][0].viewkeys()
#
one_k = all_k['scalar'][0]     # this contains only the scalar perturbations for the requested k values
#
tau = one_k['tau [Mpc]']
Theta0 = 0.25*one_k['delta_g']
phi = one_k['phi']
psi = one_k['psi']
theta_b = one_k['theta_b']
a = one_k['a']
# compute related quantitites
R = 3./4.*M.Omega_b()/M.Omega_g()*a    # R = 3/4 * (rho_b/rho_gamma)
zero_point = -(1.+R)*psi               # zero point of oscillations: -(1.+R)*psi
#
# get Theta0 oscillation amplitude (for vertical scale of plot)