python类deg()的实例源码

def do_annual_parallax_test(filename):
    """testing functions called by a few unit tests"""
    with open(filename) as data_file:
        lines = data_file.readlines()
    ulens_params = lines[3].split()
    event_params = lines[4].split()
    data = np.loadtxt(filename, dtype=None)
    model = Model({
        't_0':float(ulens_params[1])+2450000., 
        'u_0':float(ulens_params[3]), 
        't_E':float(ulens_params[4]), 
        'pi_E_N':float(ulens_params[5]), 
        'pi_E_E':float(ulens_params[6]) }, 
        coords=SkyCoord(
            event_params[1]+' '+event_params[2], unit=(u.deg, u.deg)))
    model.parameters.t_0_par = float(ulens_params[2])+2450000.

    time = data[:,0]
    dataset = MulensData([time, 20.+time*0., 0.1+time*0.,], add_2450000=True)
    model.set_datasets([dataset])
    model.parallax(satellite=False, earth_orbital=True, topocentric=False)
    return np.testing.assert_almost_equal(
        model.data_magnification[0] / data[:,1], 1.0, decimal=4)

def main():
    parser = argparse.ArgumentParser(
        description="Convert among multiple coordinate formats")
    parser.add_argument("coord", nargs="+")
    args = parser.parse_args()

    ra, dec = parse_coord(args.coord)
    info = (
        "%-14s  %-14s\n" % ("R.A.", "Dec.") +
        "%s--%s\n" % ("-"*14, "-"*14) +
        "%-14.3f  %-+14.3f\n" % (ra.deg, dec.deg) +
        "%-14s  %-14s\n" % (
            ra.to_string(unit=au.hourangle, precision=4),
            dec.to_string(unit=au.deg, alwayssign=True, precision=3)) +
        "%-14s  %-14s\n" % (
            ra.to_string(unit=au.hourangle, sep=":", precision=4),
            dec.to_string(unit=au.deg, alwayssign=True,
                          sep=":", precision=3)) +
        "%-14s  %-14s\n" % (
            ra.to_string(unit=au.hourangle, sep=" ", precision=4),
            dec.to_string(unit=au.deg, alwayssign=True,
                          sep=" ", precision=3))
    )
    print(info)

def test_meta_from_position(igmsp):
    # One match
    meta = igmsp.meta_from_position((0.0019,17.7737), 1*u.arcsec)
    assert len(meta) == 1
    # Blank
    meta2 = igmsp.meta_from_position((10.038604,55.298477), 1*u.arcsec)
    assert meta2 is None
    # Multiple sources (insure rank order)
    meta3 = igmsp.meta_from_position((0.0055,-1.5), 1*u.deg)
    assert len(meta3) == 2
    assert np.isclose(meta3['R'][0],meta3['R'][1])
    # Multiple meta entries (GGG)
    meta4 = igmsp.meta_from_position('001115.23+144601.8', 1*u.arcsec)
    assert len(meta4) == 2
    assert meta4['R'][0] != meta4['R'][1]
    # Multiple but grab closest source
    meta5 = igmsp.meta_from_position((0.0055,-1.5), 1*u.deg, max_match=1)
    assert len(meta5) == 1
    # Groups
    meta = igmsp.meta_from_position((2.813500,14.767200), 20*u.deg, groups=['GGG','HD-LLS_DR1'])
    for group in meta['GROUP'].data:
        assert group in ['GGG', 'HD-LLS_DR1']
    # Physical separation
    meta6 = igmsp.meta_from_position('001115.23+144601.8', 300*u.kpc)
    assert len(meta6) == 2

def test_query_position(igmsp):
    # One match
    _, _, idx = igmsp.qcat.query_position((0.0019,17.7737), 1*u.arcsec)
    assert idx >= 0
    # Blank
    _, _, idx = igmsp.qcat.query_position((10.038604,55.298477), 1*u.arcsec)
    assert len(idx) == 0
    # Multiple (insure rank order)
    icoord = SkyCoord(ra=0.0055, dec=-1.5, unit='deg')
    _, subcat, _ = igmsp.qcat.query_position(icoord, 1*u.deg)
    # Test
    coord = SkyCoord(ra=subcat['RA'], dec=subcat['DEC'], unit='deg')
    sep = icoord.separation(coord)
    isrt = np.argsort(sep)
    assert isrt[0] == 0
    assert isrt[-1] == len(subcat)-1
    # Multiple but grab only 1
    _, _, idxs = igmsp.qcat.query_position(icoord, 1*u.deg, max_match=1)
    assert len(idxs) == 1

