Python numpy.log10函数代码示例

    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        """
        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
        for spec of input and result values.
        """

        # Distance term
        R = np.sqrt(dists.rjb ** 2 + 11.29 ** 2)

        # Magnitude term
        M = rup.mag - 6

        # Site term only distinguishes between lava and ash;
        # since ash sites have Vs30 in the range 60-200m/s,
        # we use this upper value as class separator
        S = np.zeros(R.shape)
        S[sites.vs30 <= 200] = 1

        # Mean ground motion (log10)
        mean = (0.518 + 0.387*M - np.log10(R) - 0.00256*R + 0.335*S)

        # Converting to natural log
        mean /= np.log10(np.e)

        # Check for standard deviation type
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        # Constant (total) standard deviation
        stddevs = [0.237/np.log10(np.e) + np.zeros(R.shape)]

        return mean, stddevs

def cd_sphere_vector(Re):
    "Computes the drag coefficient of a sphere as a function of the Reynolds number Re."
    # Curve fitted after fig . A -56 in Evett & Liu :% " Fluid Mechanics & Hydraulics ",
    # Schaum ' s Solved Problems McGraw - Hill 1989.

    from numpy import log10,array,polyval
    CD = zeros_like(Re)
   
    CD = where(Re<0,0.0,0.0)     # condition 1
    
    CD = where((Re > 0.0) & (Re <=0.5),24/Re,CD) # condition 2

    p = array([4.22,-14.05,34.87,0.658])
    CD = where((Re > 0.5) & (Re <=100.0),polyval(p,1.0/Re),CD) #condition 3

    p = array([-30.41,43.72,-17.08,2.41])
    CD = where((Re >100.0)  & (Re <=1.0e4) ,polyval(p,1.0/log10(Re)),CD) #condition 4

    p = array([-0.1584,2.031,-8.472,11.932])
    CD = where((Re > 1.0e4)  &  (Re <=3.35e5),polyval(p,log10(Re)),CD) #condition 5

    CD = where((Re > 3.35e5) & (Re <=5.0e5),91.08*(log10(Re/4.5e5))**4 + 0.0764,CD) #condition 6

    p  = array([-0.06338,1.1905,-7.332,14.93])
    CD = where((Re > 5.05e5)  &  (Re <=8.0e6),polyval(p,log10(Re)),CD) #condition 7
    
    CD = where(Re>8.0e6,0.2,CD)  # condition 8

    return CD

def compute_angular_momentum_transfer(snapshot, Rmin, Rmax, NR = None, Nphi = None, alpha=0.1, h0=0.1):
    """
    Compute the angular momentum transfer as a function of
    radius from a simulation snapshot and a prescribed grid
    data structure containing the coordinates of grid points
    where the torque balance should be evaluated.

    """
    

    if (NR is None):
        NR = 512
    if (Nphi is None):
        Nphi = int(2 * np.pi / (10**((np.log10(Rmax) - np.log10(Rmin))/NR) - 1))

    # Create a polar grid
    grid = grid_polar(NR = NR, Nphi = Nphi, Rmin= Rmin,Rmax = Rmax,scale='log')

    mdot = mass_advection(snapshot,grid)
    
    torque_adv = angular_momentum_advection(snapshot,grid)

    torque_visc = angular_momentum_viscosity(snapshot,grid, alpha = alpha, h0 = h0)

    torque_grav = angular_momentum_gravity(snapshot,grid)
    
    return grid.R.mean(axis=0),mdot, torque_adv,torque_visc,torque_grav

  def __init__(self,fname=None,mindx=None,mlist=[],logopt = 0, Lbox=2750., massfxnfname=None):
    """
    Read in a list of masses.  If logopt==1, then input masses are understood to be logarithmic.
    """
    self.Lbox = Lbox
    if mlist != []:
      if logopt == 0:
        self.m = np.array(mlist)
        self.lg10m = np.log10(self.m)
      else:
        self.lg10m = np.array(mlist)
        self.m = 10**(self.log10m)
      self.lg10mcen, self.Nofm = None, None ## set these later with massfxn.
    elif massfxnfname is not None:
      self.lg10mcen, self.Nofm = np.loadtxt(massfxnfname,unpack=True,usecols=[0,1])


    else:
      if 0==0:
#      try:
        if logopt == 0:
          self.m = np.loadtxt(fname,usecols=[mindx],unpack=True)
          self.lg10m = np.log10(self.m)
        else:
          self.lg10m = np.loadtxt(fname,usecols=[mindx],unpack=True)
          self.m = 10**(self.lg10m)
        self.lg10mcen, self.Nofm = None, None ## set these later with massfxn.
      else:
#      except:
        print 'file read did not work.'
        self.m = None
        self.lg10m = None

def get_cs(e_0=100, z=74):
    """
    Returns a function representing the scaled bremsstrahlung cross_section.

    Args:
        e_0 (float): The electron kinetic energy, used to scale u=e_e/e_0.
        z (int): Atomic number of the material.

    Returns:
        A function representing cross_section(e_g,u) in mb/keV, with e_g in keV.

    """
    # NOTE: Data is given for E0>1keV. CS values below this level should be used with caution.
    # The default behaviour is to keep it constant
    with open(os.path.join(data_path, "cs/grid.csv"), 'r') as csvfile:
        r = csv.reader(csvfile, delimiter=' ', quotechar='|',
                       quoting=csv.QUOTE_MINIMAL)
        t = next(r)
        e_e = np.array([float(a) for a in t[0].split(",")])
        log_e_e = np.log10(e_e)
        t = next(r)
        k = np.array([float(a) for a in t[0].split(",")])
    t = []
    with open(os.path.join(data_path, "cs/%d.csv" % z), 'r') as csvfile:
        r = csv.reader(csvfile, delimiter=' ', quotechar='|',
                       quoting=csv.QUOTE_MINIMAL)
        for row in r:
            t.append([float(a) for a in row[0].split(",")])
    t = np.array(t)
    scaled = interpolate.RectBivariateSpline(log_e_e, k, t, kx=3, ky=1)
    m_electron = 511
    z2 = z * z
    return lambda e_g, u: (u * e_0 + m_electron) ** 2 * z2 / (u * e_0 * e_g * (u * e_0 + 2 * m_electron)) * (
        scaled(np.log10(u * e_0), e_g / (u * e_0)))

def plot_hist(axis, data_list, label_list, logx, logy, overlaid):
    """
    """

    if logx:
        # Setup the logarithmic scale on the X axis
        data_array = np.array(data_list)
        vmin = np.log10(data_array.min())
        vmax = np.log10(data_array.max())
        bins = np.logspace(vmin, vmax, 50) # Make a range from 10**vmin to 10**vmax
    else:
        bins = 50

    if overlaid:
        for data_array, label in zip(data_list, label_list):
            res_tuple = axis.hist(data_array,
                                  bins=bins,
                                  log=logy,           # Set log scale on the Y axis
                                  histtype=HIST_TYPE,
                                  alpha=ALPHA,
                                  label=label)
    else:
        res_tuple = axis.hist(data_list,
                              bins=bins,
                              log=logy,               # Set log scale on the Y axis
                              histtype=HIST_TYPE,
                              alpha=ALPHA,
                              label=label_list)

	def plotPSD(self, chan, time_interval):
		Npackets = np.int(time_interval * self.accum_freq)
		plot_range = (Npackets / 2) + 1
		figure = plt.figure(num= None, figsize=(12,12), dpi=80, facecolor='w', edgecolor='w')
		# I 
		plt.suptitle('Channel ' + str(chan) + ' , Freq = ' + str((self.freqs[chan] + self.LO_freq)/1.0e6) + ' MHz') 
		plot1 = figure.add_subplot(311)
		plot1.set_xscale('log')
		plot1.set_autoscale_on(True)
		plt.ylim((-160,-80))
		plt.title('I')
		line1, = plot1.plot(np.linspace(0, self.accum_freq/2., (Npackets/2) + 1), np.zeros(plot_range), label = 'I', color = 'green', linewidth = 1)
		plt.grid()
		# Q
		plot2 = figure.add_subplot(312)
		plot2.set_xscale('log')
		plot2.set_autoscale_on(True)
		plt.ylim((-160,-80))
		plt.title('Q')
		line2, = plot2.plot(np.linspace(0, self.accum_freq/2., (Npackets/2) + 1), np.zeros(plot_range), label = 'Q', color = 'red', linewidth = 1)
		plt.grid()
		# Phase
		plot3 = figure.add_subplot(313)
		plot3.set_xscale('log')
		plot3.set_autoscale_on(True)
		plt.ylim((-120,-70))
		#plt.xlim((0.0001, self.accum_freq/2.))
		plt.title('Phase')
		plt.ylabel('dBc rad^2/Hz')
		plt.xlabel('log Hz')
		line3, = plot3.plot(np.linspace(0, self.accum_freq/2., (Npackets/2) + 1), np.zeros(plot_range), label = 'Phase', color = 'black', linewidth = 1)
		plt.grid()
		plt.show(block = False)
		count = 0
		stop = 1.0e10
		while count < stop:
			Is, Qs, phases = self.get_stream(chan, time_interval)
			I_mags = np.fft.rfft(Is, Npackets)
			Q_mags = np.fft.rfft(Is, Npackets)
			phase_mags = np.fft.rfft(phases, Npackets)
			I_vals = (np.abs(I_mags)**2 * ((1./self.accum_freq)**2 / (1.0*time_interval)))
			Q_vals = (np.abs(Q_mags)**2 * ((1./self.accum_freq)**2 / (1.0*time_interval)))
			phase_vals = (np.abs(phase_mags)**2 * ((1./self.accum_freq)**2 / (1.0*time_interval)))
			phase_vals = 10*np.log10(phase_vals)
			phase_vals -= phase_vals[0]
			#line1.set_ydata(Is)
			#line2.set_ydata(Qs)
			#line3.set_ydata(phases)
			line1.set_ydata(10*np.log10(I_vals))
			line2.set_ydata(10*np.log10(Q_vals))
			line3.set_ydata(phase_vals)
			plot1.relim()
			plot1.autoscale_view(True,True,False)
			plot2.relim()
			plot2.autoscale_view(True,True,False)
			#plot3.relim()
			plot3.autoscale_view(True,True,False)
			plt.draw()
			count +=1
		return

    def __repr__(self):
        pars = ""
        if not isinstance(self.parameters, OrderedDict):
            keys = sorted(self.parameters.keys()) # to ensure the order is always the same
        else:
            keys = self.parameters.keys()

        for k in keys:
            v = self.parameters[k].value
            par_fmt = "{}"
            post = ""

            if hasattr(v, 'unit'):
                post = " {}".format(v.unit)
                v = v.value

            if isinstance(v, float):
                if v == 0:
                    par_fmt = "{:.0f}"
                elif np.log10(v) < -2 or np.log10(v) > 5:
                    par_fmt = "{:.2e}"
                else:
                    par_fmt = "{:.2f}"

            elif isinstance(v, int) and np.log10(v) > 5:
                par_fmt = "{:.2e}"

            pars += ("{}=" + par_fmt + post).format(k, v) + ", "

        if isinstance(self.units, DimensionlessUnitSystem):
            return "<{}: {} (dimensionless)>".format(self.__class__.__name__, pars.rstrip(", "))
        else:
            return "<{}: {} ({})>".format(self.__class__.__name__, pars.rstrip(", "), ",".join(map(str, self.units._core_units)))

def compatibility(par_low, par_high):
    """Quantify spectral compatibility of power-law
    measurements in two energy bands.

    Reference: 2008ApJ...679.1299F Equation (2)

    Compute spectral compatibility parameters for the
    situation where two power laws were measured in a low
    and a high spectral energy band.
    par_low and par_high are the measured parameters,
    which must be lists in the following order:
    e, f, f_err, g, g_err
    where e is the pivot energy, f is the flux density
    and g the spectral index
    """
    # Unpack power-law paramters
    e_high, f_high, f_err_high, g_high, g_err_high = par_high
    e_low, f_low, f_err_low, g_low, g_err_low = par_low

    log_delta_e = np.log10(e_high) - np.log10(e_low)
    log_delta_f = np.log10(f_high) - np.log10(f_low)
    # g_match is the index obtained by connecting the two points
    # with a power law, i.e. a straight line in the log_e, log_f plot
    g_match = -log_delta_f / log_delta_e

    # sigma is the number of standar deviations the match index
    # is different from the measured index in one band.
    # (see Funk et al. (2008ApJ...679.1299F) eqn. 2)
    sigma_low = (g_match - g_low) / g_err_low
    sigma_high = (g_match - g_high) / g_err_high
    sigma_comb = np.sqrt(sigma_low ** 2 + sigma_high ** 2)

    return g_match, sigma_low, sigma_high, sigma_comb

def _default_response_frequencies(A, n):
    """Compute a reasonable set of frequency points for bode plot.

    This function is used by `bode` to compute the frequency points (in rad/s)
    when the `w` argument to the function is None.

    Parameters
    ----------
    A : ndarray
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    w : ndarray
        The 1-D array of length `n` of frequency samples (in rad/s) at which
        the response is to be computed.
    """
    vals = linalg.eigvals(A)
    # Remove poles at 0 because they don't help us determine an interesting
    # frequency range. (And if we pass a 0 to log10() below we will crash.)
    poles = [pole for pole in vals if pole != 0]
    # If there are no non-zero poles, just hardcode something.
    if len(poles) == 0:
        minpole = 1
        maxpole = 1
    else:
        minpole = min(abs(real(poles)))
        maxpole = max(abs(real(poles)))
    # A reasonable frequency range is two orders of magnitude before the
    # minimum pole (slowest) and two orders of magnitude after the maximum pole
    # (fastest).
    w = numpy.logspace(numpy.log10(minpole) - 2, numpy.log10(maxpole) + 2, n)
    return w

 def error_vs_updates(self, output_dir=None):
     fig = plt.figure()
     ax = fig.add_subplot(111)
     ax.scatter(log10(centered(self.w_truth)) - log10(centered(self.medians)), self.updates, c=self.stds, cmap=plt.jet(), s=25, clip_on=False, lw=0.5)
     ax.set_xlabel('log10(w_truth)-log10(median w)')
     ax.set_ylabel('num updates')
     show(fig, output_dir, 'error_vs_updates.png')

