Compute hanging-wall effect (see eq. 7, 8, 9 and 10 page 319).
        Considers correct version of equation 8 as given in the erratum and not
        in the original paper.
505482
def _compute_hanging_wall_effect(self, C, rjb, rrup, dip, mag):
        
        # eq. 8 (to be noticed that the USGS-NSHMP implementation defines
        # the hanging-wall term for all rjb distances, while in the original
        # manuscript, hw is computed only for rjb < 5). Again the 'firm rock'
        # is considered
        hw = np.zeros_like(rjb)
        if dip <= 70.:
            hw = (5. - rjb) / 5.

        # eq. 9
        f_m = 1 if mag > 6.5 else mag - 5.5

#        # eq. 10
        f_rrup = C['c15'] + np.zeros_like(rrup)
        idx = rrup < 8
        f_rrup[idx] *= rrup[idx] / 8

        # eq. 7 (to be noticed that the f3 factor is not included
        # while this is defined in the original manuscript)
        f_hw = hw * f_m * f_rrup

        return f_hw
521078
def seebeck_eff_mass_from_carr(eta, n, T, Lambda):
    
    from fdint import fdk

    return (2 * np.pi**2 * abs(n) * 10 ** 6 / (fdk(0.5,eta))) ** (2. / 3)\
            / (2 * constants.m_e * constants.k * T / (constants.h/2/np.pi) ** 2)
1252120
function(t, b, c, d) {
            if(t < d / 2) {
                return Kinetic.Easings.BounceEaseIn(t * 2, 0, c, d) * 0.5 + b;
            }
            else {
                return Kinetic.Easings.BounceEaseOut(t * 2 - d, 0, c, d) * 0.5 + c * 0.5 + b;
            }
        }
159468
function erf( x ) {
	var sign;
	var ax;
	var z;
	var r;
	var s;
	var y;
	var p;
	var q;

	// Special case: NaN
	if ( x !== x ) {
		return NaN;
	}
	// Special case: +infinity
	if ( x === PINF ) {
		return 1;
	}
	// Special case: -infinity
	if ( x === NINF ) {
		return -1;
	}
	// Special case: +-0
	if ( x === 0 ) {
		return x;
	}
	if ( x < 0 ) {
		sign = true;
		ax = -x;
	} else {
		sign = false;
		ax = x;
	}
	// |x| < 0.84375
	if ( ax < 0.84375 ) {
		if ( ax < SMALL ) {
			if ( ax < VERY_TINY ) {
				// Avoid underflow:
				return 0.125 * (8.0*x + EFX8*x);
			}
			return x + EFX*x;
		}
		z = x * x;
		r = PPC + z*polyvalPP( z );
		s = QQC + z*polyvalQQ( z );
		y = r / s;
		return x + x*y;
	}
	// 0.84375 <= |x| < 1.25
	if ( ax < 1.25 ) {
		s = ax - 1;
		p = PAC + s*polyvalPA( s );
		q = QAC + s*polyvalQA( s );
		if ( sign ) {
			return -ERX - p/q;
		}
		return ERX + p/q;
	}
	// +inf > |x| >= 6
	if ( ax >= 6 ) {
		if ( sign ) {
			return TINY - 1.0; // raise inexact
		}
		return 1.0 - TINY; // raise inexact
	}
	s = 1.0 / (ax*ax);

	// |x| < 1/0.35 ~ 2.857143
	if ( ax < 2.857142857142857 ) {
		r = RAC + s*polyvalRA( s );
		s = SAC + s*polyvalSA( s );
	}
	// |x| >= 1/0.35 ~ 2.857143
	else {
		r = RBC + s*polyvalRB( s );
		s = SBC + s*polyvalSB( s );
	}
	z = setLowWord( ax, 0 ); // pseudo-single (20-bit) precision x
	r = exp( -z*z - 0.5625 ) * exp( (z-ax)*(z+ax) + r/s );
	if ( sign ) {
		return r/ax - 1;
	}
	return 1 - r/ax;
}
452358
private void calcRadiation( int i, int j, WritableRaster demWR, WritableRaster sOmbraWR, WritableRaster insolationWR,
            double[] sunVector, WritableRaster gradientWR, double mr ) {
        double z = demWR.getSampleDouble(i, j, 0);
        double pressure = ATM * Math.exp(-0.0001184 * z);
        double ma = mr * pressure / ATM;
        double temp = 273 + pLapse * (z - 4000);
        double vap_psat = Math.exp(26.23 - 5416.0 / temp);
        double wPrec = 0.493 * pRH * vap_psat / temp;
        double taur = Math.exp((-.09030 * Math.pow(ma, 0.84)) * (1.0 + ma - Math.pow(ma, 1.01)));
        double d = pCmO3 * mr;
        double tauo = 1 - (0.1611 * d * Math.pow(1.0 + 139.48 * d, -0.3035) - 0.002715 * d)
                / (1.0 + 0.044 * d + 0.0003 * Math.pow(d, 2));
        double taug = Math.exp(-0.0127 * Math.pow(ma, 0.26));
        double tauw = 1 - 2.4959 * (wPrec * mr) / (1.0 + 79.034 * (wPrec * mr) * 0.6828 + 6.385 * (wPrec * mr));
        double taua = Math.pow((0.97 - 1.265 * Math.pow(pVisibility, (-0.66))), Math.pow(ma, 0.9));

        double In = 0.9751 * SOLARCTE * taur * tauo * taug * tauw * taua;

        double cosinc = scalarProduct(sunVector, gradientWR.getPixel(i, j, new double[3]));

        if (cosinc < 0) {
            cosinc = 0;
        }
        double tmp = insolationWR.getSampleDouble(i, j, 0);
        insolationWR.setSample(i, j, 0, In * cosinc * sOmbraWR.getSampleDouble(i, j, 0) / 1000 + tmp);
    }
951229
def _get_sa_at_1180(self, C, imt, sites, rup, dists):
        
        # reference vs30 = 1180 m/s
        vs30_1180 = np.ones_like(sites.vs30) * 1180.
        # reference shaking intensity = 0
        ref_iml = np.zeros_like(sites.vs30)
        # fake Z1.0 - Since negative it will be replaced by the default Z1.0
        # for the corresponding region
        fake_z1pt0 = np.ones_like(sites.vs30) * -1
        return (self._get_basic_term(C, rup, dists) +
                self._get_faulting_style_term(C, rup) +
                self._get_site_response_term(C, imt, vs30_1180, ref_iml) +
                self._get_hanging_wall_term(C, dists, rup) +
                self._get_top_of_rupture_depth_term(C, imt, rup) +
                self._get_soil_depth_term(C, fake_z1pt0, vs30_1180) +
                self._get_regional_term(C, imt, vs30_1180, dists.rrup)
                )
968685
def postSolve(self):
        '''
        
        '''
        self.switchFuncPunk = self.solution[0].switchFuncPunk
        self.switchFuncJock = self.solution[0].switchFuncJock
        self.VfuncPunk      = self.solution[0].VfuncPunk
        self.VfuncJock      = self.solution[0].VfuncJock
555844
def  dBinaryRochedz(r, D, q, F):
    
    return -r[2]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*r[2]*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5
1227524
private double hermiteInterpolate(final double Y0, final double Y1, final double Y2, final double Y3,
                                      final double MU, final double TENSION, final double BIAS) {
        double m0;
        double m1;
        final double MU2;
        final double Mu3;
        final double A0;
        final double A1;
        final double A2;
        final double A3;

        MU2 = MU * MU;
        Mu3 = MU2 * MU;
        m0 = (Y1 - Y0) * (1 + BIAS) * (1 - TENSION) / 2;
        m0 += (Y2 - Y1) * (1 - BIAS) * (1 - TENSION) / 2;
        m1 = (Y2 - Y1) * (1 + BIAS) * (1 - TENSION) / 2;
        m1 += (Y3 - Y2) * (1 - BIAS) * (1 - TENSION) / 2;
        A0 = (2 * Mu3) - (3 * MU2) + 1;
        A1 = Mu3 - (2 * MU2) + MU;
        A2 = Mu3 - MU2;
        A3 = (-2 * Mu3) + (3 * MU2);

        return ((A0 * Y1) + (A1 * m0) + (A2 * m1) + (A3 * Y2));
    }