def test_cat_from_query_coord(igmsp):
    # Single
    coord = SkyCoord(ra=0.0019, dec=17.7737, unit='deg')
    #
    _, ccat, _ = igmsp.qcat.query_coords(coord)
    assert ccat['IGM_ID'][0] == 0
    # Multiple
    coords = SkyCoord(ra=[0.0028,0.0019], dec=[14.9747,17.7737], unit='deg')
    _, ccat2, _ = igmsp.qcat.query_coords(coords)
    assert len(ccat2) == 2
    assert ccat2['IGM_ID'][0] == 1
    assert ccat2['IGM_ID'][1] == 0
    # One miss
    coords3 = SkyCoord(ra=[9.0028,0.0019], dec=[-14.9747,17.7737], unit='deg')
    _, ccat3, _ = igmsp.qcat.query_coords(coords3)
    assert ccat3['IGM_ID'][0] == -1

def test_spectra_from_coord(igmsp):
    # No match
    specN, metaN = igmsp.spectra_from_coord((0.0019, -17.7737))
    assert specN is None
    assert metaN is None
    # One match
    spec, meta = igmsp.spectra_from_coord((0.0019, 17.7737))
    assert spec.nspec == 1
    assert meta['PLATE'][0] == 6173
    # Multiple matches and spectra
    spec, meta = igmsp.spectra_from_coord('001115.23+144601.8')#, groups=['GGG']) # Not in debug file for BOSS or SDSS
    assert spec.nspec == 2
    # Many matches; takes closest
    spec, meta = igmsp.spectra_from_coord((0.0019, 17.7737), tol=10*u.deg)
    assert spec.nspec == 1
    assert meta['PLATE'][0] == 6173

def setUp(self):
        self.dir = './test_results'
        self.lowcore = create_named_configuration('LOWBD2-CORE')
        os.makedirs(self.dir, exist_ok=True)
        self.times = (numpy.pi / (12.0)) * numpy.linspace(-3.0, 3.0, 7)
        self.frequency = numpy.array([1e8])
        self.channel_bandwidth = numpy.array([1e6])
        self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000')
        self.vis = create_visibility(self.lowcore, self.times, self.frequency,
                                     channel_bandwidth=self.channel_bandwidth,
                                     phasecentre=self.phasecentre, weight=1.0,
                                     polarisation_frame=PolarisationFrame('stokesI'))
        self.vis.data['vis'] *= 0.0

        # Create model
        self.model = create_test_image(cellsize=0.0015, phasecentre=self.vis.phasecentre,
                                       frequency=self.frequency)
        self.model.data[self.model.data > 1.0] = 1.0
        self.vis = predict_2d(self.vis, self.model)
        assert numpy.max(numpy.abs(self.vis.vis)) > 0.0
        export_image_to_fits(self.model, '%s/test_solve_skycomponent_model.fits' % (self.dir))
        self.bigmodel = create_image_from_visibility(self.vis, cellsize=0.0015, npixel=512)

def _checkcomponents(self, dirty, fluxthreshold=5.0, positionthreshold=1.0):
        comps = find_skycomponents(dirty, fwhm=1.0, threshold=fluxthreshold, npixels=5)
        assert len(comps) == len(self.components), "Different number of components found: original %d, recovered %d" % \
                                                   (len(self.components), len(comps))
        cellsize = abs(dirty.wcs.wcs.cdelt[0])
        # Check for agreement between image and DFT
        for comp in comps:
            sflux = sum_visibility(self.componentvis, comp.direction)[0]
            assert abs(comp.flux[0, 0] - sflux[0, 0]) < fluxthreshold, \
                "Fitted and DFT flux differ %s %s" % (comp.flux[0, 0], sflux[0, 0])
            # Check for agreement in direction
            ocomp = find_nearest_component(comp.direction, self.components)
            radiff = abs(comp.direction.ra.deg - ocomp.direction.ra.deg) / cellsize
            assert radiff < positionthreshold, "Component differs in dec %.3f pixels" % radiff
            decdiff = abs(comp.direction.dec.deg - ocomp.direction.dec.deg) / cellsize
            assert decdiff < positionthreshold, "Component differs in dec %.3f pixels" % decdiff

