Skip to content

Commit

Permalink
Cleaned up version of the extinction fix.
Browse files Browse the repository at this point in the history
  • Loading branch information
aprsa committed Sep 20, 2024
1 parent 1a124f4 commit 5d0c182
Showing 1 changed file with 1 addition and 32 deletions.
33 changes: 1 addition & 32 deletions phoebe/atmospheres/passbands.py
Original file line number Diff line number Diff line change
Expand Up @@ -889,14 +889,6 @@ def compute_intensities(self, atm, path, include_extinction=False, rvs=None, ebv
ext_photon_grid = np.empty(shape=(len(np.unique(teffs)), len(np.unique(loggs)), len(np.unique(abuns)), len(ebvs), len(rvs), 1))
ext_energy_grid = np.empty(shape=(len(np.unique(teffs)), len(np.unique(loggs)), len(np.unique(abuns)), len(ebvs), len(rvs), 1))

# ebv_list = np.tile(np.repeat(ebvs, len(rvs)), nmodels)
# rv_list = np.tile(rvs, nmodels*len(ebvs))

# ext_matrix = np.rollaxis(np.array([np.split(rv_list*ebv_list, nmodels), np.split(ebv_list, nmodels)]), 1)
# ext_matrix = np.ascontiguousarray(ext_matrix)

# ext_energy, ext_photon = np.empty((nmodels, len(rvs)*len(ebvs))), np.empty((nmodels, len(rvs)*len(ebvs)))

keep = (wls >= self.ptf_table['wl'][0]) & (wls <= self.ptf_table['wl'][-1])
wls = wls[keep]

Expand Down Expand Up @@ -924,16 +916,11 @@ def compute_intensities(self, atm, path, include_extinction=False, rvs=None, ebv
ax, bx = axbx[:,0], axbx[:,1]
for ei, ebv in enumerate(ebvs):
for ri, rv in enumerate(rvs):
Alam=10**(-0.4 * ebv * (rv * ax + bx))
Alam = 10**(-0.4 * ebv * (rv * ax + bx))
t = (teffs[i] == ext_axes[0], loggs[i] == ext_axes[1], abuns[i] == ext_axes[2], ebvs[ei] == ext_axes[3], rvs[ri] == ext_axes[4],0)
ext_energy_grid[t] = np.trapz(self.ptf(wls) * seds * Alam, wls)[-1] / np.trapz(self.ptf(wls) * seds, wls)[-1]
ext_photon_grid[t] = np.trapz(wls * self.ptf(wls) * seds * Alam, wls)[-1] / np.trapz(wls * self.ptf(wls) * seds, wls)[-1]


# ext_lambda = np.matmul(libphoebe.gordon_extinction(wls), ext_matrix[i])
# flux_frac = np.exp(-0.9210340371976184*ext_lambda) #10**(-0.4*ext_lambda)
# ext_energy[i], ext_photon[i] = np.dot([pbints_energy[-1]/fluxes_energy[-1], pbints_photon[-1]/fluxes_photon[-1]], flux_frac)

basic_axes = (np.unique(teffs), np.unique(loggs), np.unique(abuns))
self.ndp[atm] = ndpolator.Ndpolator(basic_axes=basic_axes)

Expand Down Expand Up @@ -963,24 +950,6 @@ def compute_intensities(self, atm, path, include_extinction=False, rvs=None, ebv
self.ndp[atm].register('ext@photon', associated_axes, ext_photon_grid)
self.ndp[atm].register('ext@energy', associated_axes, ext_energy_grid)

# teffs = np.repeat(teffs, len(ebvs)*len(rvs))
# loggs = np.repeat(loggs, len(ebvs)*len(rvs))
# abuns = np.repeat(abuns, len(ebvs)*len(rvs))

# ext_energy_grid = np.full(shape=[len(axis) for axis in axes]+[1], fill_value=np.nan)
# ext_photon_grid = np.copy(ext_energy_grid)

# flat_energy = ext_energy.flat
# flat_photon = ext_photon.flat

# for i in range(nmodels*len(ebvs)*len(rvs)):
# t = (teffs[i] == axes[0], loggs[i] == axes[1], abuns[i] == axes[2], ebv_list[i] == axes[3], rv_list[i] == axes[4])
# ext_energy_grid[t] = flat_energy[i]
# ext_photon_grid[t] = flat_photon[i]

# self.ndp[atm].register('ext@photon', associated_axes, ext_photon_grid)
# self.ndp[atm].register('ext@energy', associated_axes, ext_energy_grid)

if f'{atm}:ext' not in self.content:
self.content.append(f'{atm}:ext')

Expand Down

0 comments on commit 5d0c182

Please sign in to comment.