Commit c1356047 authored by Médéric Boquien's avatar Médéric Boquien
Browse files

Remove global variables. The arrays should really be passed as arguments.

parent 8707696f
......@@ -180,7 +180,7 @@ def save_table_best(obsid, chi2, chi2_red, variables, fluxes, filters,
format='ascii.commented_header')
def dchi2_over_ds2(s):
def dchi2_over_ds2(s, obs_fluxes, obs_errors, mod_fluxes):
"""Function used to estimate the normalization factor in the SED fitting
process when upper limits are included in the dataset to fit (from Eq. A11
in Sawicki M. 2012, PASA, 124, 1008).
......@@ -212,17 +212,17 @@ def dchi2_over_ds2(s):
# The mask "lim" selects the filter(s) for which upper limits are given
# i.e., when obs_fluxes is >=0. and obs_errors = 9990 <= obs_errors < 0.
wlim = np.where((gbl_obs_errors >= -9990.)&(gbl_obs_errors < 0.))
wdata = np.where(gbl_obs_errors>=0.)
wlim = np.where((obs_errors >= -9990.)&(obs_errors < 0.))
wdata = np.where(obs_errors>=0.)
mod_fluxes_data = gbl_mod_fluxes[wdata]
mod_fluxes_lim = gbl_mod_fluxes[wlim]
mod_fluxes_data = mod_fluxes[wdata]
mod_fluxes_lim = mod_fluxes[wlim]
obs_fluxes_data = gbl_obs_fluxes[wdata]
obs_fluxes_lim = gbl_obs_fluxes[wlim]
obs_fluxes_data = obs_fluxes[wdata]
obs_fluxes_lim = bs_fluxes[wlim]
obs_errors_data = gbl_obs_errors[wdata]
obs_errors_lim = -gbl_obs_errors[wlim]
obs_errors_data = obs_errors[wdata]
obs_errors_lim = -obs_errors[wlim]
dchi2_over_ds_data = np.sum(
(obs_fluxes_data-s*mod_fluxes_data) *
......
......@@ -207,8 +207,6 @@ def analysis(idx, obs):
# Tolerance threshold under which any flux or error is considered as 0.
tolerance = 1e-12
global gbl_mod_fluxes, gbl_obs_fluxes, gbl_obs_errors
# We pick the indices of the models with closest redshift assuming we have
# limited the number of decimals (usually set to 2 decimals).
w = np.where(gbl_w_redshifts[gbl_redshifts[np.abs(obs['redshift'] -
......@@ -240,9 +238,6 @@ def analysis(idx, obs):
lim_flag = gbl_lim_flag*np.logical_and(obs_errors >= -9990.,
obs_errors < tolerance)
gbl_obs_fluxes = obs_fluxes
gbl_obs_errors = obs_errors
# Normalisation factor to be applied to a model fluxes to best fit
# an observation fluxes. Normalised flux of the models. χ² and
# likelihood of the fitting. Reduced χ² (divided by the number of
......@@ -255,9 +250,10 @@ def analysis(idx, obs):
if lim_flag.any() is True:
norm_init = norm_facts
for imod in range(len(model_fluxes)):
gbl_mod_fluxes = model_fluxes[imod, :]
norm_facts[imod] = optimize.newton(dchi2_over_ds2, norm_init[imod],
tol=1e-16)
tol=1e-16,
args=(obs_fluxes, obs_errors,
model_fluxes[mod, :]))
model_fluxes *= norm_facts[:, np.newaxis]
# χ² of the comparison of each model to each observation.
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment