Returns asteroseismic nu_max in uHz

        reference: https://arxiv.org/pdf/1312.3853v1.pdf, Eq (3)
484765
def nu_max(self, *args):
        
        return 3120.* (self.mass(*args) /
                        (self.radius(*args)**2 * np.sqrt(self.Teff(*args)/5777.)))
254166
def rotfreq_to_omega(rotfreq, scale=c.R_sun.si.value, solar_units=False):
    
    if solar_units:
        omega = rotfreq / (2*np.pi) / np.sqrt(c.GM_sun.to(u.solRad**3/u.d**2).value/scale**3)
    else:
        # then SI units
        omega = rotfreq / (2*np.pi) / np.sqrt(c.GM_sun.value/scale**3)

    # print "*** rotstar.rotfreq_to_omega", rotfreq, scale, solar_units, omega

    return omega
248763
def max(self, axis=None, skipna=True):
        """"""
        nv.validate_minmax_axis(axis)
        return self._minmax('max')
443128
def _get_2e_asym(self):
        
        
        # get needed PPG parameters for splitting 1D histos into 2D histos
        try:        
            # get frequency vector
            freq = np.arange(self.ppg['freq_start'].mean,\
                             self.ppg['freq_stop'].mean+\
                                    self.ppg['freq_incr'].mean,\
                             self.ppg['freq_incr'].mean)
                             
            # number of dwelltimes per frequency bin 
            ndwell = 2*int(self.ppg['ndwell_per_f'].mean)-1
            
            # number of RF on delays for the start bin. 
            start_bin = int(self.ppg['rf_on_delay'].mean)
            
            # get bin centers in ms
            time = self.ppg['rf_on_ms'].mean*(np.arange(ndwell)+0.5-ndwell/2.)
            
            # get the time and index of the middle time 
            mid_time_i = int(np.floor(ndwell/2.))
            mid_time = time[mid_time_i]
        
            # beam off time after pulse in ms
            beam_off = int(self.ppg['beam_off_ms'].mean)
        
        except KeyError:
            raise RuntimeError("Not all dictionary variables read out to "+\
                               "proper locations")
            
        # setup output
        out = bdict()
        out['freq'] = freq
        out['time'] = time
            
        # get data
        data = np.array(self._get_area_data()) # [[fp], [bfm], [bp], [bm]]
        
        # discared initial bad bins, and beam-off trailing bins
        data = data[:,start_bin:len(freq)*ndwell+start_bin]
        
        # split data by frequency
        nsplit = len(data[0])/ndwell
        
        fp = np.array(np.split(data[0],nsplit))
        fm = np.array(np.split(data[1],nsplit))
        bp = np.array(np.split(data[2],nsplit))
        bm = np.array(np.split(data[3],nsplit))
        
        # get raw asymmetries 
        asym_p_2cntr = (bp-fp)/(bp+fp)      # two counter
        asym_m_2cntr = (bm-fm)/(bm+fm)      # two counter
        r = np.sqrt(bp*fm/(bm*fp))
        asym_4cntr = (r-1)/(r+1)            # four counter
        
        # get raw asymmetry errors
        asym_p_2cntr_err = 2*np.sqrt(bp*fp)/((bp+fp)**1.5)
        asym_m_2cntr_err = 2*np.sqrt(bm*fm)/((bm+fm)**1.5)
        asym_4cntr_err = r*np.sqrt(1./bp+1./bm+1./fp+1./fm)/((r+1)**2)

        # save to output
        out['raw_p'] = np.array([asym_p_2cntr,asym_p_2cntr_err])
        out['raw_n'] = np.array([asym_m_2cntr,asym_m_2cntr_err])
        out['raw_c'] = np.array([asym_4cntr,asym_4cntr_err])

        # wrap asymmetry arrays into one for calculations [p,m,4]
        # indexing is now [pm4][freq][time bin]
        asym     = np.array([asym_p_2cntr,    asym_m_2cntr,    asym_4cntr])
        asym_err = np.array([asym_p_2cntr_err,asym_m_2cntr_err,asym_4cntr_err])

        # compute differenced asymmetries via slopes from weighted least squares 
        # minimization.
        if ndwell >= 5:
            
            # calculate needed components element-wise
            w = asym_err**-2
            x = time
            y = asym
            
            wx = w*x
            wy = w*y
            wxy = w*x*y
            wxx = w*x*x
        
            # sum over values i < mid_time_i within each asymmetry and frequency 
            # Indexing: [pm4][freq]
            w_pre   = np.sum(w  [:,:,:mid_time_i],2)    
            wx_pre  = np.sum(wx [:,:,:mid_time_i],2)
            wy_pre  = np.sum(wy [:,:,:mid_time_i],2)
            wxy_pre = np.sum(wxy[:,:,:mid_time_i],2)
            wxx_pre = np.sum(wxx[:,:,:mid_time_i],2)
            
            # sum over values i > mid_time_i
            w_pst   = np.sum(w  [:,:,-mid_time_i:],2)
            wx_pst  = np.sum(wx [:,:,-mid_time_i:],2)
            wy_pst  = np.sum(wy [:,:,-mid_time_i:],2)
            wxy_pst = np.sum(wxy[:,:,-mid_time_i:],2)
            wxx_pst = np.sum(wxx[:,:,-mid_time_i:],2)
            
            # calculate slopes and intercepts
            delta_pre = w_pre*wxx_pre - wx_pre**2
            delta_pst = w_pst*wxx_pst - wx_pst**2
            
            sl_pre = (w_pre*wxy_pre - wx_pre*wy_pre)/delta_pre
            sl_pst = (w_pst*wxy_pst - wx_pst*wy_pst)/delta_pst
            dsl_pre = np.sqrt(w_pre/delta_pre)
            dsl_pst = np.sqrt(w_pst/delta_pst)
            
            intr_pre = (wy_pre*wxx_pre - wx_pre*wxy_pre)/delta_pre
            intr_pst = (wy_pst*wxx_pst - wx_pst*wxy_pst)/delta_pst
            dintr_pre = np.sqrt(wxx_pre/delta_pre)
            dintr_pst = np.sqrt(wxx_pst/delta_pst)
            
            # extrapolate to middle time bin
            asym_slopes = intr_pst-intr_pre+(sl_pst-sl_pre)*mid_time
            asym_slopes_err = np.sqrt(dintr_pre**2 + dintr_pst**2 + \
                                      (dsl_pre**2 + dsl_pst**2) * mid_time**2)
            
            # save to output        
            out['sl_p'] = np.array([asym_slopes[0],asym_slopes_err[0]])
            out['sl_n'] = np.array([asym_slopes[1],asym_slopes_err[1]])
            out['sl_c'] = np.array([asym_slopes[2],asym_slopes_err[2]])
        
        # calculate asymmetry using differences
        asym_diff = asym[:,:,mid_time_i+1] - asym[:,:,mid_time_i-1]
        asym_diff_err = np.sqrt(asym_err[:,:,mid_time_i+1]**2+\
                                asym_err[:,:,mid_time_i-1]**2)
        
        # save to output
        out['dif_p'] = np.array([asym_diff[0],asym_diff_err[0]])
        out['dif_n'] = np.array([asym_diff[1],asym_diff_err[1]])
        out['dif_c'] = np.array([asym_diff[2],asym_diff_err[2]])   
            
        return out
393415
def max(self):
        

        if len(self.regions) != 1:
            raise ClaripyVSAOperationError("'max()' onlly works on single-region value-sets.")

        return self.get_si(next(iter(self.regions))).max
15301
def max(self, axis=None, skipna=True, *args, **kwargs):
        
        # TODO: skipna is broken with max.
        # See https://github.com/pandas-dev/pandas/issues/24265
        nv.validate_max(args, kwargs)
        nv.validate_minmax_axis(axis)

        mask = self.isna()
        if skipna:
            values = self[~mask].asi8
        elif mask.any():
            return NaT
        else:
            values = self.asi8

        if not len(values):
            # short-circut for empty max / min
            return NaT

        result = nanops.nanmax(values, skipna=skipna)
        # Don't have to worry about NA `result`, since no NA went in.
        return self._box_func(result)
1137027
def conf_interval(x,L,conf=0.683,shortest=True,
                  conftol=0.001,return_max=False):
    
    cum = np.cumsum(L)
    cdf = cum/cum.max()
    if shortest:
        maxind = L.argmax()
        if maxind==0:   #hack alert
            maxind = 1
        if maxind==len(L)-1:
            maxind = len(L)-2
        Lval = L[maxind]

        lox = x[0:maxind]
        loL = L[0:maxind]
        locdf = cdf[0:maxind]
        hix = x[maxind:]
        hiL = L[maxind:]
        hicdf = cdf[maxind:]

        dp = 0
        s = -1
        dL = Lval
        switch = False
        last = 0
        while np.absolute(dp-conf) > conftol:
            Lval += s*dL
            if maxind==0:
                loind = 0
            else:
                loind = (np.absolute(loL - Lval)).argmin()
            if maxind==len(L)-1:
                hiind = -1
            else:
                hiind = (np.absolute(hiL - Lval)).argmin()

            dp = hicdf[hiind]-locdf[loind]
            lo = lox[loind]
            hi = hix[hiind]
            if dp == last:
                break
            last = dp
            cond = dp > conf
            if cond ^ switch:
                dL /= 2.
                s *= -1
                switch = not switch

    else:
        alpha = (1-conf)/2.
        lo = x[np.absolute(cdf-alpha).argmin()]
        hi = x[(np.absolute(cdf-(1-(alpha)))).argmin()]
        
    if return_max:
        xmaxL = x[L.argmax()]
        return xmaxL,lo,hi
    else:
        return (lo,hi)
