Python coordinates.Angle类代码示例

    def __init__(
        self, distance, inclination, azimuth, position_angle, pixels_x, pixels_y, center_x, center_y, field_x, field_y
    ):

        """
        This function ...
        :param distance:
        :param inclination:
        :param azimuth:
        :param position_angle:
        :param pixels_x:
        :param pixels_y:
        :param center_x:
        :param center_y:
        :param field_x:
        :param field_y:
        """

        self.distance = distance
        self.inclination = Angle(inclination, "deg")
        self.azimuth = Angle(azimuth, "deg")
        self.position_angle = Angle(position_angle, "deg")
        self.pixels_x = pixels_x
        self.pixels_y = pixels_y
        self.center_x = center_x
        self.center_y = center_y
        self.field_x_physical = field_x
        self.field_y_physical = field_y

def test_3rd_body_Curtis(test_params):
    # based on example 12.11 from Howard Curtis
    body = test_params['body']
    with solar_system_ephemeris.set('builtin'):
        j_date = 2454283.0 * u.day
        tof = (test_params['tof']).to(u.s).value
        body_r = build_ephem_interpolant(body, test_params['period'], (j_date, j_date + test_params['tof']), rtol=1e-2)

        epoch = Time(j_date, format='jd', scale='tdb')
        initial = Orbit.from_classical(Earth, *test_params['orbit'], epoch=epoch)
        r, v = cowell(initial, np.linspace(0, tof, 400), rtol=1e-10, ad=third_body,
                      k_third=body.k.to(u.km**3 / u.s**2).value, third_body=body_r)

        incs, raans, argps = [], [], []
        for ri, vi in zip(r, v):
            angles = Angle(rv2coe(Earth.k.to(u.km**3 / u.s**2).value, ri, vi)[2:5] * u.rad)  # inc, raan, argp
            angles = angles.wrap_at(180 * u.deg)
            incs.append(angles[0].value)
            raans.append(angles[1].value)
            argps.append(angles[2].value)

        # averaging over 5 last values in the way Curtis does
        inc_f, raan_f, argp_f = np.mean(incs[-5:]), np.mean(raans[-5:]), np.mean(argps[-5:])

        assert_quantity_allclose([(raan_f * u.rad).to(u.deg) - test_params['orbit'][3],
                                  (inc_f * u.rad).to(u.deg) - test_params['orbit'][2],
                                  (argp_f * u.rad).to(u.deg) - test_params['orbit'][4]],
                                 [test_params['raan'], test_params['inc'], test_params['argp']],
                                 rtol=1e-1)

def deg_to_sex(ra, dec):
    """
    Convert an Ra and Dec position in decimal degrees to hexadecimal
    Parameters
    ----------
    ra : float
        Ra in decimal degrees
    dec : float
        Dec in decimal degrees

    Returns
    -------
    ra : str
        Ra in hexadecimal HH:MM:SS

    dec : str
        Dec in hexadecimal DD:MM:SS
    """
    from astropy import units as u
    from astropy.coordinates import Angle
    rad = Angle(ra * u.deg)
    decd = Angle(dec * u.deg)

    ra = rad.to_string(unit=u.hour, sep=':')
    dec = decd.to_string(unit=u.deg, sep=':')

    return ra, dec

