Python numpy.rad2deg() Examples

The following are 30 code examples for showing how to use numpy.rad2deg(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: ratcave   Author: ratcave   File: test_camera.py    License: MIT License 6 votes vote down vote up
def test_look_at_updates_for_children():
    dist = 2.
    cam = StereoCameraGroup(distance=dist)
    point = np.array([0, 0, 0, 1]).reshape(-1, 1)
    point[2] = -1 #np.random.randint(1, 6)

    angle = np.arctan(point[2]/(cam.distance/2))[0]
    cam.right.rotation.y = -np.rad2deg(angle)
    cam.left.rotation.y = np.rad2deg(angle)
    point_view_mat_left = np.dot(cam.left.view_matrix, point)
    point_view_mat_right = np.dot(cam.right.view_matrix, point)
    assert (point_view_mat_left == point_view_mat_right).all()

    cam2 = StereoCameraGroup(distance=dist)
    cam2.look_at(*point[:3])
    point_view_mat_left2 = np.dot(cam2.left.view_matrix, point)
    point_view_mat_right2 = np.dot(cam2.right.view_matrix, point)
    assert (point_view_mat_left == point_view_mat_left2).all() and (point_view_mat_right == point_view_mat_right2).all() 
Example 2
Project: pyeo   Author: clcr   File: test3.py    License: GNU General Public License v3.0 6 votes vote down vote up
def sample_data_3d(shape):
    """Returns `lons`, `lats`, and fake `data`

    adapted from:
    http://scitools.org.uk/cartopy/docs/v0.15/examples/axes_grid_basic.html
    """

    nlons, nlats = shape
    lats = np.linspace(-np.pi / 2, np.pi / 2, nlats)
    lons = np.linspace(0, 2 * np.pi, nlons)
    lons, lats = np.meshgrid(lons, lats)
    wave = 0.75 * (np.sin(2 * lats) ** 8) * np.cos(4 * lons)
    mean = 0.5 * np.cos(2 * lats) * ((np.sin(2 * lats)) ** 2 + 2)

    lats = np.rad2deg(lats)
    lons = np.rad2deg(lons)
    data = wave + mean

    return lons, lats, data


# get data 
Example 3
Project: visual_foresight   Author: SudeepDasari   File: get_points.py    License: MIT License 6 votes vote down vote up
def main_sawyer():
    import intera_interface
    from visual_mpc.envs.robot_envs.sawyer.sawyer_impedance import SawyerImpedanceController
    
    controller = SawyerImpedanceController('sawyer', True, gripper_attached='none')       # doesn't initial gripper object even if gripper is attached

    def print_eep(value):
        if not value:
            return
        xyz, quat = controller.get_xyz_quat()
        yaw, roll, pitch = [np.rad2deg(x) for x in controller.quat_2_euler(quat)]
        logging.getLogger('robot_logger').info("XYZ IS: {}, ROTATION IS: yaw={} roll={} pitch={}".format(xyz, yaw, roll, pitch))
    
    navigator = intera_interface.Navigator()
    navigator.register_callback(print_eep, 'right_button_show')
    rospy.spin() 
Example 4
Project: Mastering-Elasticsearch-7.0   Author: PacktPublishing   File: polar.py    License: MIT License 6 votes vote down vote up
def __call__(self, x, pos=None):
        vmin, vmax = self.axis.get_view_interval()
        d = np.rad2deg(abs(vmax - vmin))
        digits = max(-int(np.log10(d) - 1.5), 0)

        if rcParams['text.usetex'] and not rcParams['text.latex.unicode']:
            format_str = r"${value:0.{digits:d}f}^\circ$"
            return format_str.format(value=np.rad2deg(x), digits=digits)
        else:
            # we use unicode, rather than mathtext with \circ, so
            # that it will work correctly with any arbitrary font
            # (assuming it has a degree sign), whereas $5\circ$
            # will only work correctly with one of the supported
            # math fonts (Computer Modern and STIX)
            format_str = "{value:0.{digits:d}f}\N{DEGREE SIGN}"
            return format_str.format(value=np.rad2deg(x), digits=digits) 
Example 5
Project: BIRL   Author: Borda   File: registration.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def get_affine_components(matrix):
    """ get the main components of Affine transform

    :param ndarray matrix: affine transformation matrix for 2d
    :return dict: dictionary of float values

    >>> mtx = np.array([[ -0.95,   0.1,  65.], [  0.1,   0.95, -60.], [  0.,   0.,   1.]])
    >>> import  pandas as pd
    >>> aff = pd.Series(get_affine_components(mtx)).sort_index()
    >>> aff  # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
    rotation                           173.9...
    scale                    (0.95..., 0.95...)
    shear                              -3.14...
    translation                   (65.0, -60.0)
    dtype: object
    """
    aff = AffineTransform(matrix)
    norm_rotation = norm_angle(np.rad2deg(aff.rotation), deg=True)
    comp = {
        'rotation': float(norm_rotation),
        'translation': tuple(aff.translation.tolist()),
        'scale': aff.scale,
        'shear': aff.shear,
    }
    return comp 
Example 6
Project: kite   Author: pyrocko   File: base.py    License: GNU General Public License v3.0 6 votes vote down vote up
def orientArrow(self):
        phi = num.nanmedian(self.model.scene.phi)
        theta = num.nanmedian(self.model.scene.theta)

        angle = 180. - num.rad2deg(phi)

        theta_f = theta / (num.pi/2)

        tipAngle = 30. + theta_f * 20.
        tailLen = 15 + theta_f * 15.

        self.arrow.setStyle(
            angle=0.,
            tipAngle=tipAngle,
            tailLen=tailLen,
            tailWidth=6,
            headLen=25)
        self.arrow.setRotation(angle)

        rect_label = self.label.boundingRect()
        rect_arr = self.arrow.boundingRect()

        self.label.setPos(-rect_label.width()/2., rect_label.height()*1.33) 
Example 7
Project: pvlib-python   Author: pvlib   File: spa.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def heliocentric_latitude(jme):
    b0 = 0.0
    b1 = 0.0
    for row in range(HELIO_LAT_TABLE.shape[1]):
        b0 += (HELIO_LAT_TABLE[0, row, 0]
               * np.cos(HELIO_LAT_TABLE[0, row, 1]
                        + HELIO_LAT_TABLE[0, row, 2] * jme)
               )
        b1 += (HELIO_LAT_TABLE[1, row, 0]
               * np.cos(HELIO_LAT_TABLE[1, row, 1]
                        + HELIO_LAT_TABLE[1, row, 2] * jme)
               )

    b_rad = (b0 + b1 * jme)/10**8
    b = np.rad2deg(b_rad)
    return b 
Example 8
Project: pytorch_mpiigaze   Author: hysts   File: demo.py    License: MIT License 6 votes vote down vote up
def _draw_gaze_vector(self, face: Face) -> None:
        length = self.config.demo.gaze_visualization_length
        if self.config.mode == GazeEstimationMethod.MPIIGaze.name:
            for key in [FacePartsName.REYE, FacePartsName.LEYE]:
                eye = getattr(face, key.name.lower())
                self.visualizer.draw_3d_line(
                    eye.center, eye.center + length * eye.gaze_vector)
                pitch, yaw = np.rad2deg(eye.vector_to_angle(eye.gaze_vector))
                logger.info(
                    f'[{key.name.lower()}] pitch: {pitch:.2f}, yaw: {yaw:.2f}')
        elif self.config.mode == GazeEstimationMethod.MPIIFaceGaze.name:
            self.visualizer.draw_3d_line(
                face.center, face.center + length * face.gaze_vector)
            pitch, yaw = np.rad2deg(face.vector_to_angle(face.gaze_vector))
            logger.info(f'[face] pitch: {pitch:.2f}, yaw: {yaw:.2f}')
        else:
            raise ValueError 
Example 9
Project: GraphicDesignPatternByPython   Author: Relph1119   File: polar.py    License: MIT License 6 votes vote down vote up
def __call__(self, x, pos=None):
        vmin, vmax = self.axis.get_view_interval()
        d = np.rad2deg(abs(vmax - vmin))
        digits = max(-int(np.log10(d) - 1.5), 0)

        if rcParams['text.usetex'] and not rcParams['text.latex.unicode']:
            format_str = r"${value:0.{digits:d}f}^\circ$"
            return format_str.format(value=np.rad2deg(x), digits=digits)
        else:
            # we use unicode, rather than mathtext with \circ, so
            # that it will work correctly with any arbitrary font
            # (assuming it has a degree sign), whereas $5\circ$
            # will only work correctly with one of the supported
            # math fonts (Computer Modern and STIX)
            format_str = "{value:0.{digits:d}f}\N{DEGREE SIGN}"
            return format_str.format(value=np.rad2deg(x), digits=digits) 
Example 10
Project: NeuroKit   Author: neuropsychology   File: complexity_delay.py    License: MIT License 5 votes vote down vote up
def _embedding_delay_select(metric_values, algorithm="first local minimum"):

    if algorithm == "first local minimum":
        optimal = signal_findpeaks(-1 * metric_values, relative_height_min=0.1, relative_max=True)["Peaks"][0]
    elif algorithm == "first 1/e crossing":
        metric_values = metric_values - 1 / np.exp(1)
        optimal = signal_zerocrossings(metric_values)[0]
    elif algorithm == "first zero crossing":
        optimal = signal_zerocrossings(metric_values)[0]
    elif algorithm == "closest to 40% of the slope":
        slope = np.diff(metric_values) * len(metric_values)
        slope_in_deg = np.rad2deg(np.arctan(slope))
        optimal = np.where(slope_in_deg == find_closest(40, slope_in_deg))[0][0]
    return optimal 
Example 11
Project: simnibs   Author: simnibs   File: mesh_io.py    License: GNU General Public License v3.0 5 votes vote down vote up
def triangle_angles(self):
        ''' Calculates minimum angle of the mesh triangles
        
        Returns
        --------
        min_triangle_angles: ElementData
            ElementData field with the angles for each triangle, in degrees
        '''
        tr_indexes = self.elm.triangles
        node_tr = self.nodes[self.elm[tr_indexes, :3]]
        a = np.linalg.norm(node_tr[:, 1] - node_tr[:, 0], axis=1)
        b = np.linalg.norm(node_tr[:, 2] - node_tr[:, 0], axis=1)
        c = np.linalg.norm(node_tr[:, 2] - node_tr[:, 1], axis=1)
        cos_angles = np.vstack([
            (a**2 + b**2 - c**2) / (2*a*b),
            (a**2 + c**2 - b**2) / (2*a*c),
            (b**2 + c**2 - a**2) / (2*b*c),
        ]).T
        angles = np.rad2deg(np.arccos(cos_angles))

        angle_field  = ElementData(
            np.nan*np.zeros((self.elm.nr, 3)),
            'min_triangle_angles'
        )
        angle_field[tr_indexes] = angles
        return angle_field 
Example 12
Project: 2D-Motion-Retargeting   Author: ChrisWu1997   File: datasets.py    License: MIT License 5 votes vote down vote up
def get_cluster_data(self, nr_motions=50):
        """
        get a certain subset data for clustering visualization and scoring
        :param nr_anims:
        :return: pre-processed data of shape (nr_views, nr_anims, 30, 64)
        """
        motion_items = self.motion_names
        if nr_motions < len(self.motion_names):
            idxes = np.linspace(0, len(self.motion_names) - 1, nr_motions, dtype=int)
            motion_items = [self.motion_names[i] for i in idxes]
        motion_names = np.array([item[0] for item in motion_items])

        all_data = []
        for mot in motion_items:
            char_data = []
            for char in self.character_names:
                item = self.build_item(mot, char)
                view_data = []
                for view in self.view_angles:
                    data = self.preprocessing(item, view, None)
                    view_data.append(data)
                char_data.append(torch.stack(view_data, dim=0))
            all_data.append(torch.stack(char_data, dim=0))
        all_data = torch.stack(all_data, dim=0)

        ret = (all_data, motion_names, self.character_names, np.rad2deg(self.view_angles))
        return ret 
Example 13
Project: typhon   Author: atmtools   File: geodesy.py    License: MIT License 5 votes vote down vote up
def great_circle_distance(lat1, lon1, lat2, lon2, r=None):
    r"""Calculate the distance between two geographical positions

    This is the 'As-the-crow-flies' distance between two points, specified by
    their latitude and longitude. This uses the so-called haversine formula,
    defined by:

    .. math::
        d = 2 \arcsin \left( \sqrt{
            \sin \left(\frac{\varphi_2-\varphi_1}{2} \right)^2
            + \cos(\varphi_1) \cdot \cos(\varphi_1)
            \cdot \sin \left(\frac{\lambda_2-\lambda_1}{2} \right)^2} \right)

    Args:
        lat1: Latitude of position 1.
        lon1: Longitude of position 1.
        lat2: Latitude of position 2.
        lon2: Longitude of position 2.
        r (float): The radius (common for both points).

    Returns:
        If the optional argument *r* is given, the distance in m is returned.
        Otherwise the angular distance in degrees is returned.

    .. Taken from https://stackoverflow.com/a/29546836/9144990
    """
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat / 2.0) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(
        dlon / 2.0) ** 2

    c = 2 * np.arcsin(np.sqrt(a))

    if r is None:
        return np.rad2deg(c)
    else:
        return r * c 
