Python EarthLocation.to_geodetic方法代码示例

# 需要导入模块: from astropy.coordinates import EarthLocation [as 别名]
# 或者: from astropy.coordinates.EarthLocation import to_geodetic [as 别名]
def calculateRelativeVelocity(obs_coords,psr_coords,MJD):
    """
    obs_coords in (x,y,z) from PSRFITS file
    psr_coords in (RA,DEC) from PSRFITS file (both strings)
    Returns relative velocity in km/s to first-order
    """

    t = Time(MJD, format='mjd')
    time = t.datetime.strftime("%Y/%m/%d %H:%M:%S")
    


    el = EarthLocation(x=obs_coords[0]*units.meter,y=obs_coords[1]*units.meter,z=obs_coords[2]*units.meter)#,frame='itrs')
    lon,lat,height = el.to_geodetic()

    #spin_vel = calculateSpinVelocity(lon,lat,t)

    obs = ephem.Observer()
    obs.lon = ephem.degrees(str(lon.value))
    obs.lat = ephem.degrees(str(lat.value))
    obs.date = time
    obs.pressure = 0 #no refraction

    # From pyephem docs:
    # Both hlon and hlat have a special meaning for the Sun and Moon. For a Sun body, they give the Earth's heliocentric longitude and latitude
    # Note: Sun at vernal equinox has hlong = 0, hlat = 0, Earth will return hlong = 180 degrees! This has been verified

    sun = ephem.Sun()
    sun.compute(time,epoch=2000)

    vel_ra = ephem.degrees(sun.hlong-3*np.pi/2)
    vel_dec = ephem.degrees(-1*sun.hlat)
    pulsar_ra = ephem.hours(psr_coords[0])
    pulsar_dec = ephem.degrees(psr_coords[1])

    sep = ephem.separation((vel_ra,vel_dec),(pulsar_ra,pulsar_dec))

    return 30*np.cos(sep) 

# 需要导入模块: from astropy.coordinates import EarthLocation [as 别名]
# 或者: from astropy.coordinates.EarthLocation import to_geodetic [as 别名]