def actualSetup(self, sky_pol_frame='stokesIQUV', data_pol_frame='linear', f=None):
        if f is None:
            f = [100.0, 50.0, -10.0, 40.0]
        f = numpy.array(f)
        self.flux = numpy.array([f, 0.8 * f, 0.6 * f])

        # The phase centre is absolute and the component is specified relative (for now).
        # This means that the component should end up at the position phasecentre+compredirection
        self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
        self.compabsdirection = SkyCoord(ra=+181.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
        self.comp = Skycomponent(direction=self.compabsdirection, frequency=self.frequency, flux=self.flux,
                                 polarisation_frame=PolarisationFrame(sky_pol_frame))
        self.vis = create_blockvisibility(self.lowcore, self.times, self.frequency, phasecentre=self.phasecentre,
                                          channel_bandwidth=self.channel_bandwidth, weight=1.0,
                                          polarisation_frame=PolarisationFrame(data_pol_frame))
        self.vis = predict_skycomponent_blockvisibility(self.vis, self.comp)

def setUp(self):
        self.dir = './test_results'
        os.makedirs(self.dir, exist_ok=True)

        self.vnchan = 5
        self.lowcore = create_named_configuration('LOWBD2-CORE')
        self.times = (numpy.pi / 12.0) * numpy.linspace(-3.0, 3.0, 7)
        self.frequency = numpy.linspace(8e7, 1.2e8, self.vnchan)
        self.startfrequency = numpy.array([8e7])
        self.channel_bandwidth = numpy.array(self.vnchan * [self.frequency[1] - self.frequency[0]])
        self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000')
        self.vis = create_visibility(self.lowcore, times=self.times, frequency=self.frequency,
                                     phasecentre=self.phasecentre, weight=1.0,
                                     polarisation_frame=PolarisationFrame('stokesI'),
                                     channel_bandwidth=self.channel_bandwidth)
        self.model = create_image_from_visibility(self.vis, npixel=512, cellsize=0.001,
                                                  nchan=self.vnchan,
                                                  frequency=self.startfrequency)

def setUp(self):
        self.lowcore = create_named_configuration('LOWBD2-CORE')
        self.times = (numpy.pi / 43200.0) * numpy.arange(0.0, 300.0, 30.0)
        self.frequency = numpy.linspace(1.0e8, 1.1e8, 3)
        self.channel_bandwidth = numpy.array([1e7, 1e7, 1e7])
        # Define the component and give it some spectral behaviour
        f = numpy.array([100.0, 20.0, -10.0, 1.0])
        self.flux = numpy.array([f, 0.8 * f, 0.6 * f])

        # The phase centre is absolute and the component is specified relative (for now).
        # This means that the component should end up at the position phasecentre+compredirection
        self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
        self.compabsdirection = SkyCoord(ra=+181.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
        pcof = self.phasecentre.skyoffset_frame()
        self.compreldirection = self.compabsdirection.transform_to(pcof)
        self.comp = Skycomponent(direction=self.compreldirection, frequency=self.frequency, flux=self.flux)

def test_phase_rotation_identity(self):
        self.vis = create_visibility(self.lowcore, self.times, self.frequency,
                                     channel_bandwidth=self.channel_bandwidth,
                                     phasecentre=self.phasecentre, weight=1.0,
                                     polarisation_frame=PolarisationFrame("stokesIQUV"))
        self.vismodel = predict_skycomponent_visibility(self.vis, self.comp)
        newphasecenters = [SkyCoord(182, -35, unit=u.deg), SkyCoord(182, -30, unit=u.deg),
                           SkyCoord(177, -30, unit=u.deg), SkyCoord(176, -35, unit=u.deg),
                           SkyCoord(216, -35, unit=u.deg), SkyCoord(180, -70, unit=u.deg)]
        for newphasecentre in newphasecenters:
            # Phase rotating back should not make a difference
            original_vis = self.vismodel.vis
            original_uvw = self.vismodel.uvw
            rotatedvis = phaserotate_visibility(phaserotate_visibility(self.vismodel, newphasecentre, tangent=False),
                                                self.phasecentre, tangent=False)
            assert_allclose(rotatedvis.uvw, original_uvw, rtol=1e-7)
            assert_allclose(rotatedvis.vis, original_vis, rtol=1e-7)

def setUp(self):
        self.dir = './test_results'
        self.lowcore = create_named_configuration('LOWBD2-CORE')
        os.makedirs(self.dir, exist_ok=True)
        self.times = (numpy.pi / (12.0)) * numpy.linspace(-3.0, 3.0, 7)
        self.frequency = numpy.array([1e8])
        self.channel_bandwidth = numpy.array([1e6])
        self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000')
        self.vis = create_visibility(self.lowcore, self.times, self.frequency,
                                     channel_bandwidth=self.channel_bandwidth,
                                     phasecentre=self.phasecentre, weight=1.0,
                                     polarisation_frame=PolarisationFrame('stokesI'))
        self.vis.data['vis'] *= 0.0

        # Create model
        self.model = create_test_image(cellsize=0.0015, phasecentre=self.vis.phasecentre,
                                       frequency=self.frequency)
        self.model.data[self.model.data > 1.0] = 1.0
        self.vis = predict_2d(self.vis, self.model)
        assert numpy.max(numpy.abs(self.vis.vis)) > 0.0
        export_image_to_fits(self.model, '%s/test_solve_skycomponent_model.fits' % (self.dir))
        self.bigmodel = create_image_from_visibility(self.vis, cellsize=0.0015, npixel=512)

def setUp(self):
        self.dir = './test_results'
        self.lowcore = create_named_configuration('LOWBD2-CORE')
        os.makedirs(self.dir, exist_ok=True)
        self.times = (numpy.pi / (12.0)) * numpy.linspace(-3.0, 3.0, 7)
        self.nchan=8
        self.frequency = numpy.linspace(0.8e8, 1.2e8, self.nchan)
        self.channel_bandwidth = numpy.array(self.nchan*[self.frequency[1]-self.frequency[0]])
        self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000')
        self.vis = create_visibility(self.lowcore, self.times, self.frequency,
                                     channel_bandwidth=self.channel_bandwidth,
                                     phasecentre=self.phasecentre, weight=1.0,
                                     polarisation_frame=PolarisationFrame('stokesI'))
        self.vis.data['vis'] *= 0.0

        # Create model
        self.model = create_test_image(cellsize=0.0015, phasecentre=self.vis.phasecentre,
                                       frequency=self.frequency)
        self.model.data[self.model.data > 1.0] = 1.0
        self.vis = predict_2d(self.vis, self.model)
        assert numpy.max(numpy.abs(self.vis.vis)) > 0.0
        export_image_to_fits(self.model, '%s/test_solve_image_mm_model.fits' % (self.dir))
        self.bigmodel = create_image_from_visibility(self.vis, cellsize=0.0015, npixel=512,
                                                     frequency=self.frequency)

def random_points(self, n=1):
        """
        Generate uniformly distributed random points within the
        sky image (coverage).

        Parameters
        ----------
        n : int, optional
            The number of random points required.
            Default: 1

        Returns
        -------
        lon, lat : float, or 1D `~numpy.ndarray`
            The longitudes and latitudes (in world coordinate)
            generated.
            Unit: [deg]
        """
        raise NotImplementedError

def world2pix(self, x, y):
        """
        Convert the world coordinates (R.A., Dec.) into the pixel
        coordinates (indexes) within the sky data array.

        Parameters
        ----------
        x, y : float, `~numpy.ndarray`
            The R.A., Dec. world coordinates
            Unit: [deg]

        Returns
        -------
        ri, ci : int, `~numpy.ndarray`
            The row, column indexes within the sky data array.
        """
        pixelsize = self.pixelsize * AUC.arcsec2deg  # [deg]
        x, y = np.asarray(x), np.asarray(y)  # [deg]
        ri0, ci0 = self.ysize//2, self.xsize//2
        ri = np.round((y - self.ycenter) / pixelsize + ri0).astype(int)
        ci = np.round((x - self.xcenter) / pixelsize + ci0).astype(int)
        return (ri, ci)

def random_points(self, n=1):
        """
        Generate uniformly distributed random points within the sky patch.

        Returns
        -------
        lon : float, or 1D `~numpy.ndarray`
            Longitudes (Galactic/equatorial);
            Unit: [deg]
        lat : float, or 1D `~numpy.ndarray`
            Latitudes (Galactic/equatorial);
            Unit: [deg]
        """
        lon_min, lon_max = self.lon_limit
        lat_min, lat_max = self.lat_limit
        lon = np.random.uniform(low=lon_min, high=lon_max, size=n)
        lat = np.random.uniform(low=lat_min, high=lat_max, size=n)
        return (lon, lat)

def signal_gaussian(
            signal_location=np.array([61, 21])*u.deg,
            fov_center=np.array([60, 20])*u.deg,
            width=0.05*u.deg,
            signal_events=80,
            bins=[80, 80],
            fov=4*u.deg,
        ):

    # reshape so if signal_events = 1 the array can be indexed in the same way.
    signal = multivariate_normal.rvs(
                mean=signal_location.value,
                cov=width.value,
                size=signal_events
                ).reshape(signal_events, 2)
    r = np.array([fov_center - fov/2, fov_center + fov/2]).T

    signal_hist, _, _ = np.histogram2d(signal[:, 0], signal[:, 1], bins=bins, range=r)
    return signal_hist

def prepare(self, path, star_finder):
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', AstropyWarning)
            self.calibration.selectImage(path)

        self.names, self.vmag, alt, az = self.calibration.catalog.filter(Configuration.min_alt * u.deg, Configuration.max_mag)

        altaz = np.array([alt.radian, az.radian]).transpose()

        self.pos = np.column_stack(self.calibration.project(altaz))

        self.finder = star_finder

        self.image = cv2.imread(path)

        self.finder.setImage(self.image)

        self.altaz = np.array([alt.degree, az.degree]).transpose()

def galToCel(ll, bb):
    """
    Converts Galactic (deg) to Celestial J2000 (deg) coordinates
    """
    bb = numpy.radians(bb)
    sin_bb = numpy.sin(bb)
    cos_bb = numpy.cos(bb)

    ll = numpy.radians(ll)
    ra_gp = numpy.radians(192.85948)
    de_gp = numpy.radians(27.12825)
    lcp = numpy.radians(122.932)

    sin_lcp_ll = numpy.sin(lcp - ll)
    cos_lcp_ll = numpy.cos(lcp - ll)

    sin_d = (numpy.sin(de_gp) * sin_bb) \
            + (numpy.cos(de_gp) * cos_bb * cos_lcp_ll)
    ramragp = numpy.arctan2(cos_bb * sin_lcp_ll,
                            (numpy.cos(de_gp) * sin_bb) \
                            - (numpy.sin(de_gp) * cos_bb * cos_lcp_ll))
    dec = numpy.arcsin(sin_d)
    ra = (ramragp + ra_gp + (2. * numpy.pi)) % (2. * numpy.pi)
    return numpy.degrees(ra), numpy.degrees(dec)

def dec2dms(dec):
    """
    ADW: This should really be replaced by astropy
    """
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.

    if isinstance(dec,basestring):
        dec = float(dec)

    sign = numpy.copysign(1.0,dec)

    fdeg = np.abs(dec)
    deg = int(fdeg)

    fminute = (fdeg - deg)*MINUTE
    minute = int(fminute)

    second = (fminute - minute)*MINUTE

    deg = int(deg * sign)
    return (deg, minute, second)

def test_xarray():
    import numpy as np
    import astropy.units as u
    from ..utils import xarray, yarray

    x = xarray((3, 3))
    xx = np.array([[0, 1, 2], [0, 1, 2], [0, 1, 2]])
    assert np.all(x == xx)

    x = xarray((3, 3), yx=(1, 0))
    assert np.all(x == xx)

    x = xarray((3, 3), yx=(0, 0.5))
    xx = xx - 0.5
    assert np.all(x == xx)

    x = xarray((3, 3), th=90 * u.deg)
    y = yarray((3, 3))
    assert np.allclose(x, y)

def test_yarray():
    import numpy as np
    import astropy.units as u
    from ..utils import yarray, xarray

    y = yarray((3, 3))
    yy = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]])
    assert np.all(y == yy)

    y = yarray((3, 3), yx=(0, 1))
    assert np.all(y == yy)

    y = yarray((3, 3), yx=(0.5, 0))
    yy = yy - 0.5
    assert np.all(y == yy)

    y = yarray((3, 3), th=-90 * u.deg)
    x = xarray((3, 3))
    assert np.allclose(y, x)