Example 14
Project: typhon   Author: atmtools   File: aster.py    License: MIT License 5 votes vote down vote up
def geocentric2geodetic(latitude):
    """Translate geocentric to geodetic latitudes.

    Parameters:
        latitude (ndarray): latitude values in degree.

    Returns:
        (ndarray): geodetic latitudes.
    """

    return np.rad2deg(np.arctan(1.0067395 * np.tan(np.deg2rad(latitude)))) 
Example 15
Project: typhon   Author: atmtools   File: plot_phase.py    License: MIT License 5 votes vote down vote up
def degree_formatter(x, pos):
    """Create degree ticklabels for radian data."""
    return '{:.0f}\N{DEGREE SIGN}'.format(np.rad2deg(x))


# Read wind data. 
Example 16
Project: mars   Author: mars-project   File: rad2deg.py    License: Apache License 2.0 5 votes vote down vote up
def rad2deg(x, out=None, where=None, **kwargs):
    """
    Convert angles from radians to degrees.

    Parameters
    ----------
    x : array_like
        Angle in radians.
    out : Tensor, None, or tuple of Tensor and None, optional
        A location into which the result is stored. If provided, it must have
        a shape that the inputs broadcast to. If not provided or `None`,
        a freshly-allocated tensor is returned. A tuple (possible only as a
        keyword argument) must have length equal to the number of outputs.
    where : array_like, optional
        Values of True indicate to calculate the ufunc at that position, values
        of False indicate to leave the value in the output alone.
    **kwargs

    Returns
    -------
    y : Tensor
        The corresponding angle in degrees.

    See Also
    --------
    deg2rad : Convert angles from degrees to radians.

    Notes
    -----
    rad2deg(x) is ``180 * x / pi``.

    Examples
    --------
    >>> import mars.tensor as mt

    >>> mt.rad2deg(mt.pi/2).execute()
    90.0
    """
    op = TensorRad2deg(**kwargs)
    return op(x, out=out, where=where) 
