Given two ellipses separated by *dx* and *dy*, compute their separation in
    terms of σ. Based on Pineau et al (2011A&A...527A.126P).

    The "0" ellipse is taken to be centered at (0, 0), while the "1"
    ellipse is centered at (dx, dy).
188945
def double_ell_distance (mjr0, mnr0, pa0, mjr1, mnr1, pa1, dx, dy):
    
    # 1. We need to rotate the frame so that ellipse 1 lies on the X axis.
    theta = -np.arctan2 (dy, dx)

    # 2. We also need to express these rotated ellipses in "biv" format.
    sx0, sy0, cxy0 = ellbiv (mjr0, mnr0, pa0 + theta)
    sx1, sy1, cxy1 = ellbiv (mjr1, mnr1, pa1 + theta)

    # 3. Their convolution is:
    sx, sy, cxy = bivconvolve (sx0, sy0, cxy0, sx1, sy1, cxy1)

    # 4. The separation between the centers is still just:
    d = np.sqrt (dx**2 + dy**2)

    # 5. The effective sigma in the purely X direction, taking into account
    # the covariance term, is:
    sigma_eff = sx * np.sqrt (1 - (cxy / (sx * sy))**2)

    # 6. Therefore the answer is:
    return d / sigma_eff
1111789
def ellipse2pathd(ellipse):
    """"""

    cx = ellipse.get('cx', 0)
    cy = ellipse.get('cy', 0)
    rx = ellipse.get('rx', None)
    ry = ellipse.get('ry', None)
    r = ellipse.get('r', None)

    if r is not None:
        rx = ry = float(r)
    else:
        rx = float(rx)
        ry = float(ry)

    cx = float(cx)
    cy = float(cy)

    d = ''
    d += 'M' + str(cx - rx) + ',' + str(cy)
    d += 'a' + str(rx) + ',' + str(ry) + ' 0 1,0 ' + str(2 * rx) + ',0'
    d += 'a' + str(rx) + ',' + str(ry) + ' 0 1,0 ' + str(-2 * rx) + ',0'

    return d
490566
def XX(func, a, b, c, d, x, s, ac):
    

    res = 0L
    res = res + a + func(b, c, d)
    res = res + x 
    res = res + ac
    res = res & 0xffffffffL
    res = _rotateLeft(res, s)
    res = res & 0xffffffffL
    res = res + b

    return res & 0xffffffffL
523897
function fe(Ce,Ne){var ze=Ce.c,De=Ne.c,Me=Ce.s,Ue=Ne.s,Fe=Ce.e,He=Ne.e,je,Le;// Either NaN?
if(!Me||!Ue)return null;// Either zero?
if(je=ze&&!ze[0],Le=De&&!De[0],je||Le)return je?Le?0:-Ue:Me;// Signs differ?
if(Me!=Ue)return Me;// Either Infinity?
if(je=0>Me,Le=Fe==He,!ze||!De)return Le?0:!ze^je?1:-1;// Compare exponents.
if(!Le)return Fe>He^je?1:-1;// Compare digit by digit.
for(Ue=(Fe=ze.length)<(He=De.length)?Fe:He,Me=0;Me<Ue;Me++)if(ze[Me]!=De[Me])return ze[Me]>De[Me]^je?1:-1;// Compare lengths.
return Fe==He?0:Fe>He^je?1:-1}
1210070
public void ellipse(float x1, float y1, float x2, float y2) {
        arc(x1, y1, x2, y2, 0f, 360f);
    }
126866
function _EllipseArc($x0, $y0, $rx, $ry, $seg = 1, $part = false, $start = false)
	{
		// Anticlockwise segment 1-4 TR-TL-BL-BR (part=1 or 2)
		$s = '';

		if ($rx < 0) {
			$rx = 0;
		}

		if ($ry < 0) {
			$ry = 0;
		}

		$rx *= Mpdf::SCALE;
		$ry *= Mpdf::SCALE;

		$astart = 0;

		if ($seg == 1) { // Top Right
			$afinish = 90;
			$nSeg = 4;
		} elseif ($seg == 2) { // Top Left
			$afinish = 180;
			$nSeg = 8;
		} elseif ($seg == 3) { // Bottom Left
			$afinish = 270;
			$nSeg = 12;
		} else {   // Bottom Right
			$afinish = 360;
			$nSeg = 16;
		}

		$astart = deg2rad((float) $astart);
		$afinish = deg2rad((float) $afinish);

		$totalAngle = $afinish - $astart;
		$dt = $totalAngle / $nSeg; // segment angle
		$dtm = $dt / 3;
		$x0 *= Mpdf::SCALE;
		$y0 = ($this->h - $y0) * Mpdf::SCALE;
		$t1 = $astart;
		$a0 = $x0 + ($rx * cos($t1));
		$b0 = $y0 + ($ry * sin($t1));
		$c0 = -$rx * sin($t1);
		$d0 = $ry * cos($t1);
		$op = false;

		for ($i = 1; $i <= $nSeg; $i++) {
			// Draw this bit of the total curve
			$t1 = ($i * $dt) + $astart;
			$a1 = $x0 + ($rx * cos($t1));
			$b1 = $y0 + ($ry * sin($t1));
			$c1 = -$rx * sin($t1);
			$d1 = $ry * cos($t1);
			if ($i > ($nSeg - 4) && (!$part || ($part == 1 && $i <= $nSeg - 2) || ($part == 2 && $i > $nSeg - 2))) {
				if ($start && !$op) {
					$s .= sprintf('%.3F %.3F m ', $a0, $b0);
				}
				$s .= sprintf('%.3F %.3F %.3F %.3F %.3F %.3F c ', ($a0 + ($c0 * $dtm)), ($b0 + ($d0 * $dtm)), ($a1 - ($c1 * $dtm)), ($b1 - ($d1 * $dtm)), $a1, $b1);
				$op = true;
			}
			$a0 = $a1;
			$b0 = $b1;
			$c0 = $c1;
			$d0 = $d1;
		}

		return $s;
	}
