Conversion from Mean Anomaly to Eccentric anomaly

        Procedures for solving Kepler's Equation, A. W. Odell and  R. H. Gooding,
        Celestial Mechanics 38 (1986) 307-334
424756
def _m_to_e(cls, e, M):
        

        k1 = 3 * np.pi + 2
        k2 = np.pi - 1
        k3 = 6 * np.pi - 1
        A = 3 * k2 ** 2 / k1
        B = k3 ** 2 / (6 * k1)

        m1 = float(M)
        if abs(m1) < 1 / 6:
            E = m1 + e * (6 * m1) ** (1 / 3) - m1
        elif m1 < 0:
            w = np.pi + m1
            E = m1 + e * (A * w / (B - w) - np.pi - m1)
        else:
            w = np.pi - m1
            E = m1 + e * (np.pi - A * w / (B - w) - m1)

        e1 = 1 - e
        risk_disabler = (e1 + E ** 2 / 6) >= 0.1

        for i in range(2):
            fdd = e * sin(E)
            fddd = e * cos(E)

            if risk_disabler:
                f = (E - fdd) - m1
                fd = 1 - fddd
            else:
                f = cls._e_e_sin_e(e, E) - m1
                s = sin(E / 2)
                fd = e1 + 2 * e * s ** 2
            dee = f * fd / (0.5 * f * fdd - fd ** 2)

            w = fd + 0.5 * dee * (fdd + dee * fddd / 3)
            fd += dee * (fdd + 0.5 * dee * fddd)
            E -= (f - dee * (fd - w)) / fd

        E += M - m1

        return E
376029
def _tle_to_keplerian_mean(cls, coord, center):
        
        i, Ω, e, ω, M, n = coord
        a = (center.µ / n ** 2) ** (1 / 3)

        return np.array([a, e, i, Ω, ω, M], dtype=float)
1131772
def ph_supconj(b, orbit, solve_for=None, **kwargs):
    
    orbit_ps = _get_system_ps(b, orbit)

    # metawargs = orbit_ps.meta
    #metawargs.pop('qualifier')

    # t0_ph0 and phshift both exist by default, so we don't have to worry about creating either
    # t0_ph0 = orbit_ps.get_parameter(qualifier='t0_ph0')
    # phshift = orbit_ps.get_parameter(qualifier='phshift')
    ph_supconj = orbit_ps.get_parameter(qualifier='ph_supconj')
    per0 = orbit_ps.get_parameter(qualifier='per0')
    ecc = orbit_ps.get_parameter(qualifier='ecc')
    period = orbit_ps.get_parameter(qualifier='period')

    # true_anom_supconj = pi/2 - per0
    # mean_anom_supconj = true_anom_supconj - ecc*sin(true_anom_supconj)
    # ph_supconj = (mean_anom_supconj + per0) / (2 * pi) - 1/4

    if solve_for in [None, ph_supconj]:
        lhs = ph_supconj

        #true_anom_supconj = np.pi/2*u.rad - per0
        true_anom_supconj = -1*(per0 - 360*u.deg)

        rhs = _true_anom_to_phase(true_anom_supconj, period, ecc, per0)

    #elif solve_for in [per0]:
    #    raise NotImplementedError("phshift constraint does not support solving for per0 yet")
    else:
        raise NotImplementedError

    return lhs, rhs, {'orbit': orbit}
1353746
private float[] restoreNormals(int[] intNormals, float[] vertices, int[] indices, float normalPrecision) {

        // Calculate smooth normals (nominal normals)
        float[] smoothNormals = calcSmoothNormals(vertices, indices);
        float[] normals = new float[vertices.length];

        int vc = vertices.length / CTM_POSITION_ELEMENT_COUNT;
        int ne = CTM_NORMAL_ELEMENT_COUNT;

        for (int i = 0; i < vc; ++i) {
            // Get the normal magnitude from the first of the three normal elements
            float magn = intNormals[i * ne] * normalPrecision;

            // Get phi and theta (spherical coordinates, relative to the smooth normal).
            double thetaScale, theta;
            int intPhi = intNormals[i * ne + 1];
            double phi = intPhi * (0.5 * PI) * normalPrecision;
            if (intPhi == 0) {
                thetaScale = 0.0f;
            } else if (intPhi <= 4) {
                thetaScale = PI / 2.0f;
            } else {
                thetaScale = (2.0f * PI) / intPhi;
            }
            theta = intNormals[i * ne + 2] * thetaScale - PI;

            // Convert the normal from the angular representation (phi, theta) back to
            // cartesian coordinates
            double[] n2 = new double[3];
            n2[0] = sin(phi) * cos(theta);
            n2[1] = sin(phi) * sin(theta);
            n2[2] = cos(phi);
            float[] basisAxes = makeNormalCoordSys(smoothNormals, i * ne);
            double[] n = new double[3];
            for (int j = 0; j < 3; ++j) {
                n[j] = basisAxes[j] * n2[0]
                       + basisAxes[3 + j] * n2[1]
                       + basisAxes[6 + j] * n2[2];
            }

            // Apply normal magnitude, and output to the normals array
            for (int j = 0; j < 3; ++j) {
                normals[i * ne + j] = (float) (n[j] * magn);
            }
        }

        return normals;
    }
376163
def yield_pair_hessians(self, index1, index2):
        """"""
        d_1 = 1/self.distances[index1, index2]
        d_3 = d_1**3
        if self.charges is not None:
            c1 = self.charges[index1]
            c2 = self.charges[index2]
            yield 2*c1*c2*d_3, np.zeros((3, 3))
        if self.dipoles is not None:
            d_5 = d_1**5
            d_7 = d_1**7
            p1 = self.dipoles[index1]
            p2 = self.dipoles[index2]
            yield 12*d_5*np.dot(p1, p2), np.zeros((3, 3))
            yield -90*d_7, np.outer(p1, p2) + np.outer(p2, p1)
            if self.charges is not None:
                yield 12*c1*d_5, np.zeros((3, 3))
                yield 12*c2*d_5, np.zeros((3, 3))
329950
def point
      unless @point
        # variable names correspond roughly to symbols in the OS algorithm, lowercased:
        # n0 = northing of true origin 
        # e0 = easting of true origin 
        # f0 = scale factor on central meridian
        # phi0 = latitude of true origin 
        # lambda0 = longitude of true origin and central meridian.
        # e2 = eccentricity squared
        # a = length of polar axis of ellipsoid
        # b = length of equatorial axis of ellipsoid
        # ning & eing are the northings and eastings of the supplied gridref
        # phi and lambda are the discovered latitude and longitude
      
        ning = northing
        eing = easting

        n0 = projection.n0
        e0 = projection.e0
        phi0 = projection.phi0
        l0 = projection.lambda0
        f0 = projection.scale
      
        a = ellipsoid.a
        b = ellipsoid.b
        e2 = ellipsoid.ecc
      
        # the rest is juste a transliteration of the OS equations

        n = (a - b) / (a + b)
        m = 0
        phi = phi0
    
        # iterate to within acceptable distance of solution
      
        count = 0
        while ((ning - n0 - m) >= 0.001) do
          raise RuntimeError "Demercatorising equation has not converged. Discrepancy after #{count} cycles is #{ning - n0 - m}" if count >= @@iteration_ceiling

          phi = ((ning - n0 - m) / (a * f0)) + phi
          ma = (1 + n + (1.25 * n**2) + (1.25 * n**3)) * (phi - phi0)
          mb = ((3 * n) + (3 * n**2) + (2.625 * n**3)) * Math.sin(phi - phi0) * Math.cos(phi + phi0)
          mc = ((1.875 * n**2) + (1.875 * n**3)) * Math.sin(2 * (phi - phi0)) * Math.cos(2 * (phi + phi0))
          md = (35/24) * (n**3) * Math.sin(3 * (phi - phi0)) * Math.cos(3 * (phi + phi0))
          m = b * f0 * (ma - mb + mc - md)
          count += 1
        end
      
        # engage alphabet soup
      
        nu = a * f0 * ((1-(e2) * ((Math.sin(phi)**2))) ** -0.5)
        rho = a * f0 * (1-(e2)) * ((1-(e2)*((Math.sin(phi)**2))) ** -1.5)
        eta2 = (nu/rho - 1)
      
        # fire
      
        vii = Math.tan(phi) / (2 * rho * nu)
        viii = (Math.tan(phi) / (24 * rho * (nu ** 3))) * (5 + (3 * (Math.tan(phi) ** 2)) + eta2 - 9 * eta2 * (Math.tan(phi) ** 2) )
        ix = (Math.tan(phi) / (720 * rho * (nu ** 5))) * (61 + (90 * (Math.tan(phi) ** 2)) + (45 * (Math.tan(phi) ** 4)) )
        x = sec(phi) / nu
        xi = (sec(phi) / (6 * nu ** 3)) * ((nu/rho) + 2 * (Math.tan(phi) ** 2))
        xii = (sec(phi) / (120 * nu ** 5)) * (5 + (28 * (Math.tan(phi) ** 2)) + (24 * (Math.tan(phi) ** 4)))
        xiia = (sec(phi) / (5040 * nu ** 7)) * (61 + (662 * (Math.tan(phi) ** 2)) + (1320 * (Math.tan(phi) ** 4)) + (720 * (Math.tan(phi) ** 6)))

        d = eing-e0

        # all of which was just to populate these last two equations:
      
        phi = phi - vii*(d**2) + viii*(d**4) - ix*(d**6)
        lambda = l0 + x*d - xi*(d**3) + xii*(d**5) - xiia*(d**7)

        @point = Osgb::Point.new(phi.to_degrees, lambda.to_degrees, :osgb36, precision)
      end
      @point
    end
771400
def _minimize_c(c, z=0, a_tilde=1, b_tilde=-1,
                Ascaling=900, omega_M_0=0.25, omega_lambda_0=0.75):
    

    # Fn 1 (LHS of Eqn 18)

    Y1 = np.log(2) - 0.5
    Yc = np.log(1+c) - c/(1+c)
    f1 = Y1/Yc

    # Fn 2 (RHS of Eqn 18)

    # Eqn 14 - Define the mean inner density
    rho_2 = 200 * c**3 * Y1 / Yc

    # Eqn 17 rearranged to solve for Formation Redshift
    # essentially when universe had rho_2 density
    zf = (((1 + z)**3 + omega_lambda_0/omega_M_0) *
          (rho_2/Ascaling) - omega_lambda_0/omega_M_0)**(1/3) - 1

    # RHS of Eqn 19
    f2 = ((1 + zf - z)**a_tilde) * np.exp((zf - z) * b_tilde)

    # LHS - RHS should be zero for the correct concentration
    return(f1-f2)
1014561
private static function sinDeltaBearing( $lat1, $lng1, $lat2, $lng2, $lat3, $lng3) {
        
        $sinLat1 = sin($lat1);
        $cosLat2 = cos($lat2);
        $cosLat3 = cos($lat3);
        $lat31 = $lat3 - $lat1;
        $lng31 = $lng3 - $lng1;
        $lat21 = $lat2 - $lat1;
        $lng21 = $lng2 - $lng1;
        $a = sin($lng31) * $cosLat3;
        $c = sin($lng21) * $cosLat2;
        $b = sin($lat31) + 2 * $sinLat1 * $cosLat3 * MathUtil::hav($lng31);
        $d = sin($lat21) + 2 * $sinLat1 * $cosLat2 * MathUtil::hav($lng21);
        $denom = ($a * $a + $b * $b) * ($c * $c + $d * $d);
        return $denom <= 0 ? 1 : ($a * $d - $b * $c) / sqrt($denom);
    }
596623
def eclipseprob(self):
        
        #TODO: incorporate eccentricity/omega for exact calculation?
        s = self.stars
        return ((s['radius_1'] + s['radius_2'])*RSUN /
                (semimajor(s['P'],s['mass_1'] + s['mass_2'])*AU))