Example 17
Project: callhorizons   Author: mommermi   File: test_callhorizons.py    License: MIT License 5 votes vote down vote up
def test_pyephem():
    """ test PyEphem interface, if PyEphem is available """
    
    try:
        import ephem
    except ImportError:
        return None

    asteroid_targetname = 'Ceres'
    target = callhorizons.query(asteroid_targetname)

    target.set_discreteepochs([2451544.500000])

    # create pyephem object and calculate ra and dec for a given date
    ephobj = target.export2pyephem()[0]
    ephobj.compute('2000-01-10')
    # retrieve astrometric geocentric coordinates
    ephem_ra  = np.rad2deg(ephobj.a_ra)
    ephem_dec = np.rad2deg(ephobj.a_dec)
    ephem_mag = ephobj.mag
    
    # query horizons for exact positions and check if residuals are negligible
    # i.e., less than 0.5 arcsec
    target = callhorizons.query(asteroid_targetname)
    target.set_discreteepochs([2451553.5])
    target.get_ephemerides(500)

    assert ((target['RA'][0]-ephem_ra)*3600) < 0.1
    assert ((target['DEC'][0]-ephem_dec)*3600) < 0.1
    assert (target['V'][0]-ephem_mag) < 0.1 
Example 18
def get_glasses_info(self, face_shape, face_width):
        """
        获取当前面部的眼镜信息
        :param face_shape:
        :param face_width:
        :return:
        """
        left_eye = face_shape[36:42]
        right_eye = face_shape[42:48]

        left_eye_center = left_eye.mean(axis=0).astype("int")
        right_eye_center = right_eye.mean(axis=0).astype("int")

        y = left_eye_center[1] - right_eye_center[1]
        x = left_eye_center[0] - right_eye_center[0]
        eye_angle = np.rad2deg(np.arctan2(y, x))

        deal = self.deal.resize(
            (face_width, int(face_width * self.deal.size[1] / self.deal.size[0])),
            resample=Image.LANCZOS)

        deal = deal.rotate(eye_angle, expand=True)
        deal = deal.transpose(Image.FLIP_TOP_BOTTOM)

        left_eye_x = left_eye[0, 0] - face_width // 4
        left_eye_y = left_eye[0, 1] - face_width // 6

        return {"image": deal, "pos": (left_eye_x, left_eye_y)} 
