Python scipy.integrate.trapz() Examples
The following are 28
code examples of scipy.integrate.trapz().
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
scipy.integrate
, or try the search function
.
.
Example #1
| Source File: tools.py From burnman with GNU General Public License v2.0 | 6 votes |
|
def l2(x, funca, funcb):
"""
Computes the L2 norm for one profile(assumed to be linear between points).
:type x: array of float
:param x: depths :math:`[m]`.
:type funca: list of arrays of float
:param funca: array calculated values
:type funcb: list of arrays of float
:param funcb: array of values (observed or calculated) to compare to
:returns: L2 norm
:rtype: array of floats
"""
diff = np.array(funca - funcb)
diff = diff * diff
return integrate.trapz(diff, x)
Example #2
| Source File: thermo.py From pwtools with BSD 3-Clause "New" or "Revised" License | 6 votes |
|
def _integrate(self, y, f):
"""
Integrate `y` along axis=1, i.e. over freq axis for all T.
Parameters
----------
y : 2d array (nT, ndos) where nT = len(self.T), ndos = len(self.dos)
f : self.f, (len(self.dos),)
Returns
-------
array (nT,)
"""
mask = np.isnan(y)
if mask.any():
self._printwarn("HarmonicThermo._integrate: warning: "
" %i NaNs found in y!" %len(mask))
if self.fixnan:
self._printwarn("HarmonicThermo._integrate: warning: "
"fixing %i NaNs in y!" %len(mask))
y[mask] = self.nanfill
# this call signature works for scipy.integrate,{trapz,simps}
return self.integrator(y, x=f, axis=1)
Example #3
| Source File: thermo.py From pwtools with BSD 3-Clause "New" or "Revised" License | 6 votes |
|
def debye_func(x, nstep=100, zero=1e-8):
r"""Debye function
:math:`f(x) = 3 \int_0^1 t^3 / [\exp(t x) - 1] dt`
Parameters
----------
x : float or 1d array
nstep : int
number of points for integration
zero : float
approximate the 0 in the integral by this (small!) number
"""
x = np.atleast_1d(x)
if x.ndim == 1:
x = x[:,None]
else:
raise Exception("x is not 1d array")
tt = np.linspace(zero, 1.0, nstep)
return 3.0 * trapz(tt**3.0 / (np.exp(tt*x) - 1.0), tt, axis=1)
Example #4
| Source File: Inverse.py From PyMieScatt with MIT License | 6 votes |
|
def fastMie_SD(m, wavelength, dp, ndp):
# http://pymiescatt.readthedocs.io/en/latest/inverse.html#fastMie_SD
dp = coerceDType(dp)
ndp = coerceDType(ndp)
_length = np.size(dp)
Q_sca = np.zeros(_length)
Q_abs = np.zeros(_length)
Q_back = np.zeros(_length)
aSDn = np.pi*((dp/2)**2)*ndp*(1e-6)
for i in range(_length):
Q_sca[i],Q_abs[i],Q_back[i] = fastMieQ(m,wavelength,dp[i])
Bsca = trapz(Q_sca*aSDn,dp)
Babs = trapz(Q_abs*aSDn,dp)
Bback = trapz(Q_back*aSDn,dp)
return Bsca, Babs, Bback
Example #5
| Source File: Inverse.py From PyMieScatt with MIT License | 6 votes |
|
def fastMie_SD(m, wavelength, dp, ndp):
# http://pymiescatt.readthedocs.io/en/latest/inverse.html#fastMie_SD
dp = coerceDType(dp)
ndp = coerceDType(ndp)
_length = np.size(dp)
Q_sca = np.zeros(_length)
Q_abs = np.zeros(_length)
Q_back = np.zeros(_length)
aSDn = np.pi*((dp/2)**2)*ndp*(1e-6)
for i in range(_length):
Q_sca[i],Q_abs[i],Q_back[i] = fastMieQ(m,wavelength,dp[i])
Bsca = trapz(Q_sca*aSDn,dp)
Babs = trapz(Q_abs*aSDn,dp)
Bback = trapz(Q_back*aSDn,dp)
return Bsca, Babs, Bback
Example #6
| Source File: sandbox.py From pyphot with MIT License | 6 votes |
|
def leff(self):
""" Unitwise Effective wavelength
leff = int (lamb * T * Vega dlamb) / int(T * Vega dlamb)
"""
with Vega() as v:
s = self.reinterp(v.wavelength)
w = s._wavelength
if s.transmit.max() > 0:
leff = np.trapz(w * s.transmit * v.flux.value, w, axis=-1)
leff /= np.trapz(s.transmit * v.flux.value, w, axis=-1)
else:
leff = float('nan')
if s.wavelength_unit is not None:
leff = leff * Unit(s.wavelength_unit)
if self.wavelength_unit is not None:
return leff.to(self.wavelength_unit)
return leff
else:
return leff
Example #7
| Source File: integration.py From Conditional_Density_Estimation with MIT License | 6 votes |
|
def numeric_integation(func, n_samples=10 ** 5, bound_lower=-10**3, bound_upper=10**3):
""" Numeric integration over one dimension using the trapezoidal rule
Args:
func: function to integrate over - must take numpy arrays of shape (n_samples,) as first argument
and return a numpy array of shape (n_samples,)
n_samples: (int) number of samples
Returns:
approximated integral - numpy array of shape (ndim_out,)
"""
# proposal distribution
y_samples = np.squeeze(np.linspace(bound_lower, bound_upper, num=n_samples))
values = func(y_samples)
integral = integrate.trapz(values, y_samples)
return integral
Example #8
| Source File: sandbox.py From pyphot with MIT License | 6 votes |
|
def __init__(self, wavelength, transmit, name='', dtype="photon",
unit=None):
"""Constructor"""
self.name = name
self.set_dtype(dtype)
try: # get units from the inputs
self._wavelength = wavelength.value
unit = str(wavelength.unit)
except AttributeError:
self._wavelength = wavelength
self.set_wavelength_unit(unit)
self.transmit = np.clip(transmit, 0., np.nanmax(transmit))
self.norm = trapz(self.transmit, self._wavelength)
self._lT = trapz(self._wavelength * self.transmit, self._wavelength)
self._lpivot = self._calculate_lpivot()
if self.norm > 0:
self._cl = self._lT / self.norm
else:
self._cl = 0.
Example #9
| Source File: sandbox.py From pyphot with MIT License | 6 votes |
|
def leff(self):
""" Unitwise Effective wavelength
leff = int (lamb * T * Vega dlamb) / int(T * Vega dlamb)
"""
with Vega() as v:
s = self.reinterp(v.wavelength)
w = s._wavelength
if s.transmit.max() > 0:
leff = np.trapz(w * s.transmit * v.flux.value, w, axis=-1)
leff /= np.trapz(s.transmit * v.flux.value, w, axis=-1)
else:
leff = float('nan')
if s.wavelength_unit is not None:
leff = leff * Unit(s.wavelength_unit)
if self.wavelength_unit is not None:
return leff.to(self.wavelength_unit)
return leff
else:
return leff
Example #10
| Source File: sandbox.py From pyphot with MIT License | 6 votes |
|
def __init__(self, wavelength, transmit, name='', dtype="photon",
unit=None):
"""Constructor"""
self.name = name
self.set_dtype(dtype)
try: # get units from the inputs
self._wavelength = wavelength.value
unit = str(wavelength.unit)
except AttributeError:
self._wavelength = wavelength
self.set_wavelength_unit(unit)
self.transmit = np.clip(transmit, 0., np.nanmax(transmit))
self.norm = trapz(self.transmit, self._wavelength)
self._lT = trapz(self._wavelength * self.transmit, self._wavelength)
self._lpivot = self._calculate_lpivot()
if self.norm > 0:
self._cl = self._lT / self.norm
else:
self._cl = 0.
Example #11
| Source File: fc.py From CO2MPAS-TA with European Union Public License 1.1 | 5 votes |
|
def _rescaling_score(times, rescaling_matrix, k):
x = np.dot(rescaling_matrix, k)
m = np.trapz(x, times) / (times[-1] - times[0])
std = np.sqrt(np.trapz((x - m) ** 2, times) / (times[-1] - times[0]))
return m, std
Example #12
| Source File: slice.py From cgpm with Apache License 2.0 | 5 votes |
|
def main(num_samples, burn, lag, w):
alpha = 1.0
N = 25
Z = gu.simulate_crp(N, alpha)
K = max(Z) + 1
# CRP with gamma prior.
log_pdf_fun = lambda alpha : gu.logp_crp_unorm(N, K, alpha) - alpha
proposal_fun = lambda : np.random.gamma(1.0, 1.0)
D = (0, float('Inf'))
samples = su.slice_sample(proposal_fun, log_pdf_fun, D,
num_samples=num_samples, burn=burn, lag=lag, w=w)
minval = min(samples)
maxval = max(samples)
xvals = np.linspace(minval, maxval, 100)
yvals = np.array([math.exp(log_pdf_fun(x)) for x in xvals])
yvals /= trapz(xvals, yvals)
fig, ax = plt.subplots(2,1)
ax[0].hist(samples, 50, normed=True)
ax[1].hist(samples, 100, normed=True)
ax[1].plot(xvals,-yvals, c='red', lw=3, alpha=.8)
ax[1].set_xlim(ax[0].get_xlim())
ax[1].set_ylim(ax[0].get_ylim())
plt.show()
Example #13
| Source File: co2.py From CO2MPAS-TA with European Union Public License 1.1 | 5 votes |
|
def calculate_phases_co2_emissions(
times, phases_indices, co2_emissions, phases_distances=1.0):
"""
Calculates CO2 emission or cumulative CO2 of cycle phases [CO2g/km or CO2g].
If phases_distances is not given the result is the cumulative CO2 of cycle
phases [CO2g] otherwise it is CO2 emission of cycle phases [CO2g/km].
:param times:
Time vector [s].
:type times: numpy.array
:param phases_indices:
Indices of the cycle phases [-].
:type phases_indices: numpy.array
:param co2_emissions:
CO2 instantaneous emissions vector [CO2g/s].
:type co2_emissions: numpy.array
:param phases_distances:
Cycle phases distances [km].
:type phases_distances: numpy.array | float, optional
:return:
CO2 emission or cumulative CO2 of cycle phases [CO2g/km or CO2g].
:rtype: numpy.array
"""
from scipy.integrate import trapz
co2 = []
for i, j in phases_indices:
co2.append(trapz(co2_emissions[i:j], times[i:j]))
with np.errstate(divide='ignore', invalid='ignore'):
return np.nan_to_num(np.array(co2) / phases_distances)
Example #14
| Source File: phot.py From pyphot with MIT License | 5 votes |
|
def _calculate_lpivot(self):
if 'photon' in self.dtype:
lpivot2 = self._lT / trapz(self.transmit / self._wavelength, self._wavelength)
else:
lpivot2 = self.norm / trapz(self.transmit / self._wavelength ** 2, self._wavelength)
return np.sqrt(lpivot2)
Example #15
| Source File: phot.py From pyphot with MIT License | 5 votes |
|
def __init__(self, wavelength, transmit, name='', dtype="photon", unit=None):
"""Constructor"""
self.name = name
self.set_dtype(dtype)
try: # get units from the inputs
self._wavelength = wavelength.magnitude
unit = str(wavelength.units)
except AttributeError:
self._wavelength = wavelength
self.set_wavelength_unit(unit)
self.transmit = np.clip(transmit, 0., np.nanmax(transmit))
self.norm = trapz(self.transmit, self._wavelength)
self._lT = trapz(self._wavelength * self.transmit, self._wavelength)
self._lpivot = self._calculate_lpivot()
self._cl = self._lT / self.norm
Example #16
| Source File: sandbox.py From pyphot with MIT License | 5 votes |
|
def _get_indice(cls, w, flux, blue, red, band=None, unit='ew', degree=1,
**kwargs):
""" compute spectral index after continuum subtraction
Parameters
----------
w: ndarray (nw, )
array of wavelengths in AA
flux: ndarray (N, nw)
array of flux values for different spectra in the series
blue: tuple(2)
selection for blue continuum estimate
red: tuple(2)
selection for red continuum estimate
band: tuple(2), optional
select region in this band only.
default is band = (min(blue), max(red))
unit: str
`ew` or `mag` wether equivalent width or magnitude
degree: int (default 1)
degree of the polynomial fit to the continuum
Returns
-------
ew: ndarray (N,)
equivalent width array
"""
wi, fi = cls.continuum_normalized_region_around_line(w, flux, blue,
red, band=band,
degree=degree)
if unit in (0, 'ew', 'EW'):
return np.trapz(1. - fi, wi, axis=-1)
else:
m = np.trapz(fi, wi, axis=-1)
m = -2.5 * np.log10(m / np.ptp(wi))
return m
Example #17
| Source File: sandbox.py From pyphot with MIT License | 5 votes |
|
def _calculate_lpivot(self):
if self.transmit.max() <= 0:
return 0.
if 'photon' in self.dtype:
lpivot2 = self._lT / trapz(self.transmit / self._wavelength,
self._wavelength)
else:
lpivot2 = self.norm / trapz(self.transmit / self._wavelength ** 2,
self._wavelength)
return np.sqrt(lpivot2)
Example #18
| Source File: Mie.py From PyMieScatt with MIT License | 5 votes |
|
def SF_SD(m, wavelength, dp, ndp, nMedium=1.0, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#SF_SD
nMedium = nMedium.real
m /= nMedium
wavelength /= nMedium
_steps = int(1+(maxAngle-minAngle)/angularResolution)
ndp = coerceDType(ndp)
dp = coerceDType(dp)
SL = np.zeros(_steps)
SR = np.zeros(_steps)
SU = np.zeros(_steps)
kwargs = {'minAngle':minAngle,
'maxAngle':maxAngle,
'angularResolution':angularResolution,
'space':space,
'normalization':None}
for n,d in zip(ndp,dp):
measure,l,r,u = ScatteringFunction(m,wavelength,d,**kwargs)
SL += l*n
SR += r*n
SU += u*n
if normalization in ['n','N','number','particles']:
_n = trapz(ndp,dp)
SL /= _n
SR /= _n
SU /= _n
elif normalization in ['m','M','max','MAX']:
SL /= np.max(SL)
SR /= np.max(SR)
SU /= np.max(SU)
elif normalization in ['t','T','total','TOTAL']:
SL /= trapz(SL,measure)
SR /= trapz(SR,measure)
SU /= trapz(SU,measure)
return measure,SL,SR,SU
Example #19
| Source File: Mie.py From PyMieScatt with MIT License | 5 votes |
|
def SF_SD(m, wavelength, dp, ndp, nMedium=1.0, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#SF_SD
nMedium = nMedium.real
m /= nMedium
wavelength /= nMedium
_steps = int(1+(maxAngle-minAngle)/angularResolution)
ndp = coerceDType(ndp)
dp = coerceDType(dp)
SL = np.zeros(_steps)
SR = np.zeros(_steps)
SU = np.zeros(_steps)
kwargs = {'minAngle':minAngle,
'maxAngle':maxAngle,
'angularResolution':angularResolution,
'space':space,
'normalization':None}
for n,d in zip(ndp,dp):
measure,l,r,u = ScatteringFunction(m,wavelength,d,**kwargs)
SL += l*n
SR += r*n
SU += u*n
if normalization in ['n','N','number','particles']:
_n = trapz(ndp,dp)
SL /= _n
SR /= _n
SU /= _n
elif normalization in ['m','M','max','MAX']:
SL /= np.max(SL)
SR /= np.max(SR)
SU /= np.max(SU)
elif normalization in ['t','T','total','TOTAL']:
SL /= trapz(SL,measure)
SR /= trapz(SR,measure)
SU /= trapz(SU,measure)
return measure,SL,SR,SU
Example #20
| Source File: Mie.py From PyMieScatt with MIT License | 5 votes |
|
def Mie_SD(m, wavelength, dp, ndp, nMedium=1.0, interpolate=False, asDict=False):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_SD
nMedium = nMedium.real
m /= nMedium
wavelength /= nMedium
dp = coerceDType(dp)
ndp = coerceDType(ndp)
_length = np.size(dp)
Q_ext = np.zeros(_length)
Q_sca = np.zeros(_length)
Q_abs = np.zeros(_length)
Q_pr = np.zeros(_length)
Q_back = np.zeros(_length)
Q_ratio = np.zeros(_length)
g = np.zeros(_length)
# scaling of 1e-6 to cast in units of inverse megameters - see docs
aSDn = np.pi*((dp/2)**2)*ndp*(1e-6)
# _logdp = np.log10(dp)
for i in range(_length):
Q_ext[i], Q_sca[i], Q_abs[i], g[i], Q_pr[i], Q_back[i], Q_ratio[i] = AutoMieQ(m,wavelength,dp[i],nMedium)
Bext = trapz(Q_ext*aSDn,dp)
Bsca = trapz(Q_sca*aSDn,dp)
Babs = Bext-Bsca
Bback = trapz(Q_back*aSDn,dp)
Bratio = trapz(Q_ratio*aSDn,dp)
bigG = trapz(g*Q_sca*aSDn,dp)/trapz(Q_sca*aSDn,dp)
Bpr = Bext - bigG*Bsca
if asDict:
return dict(Bext=Bext, Bsca=Bsca, Babs=Babs, G=bigG, Bpr=Bpr, Bback=Bback, Bratio=Bratio)
else:
return Bext, Bsca, Babs, bigG, Bpr, Bback, Bratio
Example #21
| Source File: num.py From pwtools with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
def norm_int(y, x, area=1.0, scale=True, func=simps):
"""Normalize integral area of y(x) to `area`.
Parameters
----------
x,y : numpy 1d arrays
area : float
scale : bool, optional
Scale x and y to the same order of magnitude before integration.
This may be necessary to avoid numerical trouble if x and y have very
different scales.
func : callable
Function to do integration (like scipy.integrate.{simps,trapz,...}
Called as ``func(y,x)``. Default: simps
Returns
-------
scaled y
Notes
-----
The argument order y,x might be confusing. x,y would be more natural but we
stick to the order used in the scipy.integrate routines.
"""
if scale:
fx = np.abs(x).max()
fy = np.abs(y).max()
sx = x / fx
sy = y / fy
else:
fx = fy = 1.0
sx, sy = x, y
# Area under unscaled y(x).
_area = func(sy, sx) * fx * fy
return y*area/_area
Example #22
| Source File: MODEL_AGNfitter.py From AGNfitter with MIT License | 4 votes |
|
def filters1( model_nus, model_fluxes, filterdict, z ):
"""
Projects the model SEDs into the filter curves of each photometric band.
##input:
- model_nus: template frequencies [log10(nu)]
- model_fluxes: template fluxes [F_nu]
- filterdict: dictionary with all band filter curves' information.
To change this, add one band and filter curve, etc,
look at DICTIONARIES_AGNfitter.py
- z: redshift
##output:
- bands [log10(nu)]
- Filtered fluxes at these bands [F_nu]
"""
bands, files_dict, lambdas_dict, factors_dict = filterdict
filtered_model_Fnus = []
# Costumize model frequencies and fluxes [F_nu]
# to same units as filter curves (to wavelengths [angstrom] and F_lambda)
model_lambdas = nu2lambda_angstrom(model_nus) * (1+z)
model_lambdas = model_lambdas[::-1]
model_fluxes_nu = model_fluxes[::-1]
model_fluxes_lambda = fluxnu_2_fluxlambda(model_fluxes_nu, model_lambdas)
mod2filter_interpol = interp1d(model_lambdas, model_fluxes_lambda, kind = 'nearest', bounds_error=False, fill_value=0.)
# For filter curve at each band.
# (Vectorised integration was not possible -> different filter-curve-arrays' sizes)
for iband in bands:
# Read filter curves info for each data point
# (wavelengths [angstrom] and factors [non])
lambdas_filter = np.array(lambdas_dict[iband])
factors_filter = np.array(factors_dict[iband])
iband_angst = nu2lambda_angstrom(iband)
# Interpolate the model fluxes to
#the exact wavelengths of filter curves
modelfluxes_at_filterlambdas = mod2filter_interpol(lambdas_filter)
# Compute the flux ratios, equivalent to the filtered fluxes:
# F = int(model)/int(filter)
integral_model = trapz(modelfluxes_at_filterlambdas*factors_filter, x= lambdas_filter)
integral_filter = trapz(factors_filter, x= lambdas_filter)
filtered_modelF_lambda = (integral_model/integral_filter)
# Convert all from lambda, F_lambda to Fnu and nu
filtered_modelFnu_atfilter_i = fluxlambda_2_fluxnu(filtered_modelF_lambda, iband_angst)
filtered_model_Fnus.append(filtered_modelFnu_atfilter_i)
return bands, np.array(filtered_model_Fnus)
Example #23
| Source File: sandbox.py From pyphot with MIT License | 4 votes |
|
def get_flux(self, slamb, sflux, axis=-1):
"""getFlux
Integrate the flux within the filter and return the integrated energy
If you consider applying the filter to many spectra, you might want to
consider extractSEDs.
Parameters
----------
slamb: ndarray(dtype=float, ndim=1)
spectrum wavelength definition domain
sflux: ndarray(dtype=float, ndim=1)
associated flux
Returns
-------
flux: float
Energy of the spectrum within the filter
"""
passb = self.reinterp(slamb)
ifT = passb.transmit
_slamb = _drop_units(slamb)
_sflux = _drop_units(passb._validate_sflux(slamb, sflux))
_w_unit = str(slamb.unit)
_f_unit = str(sflux.unit)
# if the filter is null on that wavelength range flux is then 0
# ind = ifT > 0.
nonzero = np.where(ifT > 0)[0]
if nonzero.size <= 0:
return passb._get_zero_like(sflux)
# avoid calculating many zeros
nonzero_start = max(0, min(nonzero) - 5)
nonzero_end = min(len(ifT), max(nonzero) + 5)
ind = np.zeros(len(ifT), dtype=bool)
ind[nonzero_start:nonzero_end] = True
if True in ind:
try:
_sflux = _sflux[:, ind]
except Exception:
_sflux = _sflux[ind]
# limit integrals to where necessary
if 'photon' in passb.dtype:
a = np.trapz(_slamb[ind] * ifT[ind] * _sflux, _slamb[ind],
axis=axis)
b = np.trapz(_slamb[ind] * ifT[ind], _slamb[ind])
a = a * Unit('*'.join((_w_unit, _f_unit, _w_unit)))
b = b * Unit('*'.join((_w_unit, _w_unit)))
elif 'energy' in passb.dtype:
a = np.trapz(ifT[ind] * _sflux, _slamb[ind], axis=axis)
b = np.trapz(ifT[ind], _slamb[ind])
a = a * Unit('*'.join((_f_unit, _w_unit)))
b = b * Unit(_w_unit)
if (np.isinf(a.value).any() | np.isinf(b.value).any()):
print(self.name, "Warn for inf value")
return a / b
else:
return passb._get_zero_like(sflux)
Example #24
| Source File: iterative_wiener.py From pyroomacoustics with MIT License | 4 votes |
|
def compute_squared_gain(a, noise_psd, y):
"""
Estimate the squared gain of the speech power spectral density as done on
p. 204 of
J. Lim and A. Oppenheim, *All-Pole Modeling of Degraded Speech,*
IEEE Transactions on Acoustics, Speech, and Signal Processing 26.3
(1978): 197-210.
Namely solving for :math:`g^2` such that the following expression is
satisfied:
.. math::
\dfrac{N}{2\pi} \int_{-\pi}^{\pi} \dfrac{g^2}{\\left \| 1 - \sum_{k=1}^p a_k \cdot e^{-jk\omega} \\right \|^2} d\omega = \sum_{n=0}^{N-1} y^2(n) - N\cdot\sigma_d^2,
where :math:`N` is the number of noisy samples :math:`y`, :math:`a_k`
are the :math:`p` LPC coefficients, and :math:`\sigma_d^2` is the
noise variance.
Parameters
----------
a : numpy array
LPC coefficients.
noise_psd : float or numpy array
Noise variance if white noise, numpy array otherwise.
y : numpy array
Noisy time domain samples.
Returns
-------
float
Squared gain.
"""
p = len(a)
def _lpc_all_pole(omega):
k = np.arange(p) + 1
return 1 / np.abs(1 - np.dot(a, np.exp(-1j * k * omega))) ** 2
N = len(y)
# right hand side of expression
if np.isscalar(noise_psd): # white noise, i.e. flat spectrum
rhs = np.sum(y**2) - N * noise_psd
else:
rhs = np.sum(y**2) - np.sum(noise_psd)
# estimate integral
d_omega = 2 * np.pi / 1000
omega_vals = np.arange(-np.pi, np.pi, d_omega)
vec_integrand = np.vectorize(_lpc_all_pole)
integral = integrate.trapz(vec_integrand(omega_vals), omega_vals)
return rhs * 2 * np.pi / N / integral
Example #25
| Source File: sandbox.py From pyphot with MIT License | 4 votes |
|
def get_flux(self, slamb, sflux, axis=-1):
"""getFlux
Integrate the flux within the filter and return the integrated energy
If you consider applying the filter to many spectra, you might want to
consider extractSEDs.
Parameters
----------
slamb: ndarray(dtype=float, ndim=1)
spectrum wavelength definition domain
sflux: ndarray(dtype=float, ndim=1)
associated flux
Returns
-------
flux: float
Energy of the spectrum within the filter
"""
passb = self.reinterp(slamb)
ifT = passb.transmit
_slamb = _drop_units(slamb)
_sflux = _drop_units(passb._validate_sflux(slamb, sflux))
_w_unit = str(slamb.unit)
_f_unit = str(sflux.unit)
# if the filter is null on that wavelength range flux is then 0
# ind = ifT > 0.
nonzero = np.where(ifT > 0)[0]
if nonzero.size <= 0:
return passb._get_zero_like(sflux)
# avoid calculating many zeros
nonzero_start = max(0, min(nonzero) - 5)
nonzero_end = min(len(ifT), max(nonzero) + 5)
ind = np.zeros(len(ifT), dtype=bool)
ind[nonzero_start:nonzero_end] = True
if True in ind:
try:
_sflux = _sflux[:, ind]
except Exception:
_sflux = _sflux[ind]
# limit integrals to where necessary
if 'photon' in passb.dtype:
a = np.trapz(_slamb[ind] * ifT[ind] * _sflux, _slamb[ind],
axis=axis)
b = np.trapz(_slamb[ind] * ifT[ind], _slamb[ind])
a = a * Unit('*'.join((_w_unit, _f_unit, _w_unit)))
b = b * Unit('*'.join((_w_unit, _w_unit)))
elif 'energy' in passb.dtype:
a = np.trapz(ifT[ind] * _sflux, _slamb[ind], axis=axis)
b = np.trapz(ifT[ind], _slamb[ind])
a = a * Unit('*'.join((_f_unit, _w_unit)))
b = b * Unit(_w_unit)
if (np.isinf(a.value).any() | np.isinf(b.value).any()):
print(self.name, "Warn for inf value")
return a / b
else:
return passb._get_zero_like(_sflux)
Example #26
| Source File: Mie.py From PyMieScatt with MIT License | 4 votes |
|
def ScatteringFunction(m, wavelength, diameter, nMedium=1.0, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#ScatteringFunction
nMedium = nMedium.real
m /= nMedium
wavelength /= nMedium
x = np.pi*diameter/wavelength
_steps = int(1+(maxAngle-minAngle)/angularResolution) # default 361
if angleMeasure in ['radians','RADIANS','rad','RAD']:
adjust = np.pi/180
elif angleMeasure in ['gradians','GRADIANS','grad','GRAD']:
adjust = 1/200
else:
adjust = 1
if space in ['q','qspace','QSPACE','qSpace']:
# _steps *= 10
_steps += 1
if minAngle==0:
minAngle = 1e-5
#measure = np.logspace(np.log10(minAngle),np.log10(maxAngle),_steps)*np.pi/180
measure = np.linspace(minAngle, maxAngle, _steps)*np.pi/180
_q = True
else:
measure = np.linspace(minAngle,maxAngle,_steps)*adjust
_q = False
if x == 0:
return measure,0,0,0
_measure = np.linspace(minAngle,maxAngle,_steps)*np.pi/180
SL = np.zeros(_steps)
SR = np.zeros(_steps)
SU = np.zeros(_steps)
for j in range(_steps):
u = np.cos(_measure[j])
S1, S2 = MieS1S2(m,x,u)
SL[j] = (np.sum(np.conjugate(S1)*S1)).real
SR[j] = (np.sum(np.conjugate(S2)*S2)).real
SU[j] = (SR[j]+SL[j])/2
if normalization in ['m','M','max','MAX']:
SL /= np.max(SL)
SR /= np.max(SR)
SU /= np.max(SU)
elif normalization in ['t','T','total','TOTAL']:
SL /= trapz(SL,measure)
SR /= trapz(SR,measure)
SU /= trapz(SU,measure)
if _q:
measure = (4*np.pi/wavelength)*np.sin(measure/2)*(diameter/2)
return measure,SL,SR,SU
Example #27
| Source File: Mie.py From PyMieScatt with MIT License | 4 votes |
|
def ScatteringFunction(m, wavelength, diameter, nMedium=1.0, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#ScatteringFunction
nMedium = nMedium.real
m /= nMedium
wavelength /= nMedium
x = np.pi*diameter/wavelength
_steps = int(1+(maxAngle-minAngle)/angularResolution) # default 361
if angleMeasure in ['radians','RADIANS','rad','RAD']:
adjust = np.pi/180
elif angleMeasure in ['gradians','GRADIANS','grad','GRAD']:
adjust = 1/200
else:
adjust = 1
if space in ['q','qspace','QSPACE','qSpace']:
# _steps *= 10
_steps += 1
if minAngle==0:
minAngle = 1e-5
#measure = np.logspace(np.log10(minAngle),np.log10(maxAngle),_steps)*np.pi/180
measure = np.linspace(minAngle, maxAngle, _steps)*np.pi/180
_q = True
else:
measure = np.linspace(minAngle,maxAngle,_steps)*adjust
_q = False
if x == 0:
return measure,0,0,0
_measure = np.linspace(minAngle,maxAngle,_steps)*np.pi/180
SL = np.zeros(_steps)
SR = np.zeros(_steps)
SU = np.zeros(_steps)
for j in range(_steps):
u = np.cos(_measure[j])
S1, S2 = MieS1S2(m,x,u)
SL[j] = (np.sum(np.conjugate(S1)*S1)).real
SR[j] = (np.sum(np.conjugate(S2)*S2)).real
SU[j] = (SR[j]+SL[j])/2
if normalization in ['m','M','max','MAX']:
SL /= np.max(SL)
SR /= np.max(SR)
SU /= np.max(SU)
elif normalization in ['t','T','total','TOTAL']:
SL /= trapz(SL,measure)
SR /= trapz(SR,measure)
SU /= trapz(SU,measure)
if _q:
measure = (4*np.pi/wavelength)*np.sin(measure/2)*(diameter/2)
return measure,SL,SR,SU
Example #28
| Source File: phot.py From pyphot with MIT License | 4 votes |
|
def getFlux(self, slamb, sflux, axis=-1):
"""getFlux
Integrate the flux within the filter and return the integrated energy
If you consider applying the filter to many spectra, you might want to
consider extractSEDs.
Parameters
----------
slamb: ndarray(dtype=float, ndim=1)
spectrum wavelength definition domain
sflux: ndarray(dtype=float, ndim=1)
associated flux
Returns
-------
flux: float
Energy of the spectrum within the filter
"""
_wavelength = self._get_filter_in_units_of(slamb)
_slamb = _drop_units(slamb)
#clean data for inf values by interpolation
if True in np.isinf(sflux):
indinf = np.where(np.isinf(sflux))
indfin = np.where(np.isfinite(sflux))
sflux[indinf] = np.interp(_slamb[indinf], _slamb[indfin], sflux[indfin])
# reinterpolate transmission onto the same wavelength def as the data
ifT = np.interp(_slamb, _wavelength, self.transmit, left=0., right=0.)
# if the filter is null on that wavelength range flux is then 0
# ind = ifT > 0.
nonzero = np.where(ifT > 0)[0]
if nonzero.size <= 0:
return 0.
nonzero_start = max(0, min(nonzero) - 5)
nonzero_end = min(len(ifT), max(nonzero) + 5)
ind = np.zeros(len(ifT), dtype=bool)
ind[nonzero_start:nonzero_end] = True
if True in ind:
try:
_sflux = sflux[:, ind]
except:
_sflux = sflux[ind]
# limit integrals to where necessary
# ind = ifT > 0.
if 'photon' in self.dtype:
a = np.trapz(_slamb[ind] * ifT[ind] * _sflux, _slamb[ind], axis=axis)
b = np.trapz(_slamb[ind] * ifT[ind], _slamb[ind] )
elif 'energy' in self.dtype:
a = np.trapz( ifT[ind] * _sflux, _slamb[ind], axis=axis )
b = np.trapz( ifT[ind], _slamb[ind])
if (np.isinf(a).any() | np.isinf(b).any()):
print(self.name, "Warn for inf value")
return a / b
else:
return 0.