#.........这里部分代码省略.........
    # by the International Astronomical Union.
    # See http://en.wikipedia.org/wiki/Geographic_coordinate_conversion
    antx = antposition.x.value # in meters
    anty = antposition.y.value # in meters
    antz = antposition.z.value # in meters

    # Find the unit normal vector to the ellipsoid at the antenna position.
    # Note that the WGS84 oblate spheroid is defined as:
    #
    # x^2 + y^2   z^2
    # --------- + --- = 1, with:
    #    a^2      b^2
    #
    # a = 6 378 137.0 m
    # b = 6 356 752.314 245 m
    # See http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
    a = 6378137.0
    b = 6356752.314245

    normalx = 2*antx/a**2.0
    normaly = 2*anty/a**2.0
    normalz = 2*antz/b**2.0

    unitnormalx = normalx/np.sqrt(normalx**2.0+normaly**2.0+normalz**2.0)
    unitnormaly = normaly/np.sqrt(normalx**2.0+normaly**2.0+normalz**2.0)
    unitnormalz = normalz/np.sqrt(normalx**2.0+normaly**2.0+normalz**2.0)

    # Parameterize Az, El line-of-sight towards the object in the sky.
    # Note that instead of Elevation, the inclination (polar angle) is used.
    # Also, note the minus sign because of the azimuth being definined with
    # the x-direction as north and y-direction as west.
    inclination = np.pi/2 - el

    param_x = np.sin(inclination)*np.cos(az)
    param_y = -np.sin(inclination)*np.sin(az)
    param_z = np.cos(inclination)

    # Now, rotate and translate the horizontal frame in such a way that it
    # is tangent to the Earth at the antenna position and that the
    # x-axis points towards the north in the ECEF frame.
    # This comes down to rotating and translating the parameterization of
    # the line-of-sight using theta and phi from the normal vector.
    # The angle between the north in horizontal frame and the unit normal
    # vector is theta. Phi is the angle between the x-direction in the
    # horizontal frame and the projection of the normal vector on the
    # xy-plane.
    theta = np.arccos(unitnormalz)
    phi = np.arctan(unitnormaly/unitnormalx)

    rotation_matrix_z = np.matrix( [[np.cos(phi), -np.sin(phi), 0],
                    [np.sin(phi), np.cos(phi), 0],
                    [0, 0, 1]] )
    rotation_matrix_y = np.matrix( [[np.cos(theta), 0, -np.sin(theta)],
                    [0, 1, 0],
                    [np.sin(theta), 0, np.cos(theta)]] )

    param_vector = np.array([param_x, param_y, param_z])[:,None]
    intermediate_step = np.dot(rotation_matrix_y, param_vector)
    rotated_param = np.dot(rotation_matrix_z, intermediate_step)
    rotated_param = np.squeeze(np.asarray(rotated_param))

    rotated_param_x = rotated_param[0]
    rotated_param_y = rotated_param[1]
    rotated_param_z = rotated_param[2]

    # This rotated parameterization can now be translated using the
    # antenna position. This gives a parameterization of the line-
    # of-sight in the ECEF system. The intersection with a larger
    # ellipsoid (given by the screen height) gives the pierce point.
    # The parameterization is given by:
    #     | antx |       | rotated_param_x |
    # l = | anty | + t * | rotated_param_y |
    #     | anyz |       | rotated_param_x |
    # whereas the shell representing the ionosphere is given by:
    #
    # x^2 + y^2   z^2
    # --------- + --- = 1, with:
    #   c^2      d^2
    # with c = a + height in meters
    # and d = b + height in meters

    c = a + screenheight
    d = b + screenheight

    t = -(1./(d**2.*(rotated_param_x**2. + rotated_param_y**2.) +
        c**2.*rotated_param_z**2.)) * (d**2.*(rotated_param_x*antx +
        rotated_param_y*anty) + c**2.*rotated_param_z*antz - 1/2.*np.sqrt(
        4*(d**2.*(rotated_param_x*antx + rotated_param_y*anty) +
        c**2.*rotated_param_z*antz)**2. - 4*(d**2.*(rotated_param_x**2. +
        rotated_param_y**2.) + c**2.*rotated_param_z**2.)*(d**2.*(-c**2. +
        antx**2. + anty**2.) + c**2.*antz**2.)))

    piercex = antx + t*rotated_param_x
    piercey = anty + t*rotated_param_y
    piercez = antz + t*rotated_param_z

    piercepoint_xyz = EarthLocation(x=piercex*u.m, y=piercey*u.m, z=piercez*u.m)
    piercepoint_lonlatheight = piercepoint_xyz.to_geodetic()

    return piercepoint_xyz

# 需要导入模块: from astropy.coordinates import EarthLocation [as 别名]
# 或者: from astropy.coordinates.EarthLocation import to_geodetic [as 别名]
    #STT_SMJD=                36103 / [s] Start time (sec past UTC 00h) (J) 

    ANT_X = 882589.65 #* units.meter
    ANT_Y = -4924872.32# * units.meter
    ANT_Z = 3943729.348# * units.meter
    RA = '19:09:47.448'
    DEC = '-37:44:13.920'
    STT_IMJD = 55275
    STT_SMJD = 36103

    MJD = STT_IMJD+STT_SMJD/86400.0
    MJD = 57184.5+11 #this was the test date I used approximately on the solstice

    t = Time(MJD, format='mjd')
    time = t.datetime.strftime("%Y/%m/%d %H:%M:%S")


    utc_hour = t.datetime.hour+t.datetime.minute/60.0+t.datetime.second/3600.0

    el = EarthLocation(x=ANT_X*units.meter,y=ANT_Y*units.meter,z=ANT_Z*units.meter)#,frame='itrs')
    lon,lat,height = el.to_geodetic()
    
    local_hour = utc_hour + lon.hourangle


    raise SystemExit

    #print calculateRelativeVelocity((ANT_X,ANT_Y,ANT_Z),(RA,DEC),MJD)