Example 19
Project: face-detection-induction-course   Author: tomoncle   File: input_video_stream_paste_mask.py    License: MIT License 5 votes vote down vote up
def get_glasses_info(self, face_shape, face_width):
        """
        获取当前面部的眼镜信息
        :param face_shape:
        :param face_width:
        :return:
        """
        left_eye = face_shape[36:42]
        right_eye = face_shape[42:48]

        left_eye_center = left_eye.mean(axis=0).astype("int")
        right_eye_center = right_eye.mean(axis=0).astype("int")

        y = left_eye_center[1] - right_eye_center[1]
        x = left_eye_center[0] - right_eye_center[0]
        eye_angle = np.rad2deg(np.arctan2(y, x))

        deal = self.deal.resize(
            (face_width, int(face_width * self.deal.size[1] / self.deal.size[0])),
            resample=Image.LANCZOS)

        deal = deal.rotate(eye_angle, expand=True)
        deal = deal.transpose(Image.FLIP_TOP_BOTTOM)

        left_eye_x = left_eye[0, 0] - face_width // 4
        left_eye_y = left_eye[0, 1] - face_width // 6

        return {"image": deal, "pos": (left_eye_x, left_eye_y)} 
Example 20
Project: satellite_tracker   Author: PaulKlinger   File: orbit_np.py    License: GNU General Public License v3.0 5 votes vote down vote up
def get_lonlatalt(self, utc_time):
        """Calculate sublon, sublat and altitude of satellite.
        http://celestrak.com/columns/v02n03/
        """
        (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z), errors = self.get_position(
            utc_time, normalize=True)

        lon = ((np.arctan2(pos_y * XKMPER, pos_x * XKMPER) - gmst(dt2np(utc_time)))
               % (2 * np.pi))

        lon = np.where(lon > np.pi, lon - np.pi * 2, lon)
        lon = np.where(lon <= -np.pi, lon + np.pi * 2, lon)

        r = np.sqrt(pos_x ** 2 + pos_y ** 2)
        lat = np.arctan2(pos_z, r)
        e2 = F * (2 - F)
        while True:
            lat2 = lat
            c = 1 / (np.sqrt(1 - e2 * (np.sin(lat2) ** 2)))
            lat = np.arctan2(pos_z + c * e2 * np.sin(lat2), r)
            if np.all(abs(lat - lat2) < 1e-5):
                break
        alt = r / np.cos(lat) - c
        alt *= A
        return np.rad2deg(lon), np.rad2deg(lat), alt, errors