def reformat_catalog(coordinates=True, chi=True, sharp=True, bright=True):
	"""Format RA, Dec as floats in the final catalog"""
	
	t = Table.read('../data/extract.txt', format='ascii.commented_header')
	
	if coordinates:
		ra = Angle((t['rah'], t['ram'], t['ras']), unit=u.hour)
		dec = Angle((t['decd'], t['decm'], t['decs']), unit=u.deg)
		
		t.remove_columns(['rah', 'ram', 'ras', 'decd', 'decm', 'decs'])
		add_npcolumn(t, 'ra', ra.to(u.deg), index=0)
		add_npcolumn(t, 'dec', dec, index=1)
	
	if bright:
		ind = ((t['g']<18) | (t['i']<18)) & (t['ichi']<90)
		t_bright = t[ind]
		t = t[~ind]
		
	if chi:
		#ind = (t['gchi']<1.4) & (t['ichi']<1.4)
		ind = (t['gchi']<1.5) & (t['ichi']<1.5)
		t = t[ind]
		
	if sharp:
		#ind = (np.abs(t['gsharp'])<0.2) & (np.abs(t['isharp'])<0.2)
		ind = (np.abs(t['gsharp'])<0.3) & (np.abs(t['isharp'])<0.3)
		t = t[ind]
	
	if bright:
		t = astropy.table.vstack((t, t_bright))
	
	t.pprint()
	t.write('../data/ext_catalog.txt', format='ascii.commented_header')

def axial_ratio_to_inclination_with_intrinsic(ratio, intrinsic_ratio):

    """
    This function ...
    :param ratio:
    :param intrinsic_ratio:
    :return:
    """

    # Check ratios
    if ratio >= 1.: raise ValueError("Axial ratio must be ratio of minor axis to major axis length")
    if intrinsic_ratio >= 1: raise ValueError("Intrinsic axial ratio must be ratio of scale height (smaller) to scale length (larger)")

    # Calculate logq and logq0
    logq = np.log10(ratio)
    logq0 = np.log10(intrinsic_ratio)

    # Calculate numerator and denominator of the formula
    numerator = 1. - 10 ** (2. * logq)
    denominator = 1. - 10 ** (2. * logq0)

    # Calculate the inclination angle
    sin2i = numerator / denominator
    inclination_radians = np.arcsin(np.sqrt(sin2i))
    inclination = Angle(inclination_radians, unit="rad")
    return inclination.to("deg")

def axial_ratio_to_inclination_mosenkov(ratio, hubble_stage):

    """
    This function ...
    :param ratio:
    :param hubble_stage:
    :return:
    """

    # Check that axial ratio is smaller than one!
    if ratio >= 1.: raise ValueError("Axial ratio must be ratio of minor axis to major axis length")

    # From Mosenkov et al., 2017 (Appendix)

    # Calculate logq and logq0
    logq = np.log10(ratio)
    logq0 = get_logq0_mosenkov(hubble_stage)

    # Calculate numerator and denominator of the formula
    numerator = 1. - 10**(2. * logq)
    denominator = 1. - 10**(2. * logq0)

    # Calculate the inclination angle
    sin2i = numerator / denominator
    inclination_radians = np.arcsin(np.sqrt(sin2i))
    inclination = Angle(inclination_radians, unit="rad")
    return inclination.to("deg")

def plot_vhs_tiles(table=None,ra=None, dec=None,
    wrap_ra24hr=False, PA=90.0, overplot=True, savefig=True):
    """

    table is astropy table
    ra, dec in degrees

    """

    plot_polygon=True

    if wrap_ra24hr:
        angle=Angle(xdata * u.deg)
        angle.wrap_at('180d', inplace=True)
        xdata=angle.degree

    # assumes filename is in table.meta
    print('plotting: ', table.meta['filename'])
    filename = table.meta['filename']

    if not plot_polygon:
      plt.plot(xdata, ydata, 's', ms=7.0,
       color='yellow', markeredgecolor='yellow',
       alpha=0.1,
       label='OB Progress (submitted)\n' + filename)

    if plot_polygon:
      plt.plot(xdata, ydata, '.', ms=7.0,
       color='yellow', markeredgecolor='yellow',
       alpha=0.1,
       label='OB Progress (submitted):' + filename)

      # filled polygon
      ra=xdata
      dec=ydata
      for i in range(len(ra)):
        if i==0 or i == len(ra): print('i: ', i)
        ra_poly, dec_poly=plt_tile.mk_polytile(ra_cen=ra[i], dec_cen=dec[i],
         coverage='twice', PA=PA)

        xypolygon=np.column_stack((ra_poly, dec_poly))
        if i==0 or i == len(ra):
          print('xypolygon.shape: ',   xypolygon.shape)
          print('xypolygon: ',xypolygon)

        polygon = Polygon(xypolygon, True, color='green', alpha=0.1)
        plt.gca().add_patch(polygon)

        #print(ob_table.meta)

      plt.suptitle(dqcfile_vhs)
      plt.legend(fontsize='small')

      figname='vhs_des_check_progress_vhs_obprogress_'+ datestamp + '.png'
      print('Saving: '+figname)
      plt.savefig(plotdir+figname)
      plt.suptitle('')

