def xyz2eq(xin,yin,zin, units='deg'):
    """
    Convert x,y,z on the unit sphere to RA DEC.

    parameters
    ----------
    x,y,z: 
        scalars or arrays as given by eq2xyz
    units: string, optional
        'deg' if the output is to be degrees, 'rad' if it is to be radians.
        Default is degrees.

    Notes:
        This follows the same convention as the STOMP package.
    """

    x = numpy.array(xin, ndmin=1, copy=False)
    y = numpy.array(yin, ndmin=1, copy=False)
    z = numpy.array(zin, ndmin=1, copy=False)

    theta,phi = _xyz2thetaphi(x,y,z)
    theta += _sdsspar['node']

    if units == 'deg':
        numpy.rad2deg(theta,theta)
        numpy.rad2deg(phi,phi)

    atbound(theta, 0.0, 360.0)

    # theta->ra, phi->dec
    return theta,phi

    def load_star_data(self, name):
        
#         stetson_df = pandas.read_pickle('stetson.pd')
        files = glob.glob('cam1/*.fits')
        for thisFile in files:
            thisFits = pf.open(thisFile)
            fibreTable = thisFits['FIBRES'].data
            idx = fibreTable.field('NAME').strip()==name
            if fibreTable[idx].shape[0] >0: #if there is any star in this file...
                starInfo = fibreTable[idx][0] 

                self.name = starInfo.field('NAME').strip()
                self.RA = np.rad2deg(starInfo.field('RA'))
                self.Dec = np.rad2deg(starInfo.field('DEC'))        
                self.RA_deg, self.RA_min, self.RA_sec = toolbox.dec2sex(self.RA)   
                self.Dec_deg, self.Dec_min, self.Dec_sec = toolbox.dec2sex(self.Dec)
                self.Vmag = starInfo.field('MAGNITUDE')
            
#                 self.B = stetson_df[stetson_df.target == self.name].B.values[0]
#                 self.I = stetson_df[stetson_df.target == self.name].I.values[0]
#                 self.R = stetson_df[stetson_df.target == self.name].R.values[0]
#                 self.U = stetson_df[stetson_df.target == self.name].U.values[0]
#                 self.V = stetson_df[stetson_df.target == self.name].mag.values[0]
#                 self.BV = stetson_df[stetson_df.target == self.name].BV.values[0]
                print self.name,'star created'
                break

def plot_poses_vels_theta_omega(pylab, poses, vels):
    thetas = np.array([angle_from_SE2(p) for p in poses])
    omegas = np.array([linear_angular_from_se2(vel)[1] for vel in vels])
    positive = omegas > 0
    negative = np.logical_not(positive)
    pylab.plot(np.rad2deg(thetas)[positive], omegas[positive], 'r.')
    pylab.plot(np.rad2deg(thetas)[negative], omegas[negative], 'b.')

    def __reverse_lambert(self,x,y):

        p     = math.pi
        ln    = np.log
        power = np.power

        sin = np.sin
        cos = np.cos
        tan = np.tan
        cot = lambda x: 1./tan(x)
        sec = lambda x: 1./cos(x)

        ra0  = np.deg2rad(self.center[0])
        dec0 = np.deg2rad(self.center[1])

        dec1 = np.deg2rad(self.ref_dec1)
        dec2 = np.deg2rad(self.ref_dec2)

        n = ln(cos(dec1)*sec(dec2))/ln(tan(p/4+dec2/2)*cot(p/4+dec1/2))
        F = cos(dec1)*power(tan(p/4+dec1/2),n)/n
        rho0 = F*power(cos(p/4+dec0/2)/sin(p/4+dec0/2),n)
        
        rho = np.sign(n)*np.sqrt(x**2+(rho0-y)**2)
        theta = np.arctan(x/(rho0-y))
        dec = 2.*np.arctan(np.power((F/rho),(1./n))) - p/2
        ra = ra0 + theta/n
        return np.rad2deg(ra), np.rad2deg(dec)

    def test_degrees(self):

        # the following doesn't make much sense in terms of the name of the
        # routine, but we check it gives the correct result.
        q1 = np.rad2deg(60. * u.degree)
        assert_allclose(q1.value, 60.)
        assert q1.unit == u.degree

        q2 = np.degrees(60. * u.degree)
        assert_allclose(q2.value, 60.)
        assert q2.unit == u.degree

        q3 = np.rad2deg(np.pi * u.radian)
        assert_allclose(q3.value, 180.)
        assert q3.unit == u.degree

        q4 = np.degrees(np.pi * u.radian)
        assert_allclose(q4.value, 180.)
        assert q4.unit == u.degree

        with pytest.raises(TypeError):
            np.rad2deg(3. * u.m)

        with pytest.raises(TypeError):
            np.degrees(3. * u.m)

 def __call__(self, phi, theta):
     phi = np.deg2rad(phi)
     theta = np.deg2rad(theta)
     rtheta = self._compute_rtheta(theta)
     x = np.rad2deg(rtheta * np.sin(phi))
     y = - np.rad2deg(rtheta * np.cos(phi))
     return x, y