# t1 = Tle("a", line1="1     5U 58002B   18231.96948527 -.00000019 +00000-0 -48070-4 0  9990", line2="2     5 034.2557 124.2521 1846373 229.0197 113.5561 10.84817833132766")
# t2 = Tle("b", line1="1    11U 59001A   18234.91021469 +.00000021 +00000-0 +14392-4 0  9993", line2="2    11 032.8682 030.1273 1466645 100.3180 276.5496 11.85533996196087")
# ts = np.array([t1,t2])

# o=OrbitElements(ts) 
Example 21
Project: D-VAE   Author: muhanzhang   File: test_var.py    License: MIT License 5 votes vote down vote up
def test_numpy_method():
    # This type of code is used frequently by PyMC3 users
    x = tt.dmatrix('x')
    data = np.random.rand(5, 5)
    x.tag.test_value = data
    for fct in [np.arccos, np.arccosh, np.arcsin, np.arcsinh,
                np.arctan, np.arctanh, np.ceil, np.cos, np.cosh, np.deg2rad,
                np.exp, np.exp2, np.expm1, np.floor, np.log,
                np.log10, np.log1p, np.log2, np.rad2deg,
                np.sin, np.sinh, np.sqrt, np.tan, np.tanh, np.trunc]:
        y = fct(x)
        f = theano.function([x], y)
        utt.assert_allclose(np.nan_to_num(f(data)),
                            np.nan_to_num(fct(data))) 