def deg_to_sex(ra, dec):
    from astropy import units as u
    from astropy.coordinates import Angle
    rad = Angle(ra * u.deg)
    decd = Angle(dec * u.deg)

    ra = rad.to_string(unit=u.hour, sep=':')
    dec = decd.to_string(unit=u.deg, sep=':')

    return ra, dec

def test_3rd_body_Curtis(test_params):
    # based on example 12.11 from Howard Curtis
    body = test_params["body"]
    with solar_system_ephemeris.set("builtin"):
        j_date = 2454283.0 * u.day
        tof = (test_params["tof"]).to(u.s).value
        body_r = build_ephem_interpolant(
            body,
            test_params["period"],
            (j_date, j_date + test_params["tof"]),
            rtol=1e-2,
        )

        epoch = Time(j_date, format="jd", scale="tdb")
        initial = Orbit.from_classical(Earth, *test_params["orbit"], epoch=epoch)
        rr, vv = cowell(
            Earth.k,
            initial.r,
            initial.v,
            np.linspace(0, tof, 400) * u.s,
            rtol=1e-10,
            ad=third_body,
            k_third=body.k.to(u.km ** 3 / u.s ** 2).value,
            third_body=body_r,
        )

        incs, raans, argps = [], [], []
        for ri, vi in zip(rr.to(u.km).value, vv.to(u.km / u.s).value):
            angles = Angle(
                rv2coe(Earth.k.to(u.km ** 3 / u.s ** 2).value, ri, vi)[2:5] * u.rad
            )  # inc, raan, argp
            angles = angles.wrap_at(180 * u.deg)
            incs.append(angles[0].value)
            raans.append(angles[1].value)
            argps.append(angles[2].value)

        # averaging over 5 last values in the way Curtis does
        inc_f, raan_f, argp_f = (
            np.mean(incs[-5:]),
            np.mean(raans[-5:]),
            np.mean(argps[-5:]),
        )

        assert_quantity_allclose(
            [
                (raan_f * u.rad).to(u.deg) - test_params["orbit"][3],
                (inc_f * u.rad).to(u.deg) - test_params["orbit"][2],
                (argp_f * u.rad).to(u.deg) - test_params["orbit"][4],
            ],
            [test_params["raan"], test_params["inc"], test_params["argp"]],
            rtol=1e-1,
        )

    def formatter(self, values, spacing):

        if not isinstance(values, u.Quantity) and values is not None:
            raise TypeError("values should be a Quantities array")

        if len(values) > 0:
            if self.format is None:
                spacing = spacing.to(u.arcsec).value
                if spacing > 3600:
                    fields = 1
                    precision = 0
                elif spacing > 60:
                    fields = 2
                    precision = 0
                elif spacing > 1:
                    fields = 3
                    precision = 0
                else:
                    fields = 3
                    precision = -int(np.floor(np.log10(spacing)))
                decimal = False
                unit = u.degree
            else:
                fields = self._fields
                precision = self._precision
                decimal = self._decimal
                unit = self._unit

            if decimal:
                sep = None
            elif self._sep is not None:
                sep = self._sep
            else:
                if unit == u.degree:
                    if rcParams['text.usetex']:
                        deg = r'$^\circ$'
                    else:
                        deg = six.u('\xb0')
                    sep = (deg, "'", '"')
                else:
                    sep = ('h', 'm', 's')

            angles = Angle(values)
            string = angles.to_string(unit=unit,
                                      precision=precision,
                                      decimal=decimal,
                                      fields=fields,
                                      sep=sep).tolist()
            return string
        else:
            return []

 def __init__(self, energy, offset, data=None, data_units=""):
     self.energy = EnergyBounds(energy)
     self.offset = Angle(offset)
     if data is None:
         self.data = Quantity(np.zeros((len(energy) - 1, len(offset) - 1)), data_units)
     else:
         self.data = Quantity(data, data_units)

 def __init__(self, skyevent_dict):
     d = skyevent_dict
     self.position = SkyCoord(d['ra'], d['dec'], unit='deg')
     self.position_error = Angle(d['error'], unit='deg')
     self.timestamp = None
     if d.get('time'):
         self.timestamp = iso8601.parse_date(d['time'])