862692
def spherical_galaxy_orbit(
    orbit_x,
    orbit_y,
    orbit_z,
    N_stars=100,
    sigma_r=1,
    orbit_visible=False,
    orbit_line_interpolate=5,
    N_star_orbits=10,
    color=[255, 220, 200],
    size_star=1,
    scatter_kwargs={},
):
    """"""
    if orbit_line_interpolate > 1:
        x = np.linspace(0, 1, len(orbit_x))
        x_smooth = np.linspace(0, 1, len(orbit_x) * orbit_line_interpolate)
        kind = 'quadratic'
        orbit_x_line = scipy.interpolate.interp1d(x, orbit_x, kind)(x_smooth)
        orbit_y_line = scipy.interpolate.interp1d(x, orbit_y, kind)(x_smooth)
        orbit_z_line = scipy.interpolate.interp1d(x, orbit_z, kind)(x_smooth)
    else:
        orbit_x_line = orbit_x
        orbit_y_line = orbit_y
        orbit_z_line = orbit_z
    line = ipv.plot(orbit_x_line, orbit_y_line, orbit_z_line, visible=orbit_visible)
    x = np.repeat(orbit_x, N_stars).reshape((-1, N_stars))
    y = np.repeat(orbit_y, N_stars).reshape((-1, N_stars))
    z = np.repeat(orbit_z, N_stars).reshape((-1, N_stars))
    xr, yr, zr = np.random.normal(0, scale=sigma_r, size=(3, N_stars))  # +
    r = np.sqrt(xr ** 2 + yr ** 2 + zr ** 2)

    for i in range(N_stars):
        a = np.linspace(0, 1, x.shape[0]) * 2 * np.pi * N_star_orbits
        xo = r[i] * np.sin(a)
        yo = r[i] * np.cos(a)
        zo = a * 0
        xo, yo, zo = np.dot(_randomSO3(), [xo, yo, zo])
        # print(x.shape, xo.shape)
        x[:, i] += xo
        y[:, i] += yo
        z[:, i] += zo

    sprite = ipv.scatter(
        x, y, z, texture=radial_sprite((64, 64), color), marker='square_2d', size=size_star, **scatter_kwargs
    )
    with sprite.material.hold_sync():
        sprite.material.blending = pythreejs.BlendingMode.CustomBlending  # pylint: disable=no-member
        sprite.material.blendSrc = pythreejs.BlendFactors.SrcColorFactor  # pylint: disable=no-member
        sprite.material.blendDst = pythreejs.BlendFactors.OneFactor  # pylint: disable=no-member
        sprite.material.blendEquation = 'AddEquation'
        sprite.material.transparent = True
        sprite.material.depthWrite = False
        sprite.material.alphaTest = 0.1
    return sprite, line
1263626
private void orthes(Matrix h, Matrix v, Vector ort) {
        // This is derived from the Algol procedures orthes and ortran,
        // by Martin and Wilkinson, Handbook for Auto. Comp.,
        // Vol.ii-Linear Algebra, and the corresponding
        // Fortran subroutines in EISPACK.

        int n = ort.length();
        int low = 0;
        int high = n - 1;

        for (int m = low + 1; m <= high - 1; m++) {

            // Scale column.

            double scale = 0.0;

            for (int i = m; i <= high; i++) {
                scale = scale + Math.abs(h.get(i, m - 1));
            }

            if (scale != 0.0) {

                // Compute Householder transformation.

                double hh = 0.0;
                for (int i = high; i >= m; i--) {
                    ort.set(i, h.get(i, m - 1) / scale);
                    hh += ort.get(i) * ort.get(i);
                }

                double g = Math.sqrt(hh);

                if (ort.get(m) > Matrices.EPS) {
                    g = -g;
                }

                hh = hh - ort.get(m) * g;
                ort.updateAt(m, Vectors.asMinusFunction(g));

                // Apply Householder similarity transformation
                // H = (I-u*u'/h)*H*(I-u*u')/h)

                for (int j = m; j < n; j++) {
                    double f = 0.0;
                    for (int i = high; i >= m; i--) {
                        f += ort.get(i) * h.get(i, j);
                    }
                    f = f / hh;
                    for (int i = m; i <= high; i++) {
                        h.updateAt(i, j, Matrices.asMinusFunction(f * ort.get(i)));
                    }
                }

                for (int i = 0; i <= high; i++) {
                    double f = 0.0;
                    for (int j = high; j >= m; j--) {
                        f += ort.get(j) * h.get(i, j);
                    }
                    f = f / hh;
                    for (int j = m; j <= high; j++) {
                        h.updateAt(i, j, Matrices.asMinusFunction(f * ort.get(j)));
                    }
                }
                ort.set(m, scale * ort.get(m));
                h.set(m,  m - 1, scale * g);
            }
        }

        // Accumulate transformations (Algol's ortran).

        for (int m = high - 1; m >= low + 1; m--) {
            if (Math.abs(h.get(m, m - 1)) > Matrices.EPS) {
                for (int i = m + 1; i <= high; i++) {
                    ort.set(i, h.get(i, m - 1));
                }
                for (int j = m; j <= high; j++) {
                    double g = 0.0;
                    for (int i = m; i <= high; i++) {
                        g += ort.get(i) * v.get(i, j);
                    }
                    // Double division avoids possible underflow
                    g = (g / ort.get(m)) / h.get(m, m - 1);
                    for (int i = m; i <= high; i++) {
                        v.updateAt(i, j, Matrices.asPlusFunction(g * ort.get(i)));
                    }
                }
            }
        }
    }