Example 22
Project: D-VAE   Author: muhanzhang   File: basic.py    License: MIT License 5 votes vote down vote up
def impl(self, x):
        # If x is an int8 or uint8, numpy.rad2deg will compute the result in
        # half-precision (float16), where we want float32.
        x_dtype = str(getattr(x, 'dtype', ''))
        if x_dtype in ('int8', 'uint8'):
            return numpy.rad2deg(x, sig='f')
        return numpy.rad2deg(x) 
Example 23
Project: TheCannon   Author: annayqho   File: simpletable.py    License: MIT License 5 votes vote down vote up
def sphdist(ra1, dec1, ra2, dec2):
        """measures the spherical distance between 2 points

        Parameters
        ----------
        ra1: float or sequence
            first right ascensions in degrees

        dec1: float or sequence
            first declination in degrees
        ra2: float or sequence
            second right ascensions in degrees
        dec2: float or sequence
            first declination in degrees

        Returns
        -------
        Outputs: float or sequence
            returns a distance in degrees
        """
        dec1_r = deg2rad(dec1)
        dec2_r = deg2rad(dec2)
        return 2. * rad2deg(arcsin(sqrt((sin((dec1_r - dec2_r) / 2)) ** 2 +
                                        cos(dec1_r) * cos(dec2_r) * (
                                            sin((deg2rad(ra1 - ra2)) / 2)) **
                                        2))) 
Example 24
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def heliocentric_longitude(jme):
        """Compute the Earth Heliocentric Longitude (L) in degrees given the Julian Ephemeris Millennium"""
        #L5, ..., L0
        Li = [sum(a*np.cos(b + c*jme) for a,b,c in abcs) for abcs in reversed(_sp._EHL_)]
        L = np.polyval(Li, jme) / 1e8
        L = np.rad2deg(L) % 360
        return L 
Example 25
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def heliocentric_latitude(jme):
        """Compute the Earth Heliocentric Latitude (B) in degrees given the Julian Ephemeris Millennium"""
        Bi = [sum(a*np.cos(b + c*jme) for a,b,c in abcs) for abcs in reversed(_sp._EHB_)]
        B = np.polyval(Bi, jme) / 1e8
        B = np.rad2deg(B) % 360
        return B 
Example 26
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def nutation_obliquity(jce):
        """compute the nutation in longitude (delta_psi) and the true obliquity (epsilon) given the Julian Ephemeris Century"""
        
        #mean elongation of the moon from the sun, in radians:
        #x0 = 297.85036 + 445267.111480*jce - 0.0019142*(jce**2) + (jce**3)/189474
        x0 = np.deg2rad(np.polyval([1./189474, -0.0019142, 445267.111480, 297.85036],jce))
        #mean anomaly of the sun (Earth), in radians:
        x1 = np.deg2rad(np.polyval([-1/3e5, -0.0001603, 35999.050340, 357.52772], jce))
        #mean anomaly of the moon, in radians:
        x2 = np.deg2rad(np.polyval([1./56250, 0.0086972, 477198.867398, 134.96298], jce))
        #moon's argument of latitude, in radians:
        x3 = np.deg2rad(np.polyval([1./327270, -0.0036825, 483202.017538, 93.27191], jce))
        #Longitude of the ascending node of the moon's mean orbit on the ecliptic
        # measured from the mean equinox of the date, in radians
        x4 = np.deg2rad(np.polyval([1./45e4, 0.0020708, -1934.136261, 125.04452], jce))

        x = (x0, x1, x2, x3, x4)

        dp = 0.0
        for y, ab in zip(_sp._NLOY_, _sp._NLOab_):
            a,b = ab
            dp += (a + b*jce)*np.sin(np.dot(x, y))
        dp = np.rad2deg(dp)/36e6

        de = 0.0
        for y, cd in zip(_sp._NLOY_, _sp._NLOcd_):
            c,d = cd
            de += (c + d*jce)*np.cos(np.dot(x, y))
        de = np.rad2deg(de)/36e6

        e = _sp.ecliptic_obliquity(_sp.julian_millennium(jce), de)

        return dp, e 
Example 27
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def sun_ra_decl(llambda, epsilon, beta):
        """Calculate the sun's geocentric right ascension (alpha, in degrees) and declination (delta, in degrees)"""
        l, e, b = map(np.deg2rad, (llambda, epsilon, beta))
        alpha = np.arctan2(np.sin(l)*np.cos(e) - np.tan(b)*np.sin(e), np.cos(l)) #x1 / x2
        alpha = np.rad2deg(alpha) % 360
        delta = np.arcsin(np.sin(b)*np.cos(e) + np.cos(b)*np.sin(e)*np.sin(l))
        delta = np.rad2deg(delta)
        return alpha, delta 