def main():
    if args["<date>"]:
        date = args["<date>"]
        time = args["<timeutc>"]
    else:
        date, time = enter_datetime()

    valid = False
    while not valid:
        try:
            date = datetime.strptime(date + " " + time, "%Y-%m-%d %H:%M")
            valid = True
        except ValueError:
            print("Could not parse date/time, please use the given notation\n")
            date, time = enter_datetime()

    fact = fact_setup(date)
    stars = create_dataframe(fact)
    az, alt, ra, dec = dark_spot_gridsearch(fact, stars, min_altitude, 0.25, 2)
    stars_fov = get_stars_in_fov(az, alt, stars, fact)

    print(u"best ratescan position:")
    print(u"RA: {:1.3f} h".format(np.rad2deg(ra) * 24 / 360))
    print(u"DEC: {:1.3f}°".format(np.rad2deg(dec)))
    print(u"Az: {:1.3f}°".format(np.rad2deg(az)))
    print(u"Alt: {:1.3f}°".format(np.rad2deg(alt)))
    print(u"Brightest star in FOV: {:1.2f} mag".format(stars_fov.vmag.min()))

    if args["--plot"]:
        plot_dark_spot(stars, az, alt, min_altitude)

def get_lonlatalt(pos, utc_time):
    """Calculate sublon, sublat and altitude of satellite, considering the
    earth an ellipsoid.

    http://celestrak.com/columns/v02n03/
    """
    (pos_x, pos_y, pos_z) = pos / XKMPER
    lon = ((np.arctan2(pos_y * XKMPER, pos_x * XKMPER) - astronomy.gmst(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-10):
            break
    alt = r / np.cos(lat) - c
    alt *= A
    return np.rad2deg(lon), np.rad2deg(lat), alt

 def interpolate_trajectory(self, t, return_degrees=True, **kw):
     lat = self.interpolate_timeseries('latitude', t, **kw)
     lon = self.interpolate_timeseries('longitude', t, **kw)
     dep = self.interpolate_timeseries('depth', t, **kw)
     if return_degrees and self.units('latitude') is 'radian':
         lat, lon = np.rad2deg(lat), np.rad2deg(lon)
     return lon, lat, dep

    def evaluate(cls, x, y, mu, gamma):
        phi = np.arctan2(x / np.cos(gamma), -y)
        r = cls._compute_r_theta(x, y, gamma)
        pho = r / (cls.r0 * (mu + 1) + y * np.sin(gamma))
        psi = np.arctan2(1, pho)
        omega = np.arcsin((pho * mu) / np.sqrt(pho ** 2 + 1))

        theta1 = np.rad2deg(psi - omega)
        theta2 = np.rad2deg(psi + omega) + 180

        if np.abs(mu) < 1:
            if theta1 < 90 and theta1 > -90:
                theta = theta1
            else:
                theta = theta2
        else:
            # theta1dif = 90 - theta1
            # theta2dif = 90 - theta2
            if theta1 < theta2:
                theta = theta1
            else:
                theta = theta2

        phi = np.rad2deg(phi)
        return phi, theta

    def getJointLimits(self, input_urdf, use_deg=True):
        import xml.etree.ElementTree as ET
        tree = ET.parse(input_urdf)
        limits = {}
        for j in tree.findall('joint'):
            name = j.attrib['name']
            torque = 0
            lower = 0
            upper = 0
            velocity = 0
            if j.attrib['type'] == 'revolute':
                l = j.find('limit')
                if l != None:
                    torque = l.attrib['effort']   #this is not really the physical limit but controller limit, but well
                    lower = l.attrib['lower']
                    upper = l.attrib['upper']
                    velocity = l.attrib['velocity']

                    limits[name] = {}
                    limits[name]['torque'] = float(torque)
                    if use_deg:
                        limits[name]['lower'] = np.rad2deg(float(lower))
                        limits[name]['upper'] = np.rad2deg(float(upper))
                        limits[name]['velocity'] = np.rad2deg(float(velocity))
                    else:
                        limits[name]['lower'] = float(lower)
                        limits[name]['upper'] = float(upper)
                        limits[name]['velocity'] = float(velocity)
        return limits

def compute_fiducial(wcslist, domain=None):
    """
    For a celestial footprint this is the center.
    For a spectral footprint, it is the beginning of the range.

    This function assumes all WCSs have the same output coordinate frame.
    """
    output_frame = getattr(wcslist[0], 'output_frame')
    axes_types = getattr(wcslist[0], output_frame).axes_type
    spatial_axes = np.array(axes_types) == 'SPATIAL'
    spectral_axes = np.array(axes_types) == 'SPECTRAL'
    footprints = np.hstack([w.footprint(domain=domain) for w in wcslist])
    spatial_footprint = footprints[spatial_axes]
    spectral_footprint = footprints[spectral_axes]

    fiducial = np.empty(len(axes_types))
    if (spatial_footprint).any():
        lon, lat = spatial_footprint
        lon, lat = np.deg2rad(lon), np.deg2rad(lat)
        x_mean = np.mean(np.cos(lat) * np.cos(lon))
        y_mean = np.mean(np.cos(lat) * np.sin(lon))
        z_mean = np.mean(np.sin(lat))
        lon_fiducial = np.rad2deg(np.arctan2(y_mean, x_mean)) % 360.0
        lat_fiducial = np.rad2deg(np.arctan2(z_mean, np.sqrt(x_mean**2 + y_mean\
** 2)))
        fiducial[spatial_axes] = lon_fiducial, lat_fiducial
    if (spectral_footprint).any():
        fiducial[spectral_axes] = spectral_footprint.min()
    return fiducial

def check_visible(src,aa,t=None,check27m=True,check2m=True):
    ''' Check whether a source is observable from the location of aa
        given the limits on az/el of the 2-m's and HA/dec of the 27-m's, at
        date/time dt. Returns True if visible, False otherwise.
        
        If no time is given, it checks for the current time.
        
        By default, returns True only if it is visible from both the 2-m's and the
        27-m's; use check27m=False or check2m=False to ignore either the 27-m's or
        the 2-m's, respectively.
    '''
    
    if t is None:
        t = util.Time.now()
    
    # get alt, az, ha, dec at given date/time    
    aa.set_jultime(t.jd)
    src.compute(aa)
    alt = rad2deg(src.alt)
    az  = rad2deg(src.az)
    ha  = rad2deg(eovsa_ha(src,t))
    dec = rad2deg(src.dec)
    
    if check27m and check2m:
        return check_27m_visible(ha,dec) and check_2m_visible(az,alt) and check_2meq_visible(ha,dec)
    elif check27m:
        return check_27m_visible(ha,dec)
    elif check2m:
        return check_2m_visible(az,alt) and check_2meq_visible(ha,dec)
    else:
        return True

def sphdist(ra1, dec1, ra2, dec2, units=['deg','deg']):
    """
    Get the arc length between two points on the unit sphere

    parameters
    ----------
    ra1,dec1,ra2,dec2: scalar or array
        Coordinates of two points or sets of points. 
        Must be the same length.
    units: sequence
        A sequence containing the units of the input and output.  Default
        ['deg',deg'], which means inputs and outputs are in degrees.  Units
        can be 'deg' or 'rad'
    """
    
    units_in,units_out = units

    # note x,y,z from eq2xyz always returns 8-byte float
    x1,y1,z1 = eq2xyz(ra1, dec1, units=units_in)
    x2,y2,z2 = eq2xyz(ra2, dec2, units=units_in)

    costheta = x1*x2 + y1*y2 + z1*z2
    costheta.clip(-1.0,1.0,out=costheta)

    theta = arccos(costheta)

    if units_out == 'deg':
        numpy.rad2deg(theta,theta)
    return theta

    def evaluate(cls, alpha_C, delta_C, phi, theta, psi):
        """
        Evaluate ZXZ rotation into native coordinates.

        This is like RotateNative2Celestial.evaluate except phi and psi are
        swapped in ZXZ rotation.

        Note currently the parameter values are passed in as degrees so we need
        to convert all inputs to radians and all outputs back to degrees.
        """

        alpha_C = np.deg2rad(alpha_C)
        delta_C = np.deg2rad(delta_C)
        phi = np.deg2rad(phi)
        theta = np.deg2rad(theta)
        psi = np.deg2rad(psi)

        phi_N, theta_N = cls._rotate_zxz(alpha_C, delta_C, psi, theta, phi)

        phi_N = np.rad2deg(phi_N)
        theta_N = np.rad2deg(theta_N)

        mask = phi_N > 180
        if isinstance(mask, np.ndarray):
            phi_N[mask] -= 360
        elif mask:
            phi_N -= 360

        return phi_N, theta_N

    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) = self.get_position(
            utc_time, normalize=True)

        lon = ((np.arctan2(pos_y * XKMPER, pos_x * XKMPER) - astronomy.gmst(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-10):
                break
        alt = r / np.cos(lat) - c
        alt *= A
        return np.rad2deg(lon), np.rad2deg(lat), alt

def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt):
    """Calculate observers look angle to a satellite.
    http://celestrak.com/columns/v02n02/

    utc_time: Observation time (datetime object)
    lon: Longitude of observer position on ground in degrees east
    lat: Latitude of observer position on ground in degrees north
    alt: Altitude above sea-level (geoid) of observer position on ground in km

    Return: (Azimuth, Elevation)
    """
    (pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = astronomy.observer_position(
        utc_time, sat_lon, sat_lat, sat_alt)

    (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \
        astronomy.observer_position(utc_time, lon, lat, alt)

    lon = np.deg2rad(lon)
    lat = np.deg2rad(lat)

    theta = (astronomy.gmst(utc_time) + lon) % (2 * np.pi)

    rx = pos_x - opos_x
    ry = pos_y - opos_y
    rz = pos_z - opos_z

    sin_lat = np.sin(lat)
    cos_lat = np.cos(lat)
    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)

    top_s = sin_lat * cos_theta * rx + \
        sin_lat * sin_theta * ry - cos_lat * rz
    top_e = -sin_theta * rx + cos_theta * ry
    top_z = cos_lat * cos_theta * rx + \
        cos_lat * sin_theta * ry + sin_lat * rz

    az_ = np.arctan(-top_e / top_s)

    if has_xarray and isinstance(az_, xr.DataArray):
        az_data = az_.data
    else:
        az_data = az_

    if has_dask and isinstance(az_data, da.Array):
        az_data = da.where(top_s > 0, az_data + np.pi, az_data)
        az_data = da.where(az_data < 0, az_data + 2 * np.pi, az_data)
    else:
        az_data[np.where(top_s > 0)] += np.pi
        az_data[np.where(az_data < 0)] += 2 * np.pi

    if has_xarray and isinstance(az_, xr.DataArray):
        az_.data = az_data
    else:
        az_ = az_data

    rg_ = np.sqrt(rx * rx + ry * ry + rz * rz)
    el_ = np.arcsin(top_z / rg_)

    return np.rad2deg(az_), np.rad2deg(el_)

def solarcalculator_test():
    import ephem
    import pytz
    from datetime import datetime, timedelta
    import numpy as np
    from solarcalculator import SolarCalculator

    print("Testing SolarCalculator interface.")

    sc = SolarCalculator(("65:01", "25:28")) # Passing coordinates of Oulu to constructor
    assert sc is not None, 'SolarCalculator initialization returned None'
    print("  Initialization ok.")

    hel = ephem.city('Helsinki')
    sc.observer = hel
    assert (int(sc.observer.lat/3.14159*180*1e5)/1e5) == 60.16986, 'Observer coordinates (lat) incorrect: {}'.format(sc.observer.lat/3.14159*180)
    assert (int(sc.observer.long/3.14159*180*1e5)/1e5) == 24.93826, 'Observer coordinates (long) incorrect: {}'.format(sc.observer.long/3.14159*180)
    print("  Observer setting ok.")

    sc.coords = ('37.8044', '-122.2697', 3.0)
    assert np.round(np.rad2deg(sc.observer.lat),4) == 37.8044, 'Changing observer (lat) coordinates failed: {}'.format(np.rad2deg(sc.observer.lat))
    assert np.round(np.rad2deg(sc.observer.long),4) == -122.2697, 'Changing observer (long) coordinates failed: {}'.format(np.rad2deg(sc.observer.long))
    print("  Coordinate setting ok.")

    time_str = '2012/11/16 20:34:56'
    time_fmt = '%Y/%m/%d %H:%M:%S'
    dt = datetime.strptime(time_str, time_fmt).replace(tzinfo=pytz.utc)
    assert dt.strftime(time_fmt) == time_str, 'Initializing a datetime object failed: {}'.format(dt.strftime(time_fmt))
    print("  Initialized datetime object ok.")
    sc.date = dt
    assert sc.localized_date(tz=pytz.utc).strftime(time_fmt) == time_str, 'Initializing date to solarcalculator failed: {}'.format(sc.localized_date(tz=pytz.utc).strftime(time_fmt))
    print("  Changed observer date ok.")
    tz = sc.timezone
    assert str(tz) == 'America/Los_Angeles', 'Timezone is incorrect: {}'.format(tz)
    print("  Timezone found ok.")

    assert sc.date.strftime(time_fmt) == '2012/11/16 12:34:56', 'Observer date did not return a proper result: {}'.format(sc.date.strftime(time_fmt))
    print("  Conversion to local timezone ok.")

    sdt = sc.solartime
    assert sdt.strftime(time_fmt) == '2012/11/16 12:40:55', 'Solartime did not return proper result: {}'.format(sdt.strftime(time_fmt))
    assert sdt.strftime("%z") == '-0754', 'Timezone offset incorrect: {}'.format(sdt.strftime("%z"))

    obs = ephem.city('Helsinki')
    sc.observer = obs
    assert sc.coords == (1.0501613384326405, 0.43525439939787997), 'Coordinates were not set properly: {}'.format(sc.coords)

    sc.coords = ('60:10:11.3', '24.9')

    dt = pytz.timezone('Europe/Stockholm').localize(datetime(2016, 1, 6, 20, 41, 31, 0))
    sc.date = dt
    assert sc.date.strftime(time_fmt) == '2016/01/06 21:41:31', 'Setting time from wrong local timezone failed: {}'.format(sc.date.strftime(time_fmt))
    print("  Date and time changed from wrongly localized time ok.")

    assert sc.solartime.strftime(time_fmt) == '2016/01/06 21:15:22', 'Calculating solar time in east failed: {}'.format(sc.solartime.strftime(time_fmt))
    print("  Solar time calculated ok.")

    tz = sc.timezone
    assert str(tz) == 'Europe/Helsinki', 'Timezone location failed: {}'.format(tz)
    print("  Timezone located ok.")

def calc_geoinc(trace,metric=True):
    """docstring for calc_geoinc"""
    stla=trace.stats.sac.stla
    stlo=trace.stats.sac.stlo
    stdp=trace.stats.sac.stdp
    evla=trace.stats.sac.evla
    evlo=trace.stats.sac.evlo
    evdp=trace.stats.sac.evdp
    
    
    if metric:
        baz=np.rad2deg(np.arctan2((evlo-stlo),(evla-stla)))
        EpiDist = np.sqrt((evlo-stlo)**2. +  (evla-stla)**2.)
        inc = np.rad2deg(np.arctan(EpiDist/ (evdp-stdp)))
        
        HypoDist = np.sqrt((evdp-stdp)**2. + EpiDist**2)
         
    if baz<0.:
        baz=baz+360.
        
    
    azi=np.mod(baz+180.0,360.0)
    inc=np.mod(inc+180.0,180.)

    
    return azi,inc,baz,HypoDist,EpiDist,stdp

 def calc(self, observer):
   self.compute(observer)
   self.ze = np.pi/2 - self.alt
   self.nsr = nsr(self.az, self.ze)
   self.ewr = ewr(self.az, self.ze)
   self.nsd = np.rad2deg(self.nsr)
   self.ewd = np.rad2deg(self.ewr)