239695
def get_template_amplitude_norm(template=None, **kwargs):
    
    input_params = props(template,**kwargs)
    approximant = kwargs['approximant']

    if approximant in _template_amplitude_norms:
        return _template_amplitude_norms[approximant](**input_params)
    else:
        return None
395334
def _vmomentdensity(self,R,z,n,m,o,nsigma=None,mc=False,nmc=10000,
                       _returnmc=False,_vrs=None,_vts=None,_vzs=None,
                       _rawgausssamples=False,
                       gl=False,ngl=_DEFAULTNGL,_returngl=False,_glqeval=None,
                       _return_actions=False,_jr=None,_lz=None,_jz=None,
                       _return_freqs=False,
                       _rg=None,_kappa=None,_nu=None,_Omega=None,
                       _sigmaR1=None,_sigmaz1=None,
                       **kwargs):
        """"""
        if isinstance(R,numpy.ndarray):
            return numpy.array([self._vmomentdensity(r,zz,n,m,o,nsigma=nsigma,
                                                    mc=mc,nmc=nmc,
                                                    gl=gl,ngl=ngl,**kwargs) for r,zz in zip(R,z)])
        if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                actionAngle.actionAngleAdiabaticGrid)):
            if n % 2 == 1. or o % 2 == 1.:
                return 0. #we know this must be the case
        if nsigma == None:
            nsigma= _NSIGMA
        if _sigmaR1 is None:
            sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr)
        else:
            sigmaR1= _sigmaR1
        if _sigmaz1 is None:
            sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz)
        else:
            sigmaz1= _sigmaz1
        thisvc= potential.vcirc(self._pot,R,use_physical=False)
        #Use the asymmetric drift equation to estimate va
        gamma= numpy.sqrt(0.5)
        va= sigmaR1**2./2./thisvc\
            *(gamma**2.-1. #Assume close to flat rotation curve, sigphi2/sigR2 =~ 0.5
               +R*(1./self._hr+2./self._hsr))
        if math.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center
        if gl:
            if ngl % 2 == 1:
                raise ValueError("ngl must be even")
            if not _glqeval is None and ngl != _glqeval.shape[0]:
                _glqeval= None
            #Use Gauss-Legendre integration for all
            if ngl == _DEFAULTNGL:
                glx, glw= self._glxdef, self._glwdef
                glx12, glw12= self._glxdef12, self._glwdef12
            elif ngl == _DEFAULTNGL2:
                glx, glw= self._glxdef2, self._glwdef2
                glx12, glw12= self._glxdef, self._glwdef
            else:
                glx, glw= numpy.polynomial.legendre.leggauss(ngl)
                glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2)
            #Evaluate everywhere
            if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                    actionAngle.actionAngleAdiabaticGrid)):
                vRgl= nsigma*sigmaR1/2.*(glx+1.)
                vzgl= nsigma*sigmaz1/2.*(glx+1.)
                vRglw= glw
                vzglw= glw
            else:
                vRgl= nsigma*sigmaR1/2.*(glx12+1.)
                #vRgl= 1.5/2.*(glx12+1.)
                vRgl= list(vRgl)
                vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.))
                #vRgl.extend(-1.5/2.*(glx12+1.))
                vRgl= numpy.array(vRgl)
                vzgl= nsigma*sigmaz1/2.*(glx12+1.)
                #vzgl= 1.5/2.*(glx12+1.)
                vzgl= list(vzgl)
                vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.))
                #vzgl.extend(-1.5/2.*(glx12+1.))
                vzgl= numpy.array(vzgl)
                vRglw= glw12
                vRglw= list(vRglw)
                vRglw.extend(glw12)
                vRglw= numpy.array(vRglw)
                vzglw= glw12
                vzglw= list(vzglw)
                vzglw.extend(glw12)
                vzglw= numpy.array(vzglw)
            if 'vTmax' in kwargs: vTmax = kwargs['vTmax']
            else:                 vTmax = 1.5
            vTgl= vTmax/2.*(glx+1.)
            #Tile everything
            vTgl= numpy.tile(vTgl,(ngl,ngl,1)).T
            vRgl= numpy.tile(numpy.reshape(vRgl,(1,ngl)).T,(ngl,1,ngl))
            vzgl= numpy.tile(vzgl,(ngl,ngl,1))
            vTglw= numpy.tile(glw,(ngl,ngl,1)).T #also tile weights
            vRglw= numpy.tile(numpy.reshape(vRglw,(1,ngl)).T,(ngl,1,ngl))
            vzglw= numpy.tile(vzglw,(ngl,ngl,1))
            #evaluate
            if _glqeval is None and _jr is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R+numpy.zeros(ngl*ngl*ngl),
                                           vRgl.flatten(),
                                           vTgl.flatten(),
                                           z+numpy.zeros(ngl*ngl*ngl),
                                           vzgl.flatten(),
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            elif not _jr is None and _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                                                 log=True,
                                                                 _return_actions=True,
                                                                 _return_freqs=True,
                                                                 use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            elif not _jr is None and not _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                           rg=_rg,kappa=_kappa,nu=_nu,
                                           Omega=_Omega,
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            else:
                logqeval= _glqeval
            if _returngl:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        logqeval)
            elif _return_actions and _return_freqs:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        jr,lz,jz,
                        rg,kappa,nu,Omega)
            elif _return_actions:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        jr,lz,jz)
            else:
                return numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                 *vTglw*vRglw*vzglw*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2)
        elif mc:
            mvT= (thisvc-va)/gamma/sigmaR1
            if _vrs is None:
                vrs= numpy.random.normal(size=nmc)
            else:
                vrs= _vrs
            if _vts is None:
                vts= numpy.random.normal(size=nmc)+mvT
            else:
                if _rawgausssamples:
                    vts= _vts+mvT
                else:
                    vts= _vts
            if _vzs is None:
                vzs= numpy.random.normal(size=nmc)
            else:
                vzs= _vzs
            Is= _vmomentsurfaceMCIntegrand(vzs,vrs,vts,numpy.ones(nmc)*R,
                                           numpy.ones(nmc)*z,
                                           self,sigmaR1,gamma,sigmaz1,mvT,
                                           n,m,o)
            if _returnmc:
                if _rawgausssamples:
                    return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o),
                        vrs,vts-mvT,vzs)
                else:
                    return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o),
                            vrs,vts,vzs)
            else:
                return numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
        else: #pragma: no cover because this is too slow; a warning is shown
            warnings.warn("Calculations using direct numerical integration using tplquad is not recommended and extremely slow; it has also not been carefully tested",galpyWarning)
            return integrate.tplquad(_vmomentsurfaceIntegrand,
                                     1./gamma*(thisvc-va)/sigmaR1-nsigma,
                                     1./gamma*(thisvc-va)/sigmaR1+nsigma,
                                     lambda x: 0., lambda x: nsigma,
                                     lambda x,y: 0., lambda x,y: nsigma,
                                     (R,z,self,sigmaR1,gamma,sigmaz1,n,m,o),
                                     **kwargs)[0]*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