def sky_time(coord, time, rise_set=False, limalt=0*u.deg, site=EarthLocation(0.0, 0.0, 0.0), fuse=TimeDelta(0, format='sec', scale='tai')):
    """
    """
    if type(limalt) != u.quantity.Quantity:
        limalt = limalt*u.deg
    if time.isscalar == True:
        time = Time([time.iso], format='iso', scale='utc')
    coord = coord_pack(coord)
    timeut = time - fuse
    if len(time.shape) == 1:
        timeut = Time([[i] for i in timeut.jd], format='jd', scale='utc')
    timeut.delta_ut1_utc = 0
    timeut.location = site
    ra, ts = mesh_coord(coord, timeut[:,0])
    dif_h_sid = Angle(ra-ts)
    dif_h_sid.wrap_at('180d', inplace=True)
    dif_h_sol = dif_h_sid * (23.0 + 56.0/60.0 + 4.0916/3600.0) / 24.0
    dif = TimeDelta(dif_h_sol.hour*u.h, scale='tai')
    culminacao = timeut + dif
    culminacao.delta_ut1_utc = 0
    culminacao.location = site
    if (site.latitude > 0*u.deg):
        alwaysup = np.where(coord.dec >= 90*u.deg - site.latitude + limalt)
        neverup = np.where(coord.dec <= -(90*u.deg - site.latitude - limalt))
    else:
        alwaysup = np.where(coord.dec <= -(90*u.deg + site.latitude + limalt))
        neverup = np.where(coord.dec >= 90*u.deg + site.latitude - limalt)
    if rise_set == True:
        hangle_lim = np.arccos((np.cos(90.0*u.deg-limalt) - np.sin(coord.dec)*np.sin(site.latitude)) / (np.cos(coord.dec)*np.cos(site.latitude)))
        tsg_lim = Angle(ra + hangle_lim)
        dtsg_lim = tsg_lim - culminacao.sidereal_time('mean')
        dtsg_lim.wrap_at(360 * u.deg, inplace=True)
        dtsg_lim_sol = dtsg_lim * (23.0 + 56.0/60.0 + 4.0916/3600.0) / 24.0
        a = np.where(np.isnan(dtsg_lim_sol))
        dtsg_lim_sol[a] = Angle([48.0]*len(a[0])*u.hour)
        dtsg_np = TimeDelta((dtsg_lim_sol.hour*u.h))
        sunrise = culminacao - dtsg_np
        sunset = culminacao + dtsg_np
        culminacao = culminacao + fuse
        sunrise = sunrise + fuse
        sunset = sunset + fuse
        return culminacao, sunrise, sunset, alwaysup, neverup
    culminacao = culminacao + fuse
    return culminacao, alwaysup, neverup