def match_to_database(source_list, image_list, cat, Ast, imscale, thresh=8.):
    '''
    RA, Decdetection catalogs and vizier (astroquery) match?  vo conesearch?    
    '''
    matches = []
    for sources, image in zip(source_list, image_list):
        ifuslot = sources[1]
        Ast.get_ifuslot_projection(ifuslot, imscale, image.shape[1]/2.,
                                       image.shape[0]/2.)
        for source in sources[0]:
            ra, dec = Ast.tp_ifuslot.wcs.wcs_pix2world(source['xcentroid'], 
                                                       source['ycentroid'],1)
            c = SkyCoord(ra, dec, unit=(unit.degree, unit.degree))
            idx, d2d, d3d = match_coordinates_sky(c, cat, nthneighbor=1) 
            if d2d.arcsec[0] < thresh:
                dra =  cat.ra[idx].deg - float(ra)
                ddec = cat.dec[idx].deg - float(dec)
                matches.append([int(ifuslot),int(idx), 
                                float(ra), float(dec), 
                                cat.ra[idx].deg, 
                                cat.dec[idx].deg, dra, ddec,
                                d2d.arcsec[0]])
    return matches

def get_wise(ra,dec,band):
    mission='wise'
    dataset='allwise'
    table='p3am_cdd'
    successful=False
    while not successful:
        try:
            results=Ibe.query_region(coord.SkyCoord(ra, dec, unit=(u.deg, u.deg), frame='icrs'),mission=mission,dataset=dataset,table=table)
            successful=True
        except requests.exceptions.ConnectionError:
            print 'Connection failed, retrying'
            sleep(10)

    url = 'http://irsa.ipac.caltech.edu/ibe/data/'+mission+'/'+dataset+'/'+table+'/'
    params = { 'coadd_id': results[results['band']==band]['coadd_id'][0],
           'band': band }
    params['coaddgrp'] = params['coadd_id'][:2]
    params['coadd_ra'] = params['coadd_id'][:4]
    path = str.format('{coaddgrp:s}/{coadd_ra:s}/{coadd_id:s}/{coadd_id:s}-w{band:1d}-int-3.fits',**params)
    outname=path.split('/')[-1]
    download_file(url+path,outname)
    return outname