411345
private int locate(double x) {
        int ju, jm, jl;

        boolean ascnd = (xx[n - 1] >= xx[0]);

        jl = 0;
        ju = n - 1;
        while (ju - jl > 1) {
            jm = (ju + jl) >> 1;
            if (x >= xx[jm] == ascnd) {
                jl = jm;
            } else {
                ju = jm;
            }
        }

        cor = Math.abs(jl - jsav) <= dj;
        jsav = jl;

        return Math.max(0, Math.min(n - 2, jl));
    }
298931
function (a, b) {
            if (!b.pos || !a.pos) {
                return (a.pos ? -Infinity : Infinity);
            }
            var result = b.pos.z - a.pos.z;
            return (result ? result : (b.pos.x - a.pos.x));
        }
131906
def center_of_mass(X, Y):
    
    X = np.asarray(X)
    Y = np.asarray(Y)
    return (X * Y).sum() / Y.sum()
247168
def __reverse_ellipse(axes, scalar=1):
    
    ax1 = axes.copy()[::-1]*scalar
    center = ax1[1]*N.sqrt(2)*scalar
    return ax1, center
304942
function eq (y, c) {
  const J0 = base.horner(y * 0.001, c)
  const T = base.J2000Century(J0)
  const W = 35999.373 * D2R * T - 2.47 * D2R
  const Δλ = 1 + 0.0334 * cos(W) + 0.0007 * cos(2 * W)
  let S = 0
  for (let i = terms.length - 1; i >= 0; i--) {
    const t = terms[i]
    S += t.a * cos((t.b + t.c * T) * D2R)
  }
  return J0 + 0.00001 * S / Δλ
}
725577
def point_distance_ellipsode(point1,point2):
    
    a = 6378137
    f = 1/298.25722
    b = a - a*f
    e = math.sqrt((a*a-b*b)/(a*a))
    lon1 = point1['coordinates'][0]
    lat1 = point1['coordinates'][1]
    lon2 = point1['coordinates'][0]
    lat2 = point2['coordinates'][1]
    M = a*(1-e*e)*math.pow(1-math.pow(e*math.sin(number2radius(lat1)),2),-1.5)
    N = a/(math.pow(1-math.pow(e*math.sin(number2radius(lat1)),2),0.5))
    
    distance_lat = M*number2radius(lat2-lat1)
    distance_lon = N*math.cos(number2radius(lat1))*(lon2-lon1)*3600*math.sin(1/3600*math.pi/180)
    return math.sqrt(distance_lat*distance_lat+distance_lon*distance_lon)
105292
def transform_e1e2(x, y, e1, e2, center_x=0, center_y=0):
    
    x_shift = x - center_x
    y_shift = y - center_y
    x_ = (1-e1) * x_shift - e2 * y_shift
    y_ = -e2 * x_shift + (1 + e1) * y_shift
    det = np.sqrt((1-e1)*(1+e1) + e2**2)
    return x_ / det, y_ / det
39731
def b_neg(self):
        '''
        
        '''
        nb = Point(self.center)

        if self.xAxisIsMinor:
            nb.x -= self.minorRadius
        else:
            nb.y -= self.minorRadius
        return nb
675129
private void computeDerivatives(double sigma) {
		this.Lx2 = new float[width][height];
		this.Ly2 = new float[width][height];
		this.Lxy = new float[width][height];

		float[][][] grad = new float[width][height][];
		for (int y = 0; y < height; y++)
			for (int x = 0; x < width; x++)
				grad[x][y] = sobel(x, y);

		int radius = (int) (2 * sigma);
		int window = 1 + 2 * radius;
		float[][] gaussian = new float[window][window];
		for (int j = -radius; j <= radius; j++)
			for (int i = -radius; i <= radius; i++)
				gaussian[i + radius][j + radius] = (float) gaussian(i, j, sigma);

		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {

				for (int dy = -radius; dy <= radius; dy++) {
					for (int dx = -radius; dx <= radius; dx++) {
						int xk = x + dx;
						int yk = y + dy;
						if (xk < 0 || xk >= width)
							continue;
						if (yk < 0 || yk >= height)
							continue;

						double gw = gaussian[dx + radius][dy + radius];

						this.Lx2[x][y] += gw * grad[xk][yk][0]
								* grad[xk][yk][0];
						this.Ly2[x][y] += gw * grad[xk][yk][1]
								* grad[xk][yk][1];
						this.Lxy[x][y] += gw * grad[xk][yk][0]
								* grad[xk][yk][1];
					}
				}
			}
		}
	}
627725
function (p1, p2) {
    let direction = p2.minus(p1)
    let labda = (this.w - this.normal.dot(p1)) / this.normal.dot(direction)
    if (isNaN(labda)) labda = 0
    if (labda > 1) labda = 1
    if (labda < 0) labda = 0
    let result = p1.plus(direction.times(labda))
    return result
  }
1259562
def function(self, x, y, sigma0, Rs, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        if isinstance(r, int) or isinstance(r, float):
            r = max(self._s, r)
        else:
            r[r < self._s] = self._s
        X = r / Rs
        f_ = sigma0 * Rs ** 2 * (np.log(X ** 2 / 4.) + 2 * self._F(X))
        return f_