381956
def _reciprocal(coords1, charges1, coords2, charges2, kmax, kappa, box):
    
    n = coords1.shape[0]
    m = coords2.shape[0]
    result = np.zeros(n, dtype=np.float64)
    need_self = np.zeros(n, dtype=np.uint8)

    # Reciprocal unit vectors
    g1, g2, g3 = reciprocal_vectors(box)
    
    V = box_volume(box)
    
    prefac = 1.0 / (np.pi * V)
    for i in range(n):
        q_i = charges1[i]
        r_i = coords1[i]
        
        for j in range(m):
            q_j = charges2[j]
            r_j = coords2[j]
            
            r_ij = _dist(r_i, r_j)
            if r_ij < 1e-10:
                need_self[i] = 1
        
            for k_i in range(-kmax, kmax + 1):
                for k_j in range(-kmax, kmax + 1):
                    for k_k in range(-kmax, kmax + 1):
                        if k_i == 0 and k_j == 0 and k_k == 0:
                            continue
                        # Reciprocal vector
                        k = k_i * g1 + k_j * g2 + k_k * g3
                        
                        k_sq = sqsum(k)
                        
                        result[i] += (prefac * q_i * q_j *
                                      4.0 * np.pi ** 2 / k_sq *
                                      math.exp(-k_sq / (4.0 * kappa ** 2)) *
                                      math.cos(np.dot(k, r_i - r_j)))
        
    # Self-energy correction
    # I had to do some FUCKED UP stuff because NUMBA SUCKS BALLS and 
    # breaks compatibility with simple expressions such as that one
    # apparently doing result -= something is too hard in numba
    self_energy = 2 * (need_self * kappa * charges1 ** 2) / (np.pi**0.5)
    return result - self_energy
1083599
def _Fgamma(self, x, Ep):
        
        L = np.log(Ep)
        B = 1.30 + 0.14 * L + 0.011 * L ** 2  # Eq59
        beta = (1.79 + 0.11 * L + 0.008 * L ** 2) ** -1  # Eq60
        k = (0.801 + 0.049 * L + 0.014 * L ** 2) ** -1  # Eq61
        xb = x ** beta

        F1 = B * (np.log(x) / x) * ((1 - xb) / (1 + k * xb * (1 - xb))) ** 4
        F2 = (
            1.0 / np.log(x)
            - (4 * beta * xb) / (1 - xb)
            - (4 * k * beta * xb * (1 - 2 * xb)) / (1 + k * xb * (1 - xb))
        )

        return F1 * F2
817142
def _tilt(omega, b, psi):
    
    return np.cosh(psi/2.0)**b * np.exp(-psi**2/2.0 * omega)
625274
function waveUpdate () {
    var row, col, ind
    var dum = conv2(u, coeffs)
     for (row = 0; row < height; ++row)
      for (col = 0; col < width; ++col) {
        ind = row * width + col
        un[ind] = c1 * u[ind] + c2 * up[ind] + c3 * dum[ind]
        up[ind] = u[ind] // current becomes old
        u[ind] = un[ind] // new becomes current
      }
    return u
  }
224047
def _get_hanging_wall_coeffs_mag(self, C, mag):
        
        if mag < 5.5:
            return 0.0
        elif mag > 6.5:
            return 1.0 + C["a2"] * (mag - 6.5)
        else:
            return (mag - 5.5) * (1.0 + C["a2"] * (mag - 6.5))
1220746
def hum1_wei(x,y,n=24):
    t = y-1.0j*x
    cerf=1/sqrt(pi)*t/(0.5+t**2)
    
    mask = abs(x)+y<15.0
    if any(mask):
        w24 = weideman(x[mask],y[mask],n)
        place(cerf,mask,w24)
    return cerf.real,cerf.imag