def peek(self, figsize=(10, 10), color='b', alpha=0.75,
             print_to_file=None, **kwargs):
        """
        Show extracted fieldlines overlaid on HMI image.
        """
        fig = plt.figure(figsize=figsize)
        ax = fig.gca(projection=self.hmi_map)
        self.hmi_map.plot()
        ax.set_autoscale_on(False)
        for stream, _ in self.streamlines:
            ax.plot(self._convert_angle_to_length(stream[:, 0]*u.cm,
                                                  working_units=u.arcsec).to(u.deg),
                    self._convert_angle_to_length(stream[:, 1]*u.cm,
                                                  working_units=u.arcsec).to(u.deg),
                    alpha=alpha, color=color, transform=ax.get_transform('world'))

        if print_to_file is not None:
            plt.savefig(print_to_file, **kwargs)
        plt.show()

def coords(self):
        """Coordinates of event"""
        coords = SkyCoord(
            self.event_params[1]+' '+self.event_params[2], 
            unit=(u.deg, u.deg))
        return coords

def test_model_coords():
    coords = SkyCoord('18:00:00 -30:00:00', unit=(u.hourangle, u.deg))
    model_1 = Model({'t_0':2450000, 'u_0':0.1, 't_E':100}, coords='18:00:00 -30:00:00')
    assert isinstance(model_1.coords, SkyCoord)
    assert model_1.coords.ra == coords.ra
    assert model_1.coords.dec == coords.dec
    assert  model_1.coords.dec.deg == -30.00

    model_3 = Model({'t_0':2450000, 'u_0':0.1, 't_E':100})
    model_3.coords = '17:00:00 -27:32:14'
    assert model_3.coords.to_string('hmsdms') == '17h00m00s -27d32m14s'

