def test_norm_int():
# simps(y, x) == 2.0
x = np.linspace(0, np.pi, 100)
y = np.sin(x)
for scale in [True, False]:
yy = norm_int(y, x, area=10.0, scale=scale)
assert np.allclose(simps(yy,x), 10.0)
python类simps()的实例源码
def AUCError(errors, failureThreshold, step=0.0001, showCurve=False):
nErrors = len(errors)
xAxis = list(np.arange(0., failureThreshold + step, step))
ced = [float(np.count_nonzero([errors <= x])) / nErrors for x in xAxis]
AUC = simps(ced, x=xAxis) / failureThreshold
failureRate = 1. - ced[-1]
print "AUC @ {0}: {1}".format(failureThreshold, AUC)
print "Failure rate: {0}".format(failureRate)
if showCurve:
plt.plot(xAxis, ced)
plt.show()
def emissivity(self, frequencies):
"""
Calculate the synchrotron emissivity (power emitted per volume
and per frequency) at the requested frequency.
NOTE
----
Since ``self.gamma`` and ``self.n_e`` are sampled on a logarithmic
grid, we integrate over ``ln(gamma)`` instead of ``gamma`` directly:
I = int_gmin^gmax f(g) d(g)
= int_ln(gmin)^ln(gmax) f(g) g d(ln(g))
XXX
---
Assume that the electrons have a pitch angle of ``pi/2`` with
respect to the magnetic field. (I think it is a good simplification
considering that the magnetic field is also assumed to be uniform.)
Parameters
----------
frequencies : float, or 1D `~numpy.ndarray`
The frequencies where to calculate the synchrotron emissivity.
Unit: [MHz]
Returns
-------
syncem : float, or 1D `~numpy.ndarray`
The calculated synchrotron emissivity at each specified
frequency.
Unit: [erg/s/cm^3/Hz]
"""
j_coef = np.sqrt(3) * AC.e**3 * self.B_gauss / AU.mec2
nu_c = self.frequency_crit(self.gamma)
frequencies = np.array(frequencies, ndmin=1)
syncem = np.zeros(shape=frequencies.shape)
for i, freq in enumerate(frequencies):
logger.debug("Calculating emissivity at %.2f [MHz]" % freq)
kernel = self.F(freq / nu_c)
# Integrate over energy ``gamma`` in logarithmic grid
syncem[i] = j_coef * integrate.simps(
self.n_e*kernel*self.gamma, x=np.log(self.gamma))
if len(syncem) == 1:
return syncem[0]
else:
return syncem
def KM_Cond_Average(x_min,x_max,Ref,G,Dh,Tbubble,Tdew,p,beta,satTransport=None):
"""
Returns the average pressure gradient and average heat transfer coefficient
between qualities of x_min and x_max.
for Kim&Mudawar two-phase condensation in mico-channel HX
To obtain the pressure gradient for a given value of x, pass it in as x_min and x_max
Required parameters:
* x_min : The minimum quality for the range [-]
* x_max : The maximum quality for the range [-]
* Ref : String with the refrigerant name
* G : Mass flux [kg/m^2/s]
* Dh : Hydraulic diameter of tube [m]
* Tbubble : Bubblepoint temperature of refrigerant [K]
* Tdew : Dewpoint temperature of refrigerant [K]
* beta: channel aspect ratio (=width/height)
Optional parameters:
* satTransport : A dictionary with the keys 'mu_f','mu_g,'rho_f','rho_g', 'sigma' for the saturation properties. So they can be calculated once and passed in for a slight improvement in efficiency
"""
def KMFunc(x):
dpdz, h = Kim_Mudawar_condensing_DPDZ_h(Ref,G,Dh,x,Tbubble,Tdew,p,beta,satTransport)
return dpdz , h
## Use Simpson's Rule to calculate the average pressure gradient
## Can't use adapative quadrature since function is not sufficiently smooth
## Not clear why not sufficiently smooth at x>0.9
if x_min==x_max:
return KMFunc(x_min)
else:
#Calculate the tranport properties once
satTransport={}
satTransport['rho_f']=PropsSI('D','T',Tbubble,'Q',0.0,Ref)
satTransport['rho_g']=PropsSI('D','T',Tdew,'Q',1.0,Ref)
satTransport['mu_f']=PropsSI('V','T',Tbubble,'Q',0.0,Ref)
satTransport['mu_g']=PropsSI('V','T',Tdew,'Q',1.0,Ref)
satTransport['cp_f']=PropsSI('C', 'T', Tbubble, 'Q', 0, Ref)
satTransport['k_f']=PropsSI('L', 'T', Tbubble, 'Q', 0, Ref)
#Calculate Dp and h over the range of xx
xx=np.linspace(x_min,x_max,30)
DP=np.zeros_like(xx)
h=np.zeros_like(xx)
for i in range(len(xx)):
DP[i]=KMFunc(xx[i])[0]
h[i]=KMFunc(xx[i])[1]
#Use Simpson's rule to carry out numerical integration to get average DP and average h
if abs(x_max-x_min)<5*machine_eps:
#return just one of the edge values
return DP[0], h[0]
else:
#Use Simpson's rule to carry out numerical integration to get average DP and average h
return -simps(DP,xx)/(x_max-x_min), simps(h,xx)/(x_max-x_min)
def KM_Cond_Average(x_min,x_max,Ref,G,Dh,Tbubble,Tdew,p,beta,satTransport=None):
"""
Returns the average pressure gradient and average heat transfer coefficient
between qualities of x_min and x_max.
for Kim&Mudawar two-phase condensation in mico-channel HX
To obtain the pressure gradient for a given value of x, pass it in as x_min and x_max
Required parameters:
* x_min : The minimum quality for the range [-]
* x_max : The maximum quality for the range [-]
* Ref : String with the refrigerant name
* G : Mass flux [kg/m^2/s]
* Dh : Hydraulic diameter of tube [m]
* Tbubble : Bubblepoint temperature of refrigerant [K]
* Tdew : Dewpoint temperature of refrigerant [K]
* beta: channel aspect ratio (=width/height)
Optional parameters:
* satTransport : A dictionary with the keys 'mu_f','mu_g,'rho_f','rho_g', 'sigma' for the saturation properties. So they can be calculated once and passed in for a slight improvement in efficiency
"""
def KMFunc(x):
dpdz, h = Kim_Mudawar_condensing_DPDZ_h(Ref,G,Dh,x,Tbubble,Tdew,p,beta,satTransport)
return dpdz , h
## Use Simpson's Rule to calculate the average pressure gradient
## Can't use adapative quadrature since function is not sufficiently smooth
## Not clear why not sufficiently smooth at x>0.9
if x_min==x_max:
return KMFunc(x_min)
else:
#Calculate the tranport properties once
satTransport={}
satTransport['rho_f']=PropsSI('D','T',Tbubble,'Q',0.0,Ref)
satTransport['rho_g']=PropsSI('D','T',Tdew,'Q',1.0,Ref)
satTransport['mu_f']=PropsSI('V','T',Tbubble,'Q',0.0,Ref)
satTransport['mu_g']=PropsSI('V','T',Tdew,'Q',1.0,Ref)
satTransport['cp_f']=PropsSI('C', 'T', Tbubble, 'Q', 0, Ref)
satTransport['k_f']=PropsSI('L', 'T', Tbubble, 'Q', 0, Ref)
#Calculate Dp and h over the range of xx
xx=np.linspace(x_min,x_max,30)
DP=np.zeros_like(xx)
h=np.zeros_like(xx)
for i in range(len(xx)):
DP[i]=KMFunc(xx[i])[0]
h[i]=KMFunc(xx[i])[1]
#Use Simpson's rule to carry out numerical integration to get average DP and average h
if abs(x_max-x_min)<5*machine_eps:
#return just one of the edge values
return DP[0], h[0]
else:
#Use Simpson's rule to carry out numerical integration to get average DP and average h
return -simps(DP,xx)/(x_max-x_min), simps(h,xx)/(x_max-x_min)
def energy_radiated_farfield(self, trajectory, gamma, x, y, distance):
# N = trajectory.shape[1]
N = trajectory.nb_points()
if distance == None:
# in radian :
n_chap = np.array([x, y, 1.0 - 0.5 * (x ** 2 + y ** 2)])
X = np.sqrt(x ** 2 + y ** 2 )#TODO a changer
#in meters :
else :
X = np.sqrt(x ** 2 + y ** 2 + distance ** 2)
n_chap = np.array([x, y, distance]) / X
R= X/codata.c - n_chap[0] * trajectory.x- n_chap[1] * trajectory.y - n_chap[2] * trajectory.z
E = np.zeros((3,), dtype=np.complex)
integrand = np.zeros((3, N), dtype=np.complex)
A1 = (n_chap[1] * trajectory.v_z - n_chap[2] * trajectory.v_y)
A2 = (-n_chap[0] * trajectory.v_z + n_chap[2] * trajectory.v_x)
A3 = (n_chap[0] * trajectory.v_y - n_chap[1] * trajectory.v_x)
Alpha1 = (1.0 / (1.0 - n_chap[0] * trajectory.v_x
- n_chap[1] * trajectory.v_y - n_chap[2] * trajectory.v_z)) ** 2
Alpha2 = np.exp(
0. + 1j * self.photon_frequency * (trajectory.t + R))
cst = codata.c / ((gamma ** 2) * R)
integrand[0] -= ((n_chap[1] * A3 - n_chap[2] * A2) * self.photon_frequency * 1j
+ cst * (n_chap[0] - trajectory.v_x) * Alpha1) * Alpha2
integrand[1] -= ((- n_chap[0] * A3 + n_chap[2] * A1) * self.photon_frequency * 1j
+ cst * (n_chap[1] - trajectory.v_y) * Alpha1) * Alpha2
integrand[2] -= ((n_chap[0] * A2 - n_chap[1] * A1) * self.photon_frequency * 1j
+ cst * (n_chap[2] - trajectory.v_z) * Alpha1) * Alpha2
for k in range(3):
E[k] = np.trapz(integrand[k], trajectory.t)
#E[k] = integrate.simps(integrand[k], trajectory.t)
terme_bord = np.full((3), 0. + 1j * 0., dtype=np.complex)
Alpha_1 = (1.0 / (1.0 - n_chap[0][-1] * trajectory.v_x[-1]
- n_chap[1][-1] * trajectory.v_y[-1] - n_chap[2][-1] * trajectory.v_z[-1]))
Alpha_0 = (1.0 / (1.0 - n_chap[0][0] * trajectory.v_x[0]
- n_chap[1][0] * trajectory.v_y[0] - n_chap[2][0] * trajectory.v_z[0]))
terme_bord += ((n_chap[1][-1] * A3[-1] - n_chap[2][-1] * A2[-1]) * Alpha_1 *
np.exp(1j * self.photon_frequency * (trajectory.t[-1] + R[-1] / codata.c)))
terme_bord -= ((n_chap[1][0] * A3[0] - n_chap[2][0] * A2[0]) * Alpha_0 *
np.exp(1j * self.photon_frequency * (trajectory.t[0] + R[0] / codata.c)))
E += terme_bord
return E
# exact equation for the energy radiated
# warning !!!!!!! far field approximation
def energy_radiated_near_field(self, trajectory, gamma, x, y, distance):
# N = trajectory.shape[1]
N = trajectory.nb_points()
n_chap = np.array(
[x - trajectory.x * codata.c, y - trajectory.y * codata.c, distance - trajectory.z * codata.c])
R = np.sqrt(n_chap[0] ** 2 + n_chap[1] ** 2 + n_chap[2] ** 2)
n_chap[0, :] /= R
n_chap[1, :] /= R
n_chap[2, :] /= R
E = np.zeros((3,), dtype=np.complex)
integrand = np.zeros((3, N), dtype=np.complex)
A1 = (n_chap[1] * trajectory.v_z - n_chap[2] * trajectory.v_y)
A2 = (-n_chap[0] * trajectory.v_z + n_chap[2] * trajectory.v_x)
A3 = (n_chap[0] * trajectory.v_y - n_chap[1] * trajectory.v_x)
Alpha1 = np.exp(
0. + 1j * self.photon_frequency * (trajectory.t + R/codata.c))
Alpha2 = codata.c / (gamma ** 2 * R)
integrand[0] -= ((n_chap[1] * A3 - n_chap[2] * A2) * self.photon_frequency* 1j
+ Alpha2 * (n_chap[0] - trajectory.v_x)
) * Alpha1
integrand[1] -= ((-n_chap[0] * A3 + n_chap[2] * A1) * self.photon_frequency * 1j
+ Alpha2 * (n_chap[1] - trajectory.v_y)
) * Alpha1
integrand[2] -= ((n_chap[0] * A2 - n_chap[1] * A1) * self.photon_frequency * 1j
+ Alpha2 * (n_chap[2] - trajectory.v_z)
) * Alpha1
for k in range(3):
E[k] = np.trapz(integrand[k], trajectory.t)
#E[k] = integrate.simps(integrand[k], trajectory.t)
terme_bord = np.full((3), 0. + 1j * 0., dtype=np.complex)
Alpha_1 = (1.0 / (1.0 - n_chap[0][-1] * trajectory.v_x[-1]
- n_chap[1][-1] * trajectory.v_y[-1] - n_chap[2][-1] * trajectory.v_z[-1]))
Alpha_0 = (1.0 / (1.0 - n_chap[0][0] * trajectory.v_x[0]
- n_chap[1][0] * trajectory.v_y[0] - n_chap[2][0] * trajectory.v_z[0]))
terme_bord += ((n_chap[1][-1] * A3[-1] - n_chap[2][-1] * A2[-1]) * Alpha_1*
np.exp(1j * self.photon_frequency * (trajectory.t[-1] +R[-1]/codata.c) ))
terme_bord -= ((n_chap[1][0] * A3[0] - n_chap[2][0] * A2[0]) * Alpha_0*
np.exp(1j * self.photon_frequency * (trajectory.t[0] + R[0]/codata.c)))
# E += terme_bord
return E