509319
def accel_ES(q: np.ndarray):
    

    # Number of celestial bodies
    num_bodies: int = 2
    # Number of dimensions in arrays; 3 spatial dimensions times the number of bodies
    dims = 3 * num_bodies

    # Body 0 is the sun; Body 1 is the earth
    m0 = mass[0]
    m1 = mass[1]

    # Extract position of the sun and earth as 3-vectors
    pos_0 = q[slices[0]]
    pos_1 = q[slices[1]]

    # Displacement vector from sun to earth
    dv_01: np.ndarray = pos_1 - pos_0

    # Distance from sun to earth
    r_01: float = np.linalg.norm(dv_01)
    
    # Unit vector pointing from sun to earth
    udv_01 = dv_01 / r_01

    # The force between these has magnitude G*m0*m1 / r^2
    f_01: float = (G * m0 * m1) / (r_01 ** 2)
    
    # Initialize acceleration as 6x1 array
    a: np.ndarray = np.zeros(dims)
    
    # The force vectors are attractive
    a[slices[0]] += f_01 * udv_01 / m0
    a[slices[1]] -= f_01 * udv_01 / m1

    # Return the acceleration vector
    return a
792102
def inclination(vector)
      Math.atan2(-vector.y, Math.sqrt(vector.x * vector.x + vector.z * vector.z))
    end
2344
def solver(AA, N_max, symNx = 2, throw_out_modes=False):
    

    # Find all integer component n_vectors which lie within sphere of radius N_max
    # Here we have assumed that the potential is symmetric x->-x, y->-y, z->-z
    # This can be relaxed by changing symN to 1
    # Additionally due to time reversal symmetry S_n = -S_-n so we only consider
    # "half" of the n-vector-space

    angs = unroll_angles(AA.T[3:].T,np.ones(3))

    symNz = 2
    NNx = range(-N_max, N_max+1, symNx)
    NNy = range(-N_max, N_max+1, symNz)
    NNz = range(-N_max, N_max+1, symNz)
    n_vectors = np.array([[i,j,k] for (i,j,k) in product(NNx,NNy,NNz)
                          if(not(i==0 and j==0 and k==0)            # exclude zero vector
                             and (k>0                               # northern hemisphere
                                  or (k==0 and j>0)                 # half of x-y plane
                                  or (k==0 and j==0 and i>0))       # half of x axis
                             and np.sqrt(i*i+j*j+k*k)<=N_max)])     # inside sphere

    xxx = check_each_direction(n_vectors,angs)

    if(throw_out_modes):
        n_vectors = np.delete(n_vectors,check_each_direction(n_vectors,angs),axis=0)

    n = len(n_vectors)+3
    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n,n))

    a[:3,:3]=len(AA)*np.identity(3)

    for i in AA:
        a[:3,3:]+=2.*n_vectors.T[:3]*np.cos(np.dot(n_vectors,i[3:]))
        a[3:,3:]+=4.*np.dot(n_vectors,n_vectors.T)*np.outer(np.cos(np.dot(n_vectors,i[3:])),np.cos(np.dot(n_vectors,i[3:])))
        b[:3]+=i[:3]
        b[3:]+=2.*np.dot(n_vectors,i[:3])*np.cos(np.dot(n_vectors,i[3:]))

    a[3:,:3]=a[:3,3:].T

    return np.array(solve(a,b)), n_vectors