375536
private void computePsiSmooth(GrayF32 ux , GrayF32 uy , GrayF32 vx , GrayF32 vy ,
								  GrayF32 psiSmooth ) {
		int N = derivFlowUX.width * derivFlowUX.height;

		for( int i = 0; i < N; i++ ) {
			float vux = ux.data[i];
			float vuy = uy.data[i];
			float vvx = vx.data[i];
			float vvy = vy.data[i];

			float mu = vux*vux + vuy*vuy;
			float mv = vvx*vvx + vvy*vvy;

			psiSmooth.data[i] = (float)(1.0/(2.0*Math.sqrt(mu + mv + EPSILON*EPSILON)));
		}
	}
74701
function forward$12(p) {

        /* Forward equations
         -----------------*/
        var x, y, coslam, sinlam, sinphi, q, sinb, cosb, b, cosphi;
        var lam = p.x;
        var phi = p.y;

        lam = adjust_lon(lam - this.long0);
        if (this.sphere) {
            sinphi = Math.sin(phi);
            cosphi = Math.cos(phi);
            coslam = Math.cos(lam);
            if (this.mode === this.OBLIQ || this.mode === this.EQUIT) {
                y = (this.mode === this.EQUIT) ? 1 + cosphi * coslam : 1 + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam;
                if (y <= EPSLN) {
                    return null;
                }
                y = Math.sqrt(2 / y);
                x = y * cosphi * Math.sin(lam);
                y *= (this.mode === this.EQUIT) ? sinphi : this.cosph0 * sinphi - this.sinph0 * cosphi * coslam;
            }
            else if (this.mode === this.N_POLE || this.mode === this.S_POLE) {
                if (this.mode === this.N_POLE) {
                    coslam = -coslam;
                }
                if (Math.abs(phi + this.phi0) < EPSLN) {
                    return null;
                }
                y = FORTPI - phi * 0.5;
                y = 2 * ((this.mode === this.S_POLE) ? Math.cos(y) : Math.sin(y));
                x = y * Math.sin(lam);
                y *= coslam;
            }
        }
        else {
            sinb = 0;
            cosb = 0;
            b = 0;
            coslam = Math.cos(lam);
            sinlam = Math.sin(lam);
            sinphi = Math.sin(phi);
            q = qsfnz(this.e, sinphi);
            if (this.mode === this.OBLIQ || this.mode === this.EQUIT) {
                sinb = q / this.qp;
                cosb = Math.sqrt(1 - sinb * sinb);
            }
            switch (this.mode) {
                case this.OBLIQ:
                    b = 1 + this.sinb1 * sinb + this.cosb1 * cosb * coslam;
                    break;
                case this.EQUIT:
                    b = 1 + cosb * coslam;
                    break;
                case this.N_POLE:
                    b = HALF_PI + phi;
                    q = this.qp - q;
                    break;
                case this.S_POLE:
                    b = phi - HALF_PI;
                    q = this.qp + q;
                    break;
            }
            if (Math.abs(b) < EPSLN) {
                return null;
            }
            switch (this.mode) {
                case this.OBLIQ:
                case this.EQUIT:
                    b = Math.sqrt(2 / b);
                    if (this.mode === this.OBLIQ) {
                        y = this.ymf * b * (this.cosb1 * sinb - this.sinb1 * cosb * coslam);
                    }
                    else {
                        y = (b = Math.sqrt(2 / (1 + cosb * coslam))) * sinb * this.ymf;
                    }
                    x = this.xmf * b * cosb * sinlam;
                    break;
                case this.N_POLE:
                case this.S_POLE:
                    if (q >= 0) {
                        x = (b = Math.sqrt(q)) * sinlam;
                        y = coslam * ((this.mode === this.S_POLE) ? b : -b);
                    }
                    else {
                        x = y = 0;
                    }
                    break;
            }
        }

        p.x = this.a * x + this.x0;
        p.y = this.a * y + this.y0;
        return p;
    }