743360
def max_intensity(self, time):
        
        ti = np.where(time == self.times)[0][0]
        return self.timesteps[ti].max()
1010334
def calcTemperature(self):
        
        try:
            return eq.MeanPlanetTemp(self.albedo, self.star.T, self.star.R, self.a).T_p
        except (ValueError, HierarchyError):  # ie missing value (.a) returning nan
            return np.nan
955204
def _get_pga_on_rock(self, C, rup, dists):
        
        return np.exp(self._get_magnitude_scaling_term(C, rup) +
                      self._get_path_scaling(C, dists, rup.mag))
239379
def sun_ra_dec(utc_time):
    
    jdate = jdays2000(utc_time) / 36525.0
    eps = np.deg2rad(23.0 + 26.0 / 60.0 + 21.448 / 3600.0 -
                     (46.8150 * jdate + 0.00059 * jdate * jdate -
                      0.001813 * jdate * jdate * jdate) / 3600)
    eclon = sun_ecliptic_longitude(utc_time)
    x__ = np.cos(eclon)
    y__ = np.cos(eps) * np.sin(eclon)
    z__ = np.sin(eps) * np.sin(eclon)
    r__ = np.sqrt(1.0 - z__ * z__)
    # sun declination
    declination = np.arctan2(z__, r__)
    # right ascension
    right_ascension = 2 * np.arctan2(y__, (x__ + r__))
    return right_ascension, declination