Example 28
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def sun_topo_ra_decl_hour(latitude, longitude, elevation, jd, delta_t = 0):
        """Calculate the sun's topocentric right ascension (alpha'), declination (delta'), and hour angle (H')"""
        
        jde = _sp.julian_ephemeris_day(jd, delta_t)
        jce = _sp.julian_century(jde)
        jme = _sp.julian_millennium(jce)

        helio_pos = _sp.heliocentric_position(jme)
        R = helio_pos[-1]
        phi, sigma, E = latitude, longitude, elevation
        #equatorial horizontal parallax of the sun, in radians
        xi = np.deg2rad(8.794/(3600*R)) #
        #rho = distance from center of earth in units of the equatorial radius
        #phi-prime = geocentric latitude
        #NB: These equations look like their based on WGS-84, but are rounded slightly
        # The WGS-84 reference ellipsoid has major axis a = 6378137 m, and flattening factor 1/f = 298.257223563
        # minor axis b = a*(1-f) = 6356752.3142 = 0.996647189335*a
        u = np.arctan(0.99664719*np.tan(phi)) #
        x = np.cos(u) + E*np.cos(phi)/6378140 #rho sin(phi-prime)
        y = 0.99664719*np.sin(u) + E*np.sin(phi)/6378140 #rho cos(phi-prime)

        delta_psi, epsilon = _sp.nutation_obliquity(jce) #

        llambda, beta = _sp.sun_longitude(helio_pos, delta_psi) #
        
        alpha, delta = _sp.sun_ra_decl(llambda, epsilon, beta) #

        v = _sp.greenwich_sidereal_time(jd, delta_psi, epsilon) #

        H = v + longitude - alpha #
        Hr, dr = map(np.deg2rad,(H,delta))

        dar = np.arctan2(-x*np.sin(xi)*np.sin(Hr), np.cos(dr)-x*np.sin(xi)*np.cos(Hr))
        delta_alpha = np.rad2deg(dar) #
        
        alpha_prime = alpha + delta_alpha #
        delta_prime = np.rad2deg(np.arctan2((np.sin(dr) - y*np.sin(xi))*np.cos(dar), np.cos(dr) - y*np.sin(xi)*np.cos(Hr))) #
        H_prime = H - delta_alpha #

        return alpha_prime, delta_prime, H_prime 
Example 29
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def norm_lat_lon(lat,lon):
        if lat < -90 or lat > 90:
            #convert to cartesian and back
            x = cos(deg2rad(lon))*cos(deg2rad(lat))
            y = sin(deg2rad(lon))*cos(deg2rad(lat))
            z = sin(deg2rad(lat))
            r = sqrt(x**2 + y**2 + z**2)
            lon = rad2deg(arctan2(y,x)) % 360
            lat = rad2deg(arcsin(z/r))
        elif lon < 0 or lon > 360:
            lon = lon % 360
        return lat,lon 
Example 30
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def arcdist(p0,p1,radians=False):
    """Angular distance between azimuth,zenith pairs
    
    Parameters
    ----------
    p0 : array_like, shape (..., 2)
    p1 : array_like, shape (..., 2)
        p[...,0] = azimuth angles, p[...,1] = zenith angles
    radians : boolean (default False)
        If False, angles are in degrees, otherwise in radians

    Returns
    -------
    ad :  array_like, shape is broadcast(p0,p1).shape
        Arcdistances between corresponding pairs in p0,p1
        In degrees by default, in radians if radians=True
    """
    #formula comes from translating points into cartesian coordinates
    #taking the dot product to get the cosine between the two vectors
    #then arccos to return to angle, and simplify everything assuming real inputs
    p0,p1 = np.array(p0), np.array(p1)
    if not radians:
        p0,p1 = np.deg2rad(p0), np.deg2rad(p1)
    a0,z0 = p0[...,0], p0[...,1]
    a1,z1 = p1[...,0], p1[...,1]
    d = np.arccos(np.cos(z0)*np.cos(z1)+np.cos(a0-a1)*np.sin(z0)*np.sin(z1))
    if radians:
        return d
    else:
        return np.rad2deg(d)