def test_data_coords():
    coords = SkyCoord('18:00:00 -30:00:00', unit=(u.hourangle, u.deg))
    data_1 = MulensData(
        file_name=SAMPLE_FILE_01,
        coords='18:00:00 -30:00:00')
    assert isinstance(data_1.coords, SkyCoord)
    assert data_1.coords.ra == coords.ra
    assert data_1.coords.dec == coords.dec
    assert data_1.coords.dec.deg == -30.00

    data_3 = MulensData(file_name=SAMPLE_FILE_01)
    data_3.coords = '17:00:00 -27:32:14'
    assert data_3.coords.to_string('hmsdms') == '17h00m00s -27d32m14s'

def test_BLPS_01():
    """simple binary lens with point source"""
    params = ModelParameters({
            't_0':6141.593, 'u_0':0.5425, 't_E':62.63*u.day, 
            'alpha':49.58*u.deg, 's':1.3500, 'q':0.00578})
    model = Model(parameters=params)
    t = np.array([6112.5])
    data = MulensData(data_list=[t, t*0.+16., t*0.+0.01])
    model.set_datasets([data])
    magnification = model.data_magnification[0][0]
    np.testing.assert_almost_equal(magnification, 4.691830781584699) # This value comes from early version of this code.
    # np.testing.assert_almost_equal(m, 4.710563917) # This value comes from Andy's getbinp().