771290
function forward$21(p) {

        /* Forward equations
         -----------------*/
        var lon = p.x;
        var lat = p.y;

        var delta_lon = adjust_lon(lon - this.long0);
        var theta = lat;
        var con = Math.PI * Math.sin(lat);

        /* Iterate using the Newton-Raphson method to find theta
         -----------------------------------------------------*/
        for (var i = 0; true; i++) {
            var delta_theta = -(theta + Math.sin(theta) - con) / (1 + Math.cos(theta));
            theta += delta_theta;
            if (Math.abs(delta_theta) < EPSLN) {
                break;
            }
        }
        theta /= 2;

        /* If the latitude is 90 deg, force the x coordinate to be "0 + false easting"
         this is done here because of precision problems with "cos(theta)"
         --------------------------------------------------------------------------*/
        if (Math.PI / 2 - Math.abs(lat) < EPSLN) {
            delta_lon = 0;
        }
        var x = 0.900316316158 * this.a * delta_lon * Math.cos(theta) + this.x0;
        var y = 1.4142135623731 * this.a * Math.sin(theta) + this.y0;

        p.x = x;
        p.y = y;
        return p;
    }
922125
function toLatLng()
	{
		$wgs84 = new ReferenceEllipsoid(ReferenceEllipsoid::WGS84_MAJ, ReferenceEllipsoid::WGS84_MIN);

		$UTM_F0	 = 0.9996;
		$a = $wgs84->maj;
		$eSquared = $wgs84->ecc;
		$ePrimeSquared = $eSquared / (1.0 - $eSquared);
		$e1 = (1 - sqrt(1 - $eSquared)) / (1 + sqrt(1 - $eSquared));

		$x = $this->easting - 500000.0;;
		$y = $this->northing;

		$zoneNumber = $this->lngZone;
		$zoneLetter = $this->latZone;

		$longitudeOrigin = ($zoneNumber - 1.0) * 6.0 - 180.0 + 3.0;

		// Correct y for southern hemisphere
		if ((ord($zoneLetter) - ord("N")) < 0) {
			$y -= 10000000.0;
		}

		$m = $y / $UTM_F0;

		$mu = $m
			/ ($a
			* (1.0
				- $eSquared / 4.0
				- 3.0 * $eSquared * $eSquared / 64.0
				- 5.0
				* pow($eSquared, 3.0)
				/ 256.0));

		$phi1Rad = $mu
			+ (3.0 * $e1 / 2.0 - 27.0 * pow($e1, 3.0) / 32.0) * sin(2.0 * $mu)
			+ (21.0 * $e1 * $e1 / 16.0 - 55.0 * pow($e1, 4.0) / 32.0)
			* sin(4.0 * $mu)
			+ (151.0 * pow($e1, 3.0) / 96.0) * sin(6.0 * $mu);

		$n = $a / sqrt(1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad));

		$t = tan($phi1Rad) * tan($phi1Rad);
		$c = $ePrimeSquared * cos($phi1Rad) * cos($phi1Rad);

		$r = $a
			* (1.0 - $eSquared)
			/ pow(
			1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad),
			1.5);

		$d = $x / ($n * $UTM_F0);

		$latitude = ($phi1Rad
			- ($n * tan($phi1Rad) / $r)
			* ($d * $d / 2.0
				- (5.0
				+ (3.0 * $t)
				+ (10.0 * $c)
				- (4.0 * $c * $c)
				- (9.0 * $ePrimeSquared))
				* pow($d, 4.0)
				/ 24.0
				+ (61.0
				+ (90.0 * $t)
				+ (298.0 * $c)
				+ (45.0 * $t * $t)
				- (252.0 * $ePrimeSquared)
				- (3.0 * $c * $c))
				* pow($d, 6.0)
				/ 720.0)) * (180.0 / pi());

		$longitude = $longitudeOrigin + (($d
			- (1.0 + 2.0 * $t + $c) * pow($d, 3.0) / 6.0
			+ (5.0
			- (2.0 * $c)
			+ (28.0 * $t)
			- (3.0 * $c * $c)
			+ (8.0 * $ePrimeSquared)
			+ (24.0 * $t * $t))
			* pow($d, 5.0)
			/ 120.0)
			/ cos($phi1Rad)) * (180.0 / pi());

		return new LatLng($latitude, $longitude);
	}
473475
def expected_magz(RAW_IMU, ATTITUDE, inclination, declination):
    ''''''
    v = expected_mag(RAW_IMU, ATTITUDE, inclination, declination)
    return v.z