def time_scatterplot(time, peakr, flux, proto, num, t_init, posx, posy):

	#Convert the decimal degrees to HMS strings.
	ra_deg = Angle(str(posx)+'d')
	dec_deg = Angle(str(posy)+'d')
	ra_hms = ra_deg.to_string(unit=units.hour, sep=':', precision=2)
	dec_dms = dec_deg.to_string(unit=units.degree, sep=':', precision=2)

	#Set the x-range.
	x_hi = max(time)*1.1
	x_lo = -max(time)*0.1
	plt.xlim(x_lo, x_hi)

	#Set the y-range.
	y_hi = 1.25
	y_lo = 0.75
	if np.amax(peakr) > 1.25:
		y_hi = np.amax(peakr)*1.5
	##fi
	if np.amax(peakr) < 0.75:
		y_lo = min(peakr)*0.7
	##fi
	plt.ylim(y_lo, y_hi)

	#Define the labels.
	plt.xlabel('Time (days + '+t_init+')')

	#Write important details about the source on the plot.
	plt.annotate('Source Number: '+str(num), xy = (0.05,0.2),
				 xycoords = 'axes fraction', size = 'small')
	plt.annotate('Peak Flux: '+'{:.4}'.format(flux)+' (Jy/Beam)',
				 xy = (0.05,0.15), xycoords = 'axes fraction', size = 'small')
	plt.annotate('Class: '+proto, xy = (0.05,0.1), xycoords = 'axes fraction',
				 size = 'small')
	plt.annotate('RA: '+ra_hms+',  Dec: '+dec_dms, xy = (0.05,0.05),
				 xycoords = 'axes fraction', size = 'small')

	#Plot the data.
	plt.scatter(time, peakr, alpha = 0.6, s = 60, marker = 'o', facecolor = 'r')

	#Plot reference lines.
	plt.plot((x_lo,x_hi),(1.0,1.0), '-k')
	plt.plot((x_lo,x_hi),(0.9,0.9), '--k')
	plt.plot((x_lo,x_hi),(1.1,1.1), '--k')

def averaging_params(max_decorr=0.1, frequency=(250. * 1e6 * units.Hz),
                     baseline=(870 * units.m), corr_FoV=Angle(90., units.degree),
                     latitude=hera_latitude, verbose=False, wrap_len=default_wrap_len):

    wavelength = const.c / frequency.to(1 / units.s)
    earth_rot_speed = (Angle(360, units.deg) / units.sday).to(units.arcminute / units.s)

    params = {}

    max_resolution = Angle(np.arcsin(wavelength / baseline), units.radian)
    params['max_resolution'] = max_resolution
    if verbose:
        print_desc_val('Max resolution:', max_resolution.to(units.arcminute), wrap_len)

    integration_time = (max_resolution * max_decorr).to(units.arcminute) / earth_rot_speed
    params['integration_time'] = integration_time
    if verbose:
        print_desc_val('Max integration time required to keep the decorrelation '
                       'due to time integrating under max_decorr on the longest '
                       'baselines:', integration_time, wrap_len)

    channel_width = ((const.c * max_decorr) / (baseline * np.sin(corr_FoV.to(units.rad)))).to(units.kHz)
    params['channel_width'] = channel_width
    if verbose:
        print_desc_val('Max channel width to keep the decorrelation due to channel '
                       'width under max_decorr for a {fov} degree correlator '
                       'FoV on the longest baselines:'.format(fov=corr_FoV.degree),
                       channel_width, wrap_len)

    # After fringe stopping, the rotation is from the sky rotating in the beam
    # This is slower, so we can sum in time to longer integrations
    # (this summing does cause more decorrelation but it's better than without
    # fringe stopping and it decreases data rates)
    fringe_stopped_int_time = ((max_resolution.to(units.arcminute) * max_decorr) /
                               (np.sin(corr_FoV.radian) * earth_rot_speed * abs(np.sin(hera_latitude))))
    params['fringe_stopped_int_time'] = fringe_stopped_int_time
    if verbose:
        print_desc_val('Max integration time to keep the decorrelation due to time '
                       'integrating after fringe stopping under max_decorr for a {fov} '
                       'degree correlator FoV on the longest baselines:'.
                       format(fov=corr_FoV.degree), fringe_stopped_int_time, wrap_len)

    return params