def create_center_frequencies(stt=180, stp=7000, n_bands=32, kind='log'):
    '''
    Define center frequencies for spectrograms.

    Generally this is for auditory spectrogram extraction. Most auditory
    analysis uses 180 - 7000 Hz, so for now those
    are the defaults.

    Parameters
    ----------
    stt : float | int
        The starting frequency
    stp : float | int
        The end frequency
    n_bands : int
        The number of bands to calculate
    kind : 'log' | 'erb'
        Whether to use log or erb spacing

    Returns
    -------
    freqs : array, shape (n_frequencies,)
        An array of center frequencies.
    '''
    if kind == 'log':
        freqs = np.logspace(np.log10(stt), np.log10(stp), n_bands).astype(int)
    elif kind == 'erb':
        freqs = hears.erbspace(stt * Hz, stp * Hz, n_bands)
    else:
        print("I don't know what kind of spacing that is")
    return freqs

    def fluence_dist(self):
        """ Plots the fluence distribution and gives the mean and median fluence
        values of the sample """
        fluences = []
        for i in range(0,len(self.fluences),1):
            try:
                fluences.append(float(self.fluences[i]))

            except ValueError:
                continue

        fluences = np.array(fluences)
        mean_fluence = np.mean(fluences)
        median_fluence = np.median(fluences)
        print('Mean Fluence =',mean_fluence,'(15-150 keV) [10^-7 erg cm^-2]')
        print('Median Fluence =',median_fluence,'(15-150 keV) [10^-7 erg cm^-2]')

        plt.figure()
        plt.xlabel('Fluence (15-150 keV) [$10^{-7}$ erg cm$^{-2}$]')
        plt.ylabel('Number of GRBs')
        plt.xscale('log')
        minimum, maximum = min(fluences), max(fluences)
        plt.axvline(mean_fluence,color='red',linestyle='-')
        plt.axvline(median_fluence,color='blue',linestyle='-')
        plt.hist(fluences,bins= 10**np.linspace(np.log10(minimum),np.log10(maximum),20),color='grey',alpha=0.5)
        plt.show()

    def t90_dist(self):
        """ Plots T90 distribution, gives the mean and median T90 values of the
        sample and calculates the number of short, long bursts in the sample """
        t90s = []
        for i in range(0,len(self.t90s),1):
            try:
                t90s.append(float(self.t90s[i]))

            except ValueError:
                continue

        t90s = np.array(t90s)
        mean_t90 = np.mean(t90s)
        median_t90 = np.median(t90s)
        print('Mean T90 time =',mean_t90,'s')
        print('Median T90 time=',median_t90,'s')
        mask = np.ma.masked_where(t90s < 2, t90s)
        short_t90s = t90s[mask == False]
        long_t90s = t90s[mask != False]
        print('Number of Short/Long GRBs =',len(short_t90s),'/',len(long_t90s))

        plt.figure()
        plt.xlabel('T$_{90}$ (s)')
        plt.ylabel('Number of GRBs')
        plt.xscale('log')
        minimum, maximum, = min(short_t90s), max(long_t90s)
        plt.axvline(mean_t90,color='red',linestyle='-')
        plt.axvline(median_t90,color='blue',linestyle='-')
        plt.hist(t90s,bins= 10**np.linspace(np.log10(minimum),np.log10(maximum),20),color='grey',alpha=0.5)
        plt.show()

def pltytfield():
    fig=plt.figure()
    ppy=yt.ProjectionPlot(ds, "x", "Bxy", weight_field="density") #Project X-component of B-field from z-direction
    By=ppy._frb["Bxy"]
    ax=fig.add_subplot(111)
    plt.xticks(tick_locs,tick_lbls)
    plt.yticks(tick_locs,tick_lbls)
    Bymag=ax.pcolormesh(np.log10(By))
    cbar_m=plt.colorbar(Bymag)
    cbar_m.set_label("Bxy")
    plt.title("Bxy in yz plane")

    fig=plt.figure()
    ppy=yt.ProjectionPlot(ds, "y", "Bxy", weight_field="density") #Project X-component of B-field from z-direction
    By=ppy._frb["Bxy"]
    ax=fig.add_subplot(111)
    plt.xticks(tick_locs,tick_lbls)
    plt.yticks(tick_locs,tick_lbls)
    Bymag=ax.pcolormesh(np.log10(By))
    cbar_m=plt.colorbar(Bymag)
    cbar_m.set_label("Bxy")
    plt.title("Bxy in xz plane")

    fig=plt.figure()
    ppy=yt.ProjectionPlot(ds, "z", "Bxy", weight_field="density") #Project X-component of B-field from z-direction
    By=ppy._frb["Bxy"]
    ax=fig.add_subplot(111)
    plt.xticks(tick_locs,tick_lbls)
    plt.yticks(tick_locs,tick_lbls)
    Bymag=ax.pcolormesh(np.log10(By))
    cbar_m=plt.colorbar(Bymag)
    cbar_m.set_label("Bxy")
    plt.title("Bxy in xy plane")