697161
private int getMaxOccurence(double massTo, int element_pos, int[] matrix, List<IIsotope> isoToCond_new) {
        double massIn = isoToCond_new.get(element_pos).getExactMass();
        double massToM = massTo;
        for (int i = 0; i < matrix.length; i++)
            if (i != element_pos) if (matrix[i] != 0) massToM -= isoToCond_new.get(i).getExactMass() * matrix[i];

        int value = (int) ((massToM + 1) / massIn);
        return value;
    }
1087483
private double extremum_alpha_n(int n, double[] alpha) {
    // arctan(infinity) = PI/2
    if(vec.doubleValue(n) == 0) {
      return MathUtil.HALFPI;
    }

    double tan = 0;
    for(int j = n + 1; j < vec.getDimensionality(); j++) {
      double alpha_j = j == vec.getDimensionality() - 1 ? 0 : alpha[j];
      tan += vec.doubleValue(j) * sinusProduct(n + 1, j, alpha) * FastMath.cos(alpha_j);
    }
    tan /= vec.doubleValue(n);

    // if (debug) {
    // debugFiner("tan alpha_" + (n + 1) + " = " + tan);
    // }
    double alpha_n = Math.atan(tan);
    if(alpha_n < 0) {
      alpha_n = Math.PI + alpha_n;
    }
    return alpha_n;
  }
473464
def canonical(self):
        
        Ac1 = np.hstack((self.deltah, np.zeros((self.nh, self.nz))))
        Ac2 = np.hstack((np.zeros((self.nz, self.nh)), self.a22))
        Ac = np.vstack((Ac1, Ac2))
        Bc = np.vstack((self.thetah, np.zeros((self.nz, self.nc))))
        Cc = np.vstack((np.zeros((self.nh, self.nw)), self.c2))
        Rc1 = np.hstack((self.llambda.T.dot(self.llambda), -
                         self.llambda.T.dot(self.ub)))
        Rc2 = np.hstack((-self.ub.T.dot(self.llambda), self.ub.T.dot(self.ub)))
        Rc = np.vstack((Rc1, Rc2))
        Qc = self.pih.T.dot(self.pih)
        Nc = np.hstack(
            (self.pih.T.dot(self.llambda), -self.pih.T.dot(self.ub)))

        lq_aux = LQ(Qc, Rc, Ac, Bc, N=Nc, beta=self.beta)

        P1, F1, d1 = lq_aux.stationary_values()

        self.F_b = F1[:, 0:self.nh]
        self.F_f = F1[:, self.nh:]

        self.pihat = np.linalg.cholesky(self.pih.T.dot(
            self.pih) + self.beta.dot(self.thetah.T).dot(P1[0:self.nh, 0:self.nh]).dot(self.thetah)).T
        self.llambdahat = self.pihat.dot(self.F_b)
        self.ubhat = - self.pihat.dot(self.F_f)

        return
754339
def lm_tfinal(damping_times, modes):
    

    t_max = {}
    for lmn in modes:
        l, m, nmodes = int(lmn[0]), int(lmn[1]), int(lmn[2])
        for n in range(nmodes):
            t_max['%d%d%d' %(l,m,n)] = \
            qnm_time_decay(damping_times['%d%d%d' %(l,m,n)], 1./1000)

    return max(t_max.values())