def dl_to_rphi_2d(d,l,degree=False,ro=1.,phio=0.):
    """
    NAME:

       dl_to_rphi_2d

    PURPOSE:

       convert Galactic longitude and distance to Galactocentric radius and azimuth

    INPUT:

       d - distance

       l - Galactic longitude [rad/deg if degree]

    KEYWORDS:

       degree= (False): l is in degrees rather than rad

       ro= (1) Galactocentric radius of the observer

       phio= (0) Galactocentric azimuth of the observer [rad/deg if degree]

    OUTPUT:

       (R,phi); phi in degree if degree

    HISTORY:

       2012-01-04 - Written - Bovy (IAS)

    """
    scalarOut, listOut= False, False
    if isinstance(d,(int,float)):
        d= sc.array([d])
        scalarOut= True
    elif isinstance(d,list):
        d= sc.array(d)
        listOut= True
    if isinstance(l,(int,float)):
        l= sc.array([l])
    elif isinstance(l,list):
        l= sc.array(l)
    if degree:
        l*= _DEGTORAD
    R= sc.sqrt(ro**2.+d**2.-2.*d*ro*sc.cos(l))
    phi= sc.arcsin(d/R*sc.sin(l))
    indx= (ro/sc.cos(l) < d)*(sc.cos(l) > 0.)
    phi[indx]= sc.pi-sc.arcsin(d[indx]/R[indx]*sc.sin(l[indx]))
    if degree:
        phi/= _DEGTORAD
    phi+= phio
    if scalarOut:
        return (R[0],phi[0])
    elif listOut:
        return (list(R),list(phi))
    else:
        return (R,phi)

def rphi_to_dl_2d(R,phi,degree=False,ro=1.,phio=0.):
    """
    NAME:

       rphi_to_dl_2d

    PURPOSE:

       convert Galactocentric radius and azimuth to distance and Galactic longitude

    INPUT:

       R - Galactocentric radius

       phi - Galactocentric azimuth [rad/deg if degree]

    KEYWORDS:

       degree= (False): phi is in degrees rather than rad

       ro= (1) Galactocentric radius of the observer

       phio= (0) Galactocentric azimuth of the observer [rad/deg if degree]

    OUTPUT:

       (d,l); phi in degree if degree

    HISTORY:

       2012-01-04 - Written - Bovy (IAS)

    """
    scalarOut, listOut= False, False
    if isinstance(R,(int,float)):
        R= sc.array([R])
        scalarOut= True
    elif isinstance(R,list):
        R= sc.array(R)
        listOut= True
    if isinstance(phi,(int,float)):
        phi= sc.array([phi])
    elif isinstance(phi,list):
        phi= sc.array(phi)
    phi-= phio
    if degree:
        phi*= _DEGTORAD
    d= sc.sqrt(R**2.+ro**2.-2.*R*ro*sc.cos(phi))
    l= sc.arcsin(R/d*sc.sin(phi))
    indx= (ro/sc.cos(phi) < R)*(sc.cos(phi) > 0.)
    l[indx]= sc.pi-sc.arcsin(R[indx]/d[indx]*sc.sin(phi[indx]))
    if degree:
        l/= _DEGTORAD
    if scalarOut:
        return (d[0],l[0])
    elif listOut:
        return (list(d),list(l))
    else:
        return (d,l)

def get_angle (sin_in, cos_in):
    """ Just use arctan2 - much more efficient"""
    if (sin_in >= 0.0 and cos_in >= 0.0):     out = sc.arcsin(sin_in)
    elif (sin_in >= 0.0 and cos_in < 0.0):    out = ma.pi-sc.arcsin(sin_in)
    elif (sin_in < 0.0 and cos_in < 0.0):     out = ma.pi-sc.arcsin(sin_in)
    elif (sin_in < 0.0 and cos_in >= 0.0):    out = 2.0*ma.pi + sc.arcsin(sin_in)
    else: out = float('nan')
    return out

def matrixToEuler(m,order='Aerospace',inDegrees=True):
    if order == 'Aerospace' or order == 'ZYX':
        sp = -m[2,0]
        if sp < (1-EPS):
            if sp > (-1+EPS):
                p = arcsin(sp)
                r = arctan2(m[2,1],m[2,2])
                y = arctan2(m[1,0],m[0,0])
            else:
                p = -pi/2.
                r = 0
                y = pi-arctan2(-m[0,1],m[0,2])
        else:
            p = pi/2.
            y = arctan2(-m[0,1],m[0,2])
            r = 0
        
        if inDegrees:
            return degrees((y,p,r))
        else:
            return (y,p,r)
    elif order == 'BVH' or order == 'ZXY':
        sx = m[2,1]
        if sx < (1-EPS):
            if sx > (-1+EPS):
                x = arcsin(sx)
                z = arctan2(-m[0,1],m[1,1])
                y = arctan2(-m[2,0],m[2,2])
            else:
                x = -pi/2
                y = 0
                z = -arctan2(m[0,2],m[0,0])
        else:
            x = pi/2
            y = 0
            z = arctan2(m[0,2],m[0,0])
        if inDegrees:
            return degrees((z,x,y))
        else:
            return (z,x,y)

    elif order == "ZXZ":
        x = arccos(m[2,2])
        z2 = arctan2(m[2,0],m[2,1])
        z1 = arctan2(m[0,2],-m[1,2])
        if inDegrees:
            return degrees((z1,x,z2))
        else:
            return (z1,x,z2)

def rayshape_old(self, A, d, n):
    # We'll get 4n+4 points, but there are duplicates at the four corners.
    # So, total = 4n
    rad = sp.sqrt(A/sp.pi) # the radius of the base circle
    angle = sp.arcsin(1/(sp.sqrt(2) * (1+d))) # each quadrant isn't -pi/4 < theta < pi/4. A bit less than that.
    base = sp.linspace(-angle, angle, n+1) # split up a quarter of the circumference to n pieces (omitting the point at pi/4)
    C1 = []
    for arg in base:
        x = -d + (rad+d)*sp.cos(arg)
        y = (rad+d)*sp.sin(arg)
        C1.append((x,y))
    # Construct C3, which is the image of C1 by rotation by pi.
    C3 = []
    for pt in C1:
        C3.append((-pt[0] , -pt[1]))
    # Now construct C2
    base = sp.linspace(sp.pi/2 - angle, sp.pi/2 + angle, n+1) # split up a quarter of the circumference to n pieces
    C2 = []
    for arg in base:
        x = (rad+d)*sp.cos(arg)
        y = -d + (rad+d)*sp.sin(arg)
        C2.append((x,y))
    # Construct C4 from C2 by applying rotation by pi.
    C4 = []
    for pt in C2:
        C4.append((-pt[0] , -pt[1]))
 
    return sp.vstack((C1, C2, C3, C4)) 

    def get_bl(self,ra=None,dec=None):
        """
        http://scienceworld.wolfram.com/astronomy/GalacticCoordinates.html
        """
        if ra==None: ra=self.get_column("RA")
        if dec==None: dec=self.get_column("DEC")
        if type(ra)==float: 
            ral=scipy.zeros(1)
            ral[0]=ra
            ra=ral
        if type(dec)==float: 
            decl=scipy.zeros(1)
            decl[0]=dec
            dec=decl

        c62=math.cos(62.6*D2R)
        s62=math.sin(62.6*D2R)
        
        b=scipy.sin(dec*D2R)*c62
        b-=scipy.cos(dec*D2R)*scipy.sin((ra-282.25)*D2R)*s62
        b=scipy.arcsin(b)*R2D
        
        cosb=scipy.cos(b*D2R)
        #l minus 33 degrees
        lm33=(scipy.cos(dec*D2R)/cosb)*scipy.cos((ra-282.25))
        l=scipy.arccos(lm33)*R2D+33.0
        return b,l

def rect_to_cyl(X,Y,Z):
    """
    NAME:

       rect_to_cyl

    PURPOSE:

       convert from rectangular to cylindrical coordinates

    INPUT:

       X, Y, Z - rectangular coordinates

    OUTPUT:

       [:,3] R,phi,z

    HISTORY:

       2010-09-24 - Written - Bovy (NYU)

    """
    R= sc.sqrt(X**2.+Y**2.)
    phi= sc.arcsin(Y/R)
    if isinstance(X,float) and X < 0.:
        phi= m.pi-phi
    elif isinstance(X,sc.ndarray):
        phi[(X < 0.)]= m.pi-phi[(X < 0.)]
    return (R,phi,Z)

def pix2sky(header,x,y):
	hdr_info = parse_header(header)
	x0 = x-hdr_info[1][0]+1.	# Plus 1 python->image
	y0 = y-hdr_info[1][1]+1.
	x0 = x0.astype(scipy.float64)
	y0 = y0.astype(scipy.float64)
	x = hdr_info[2][0,0]*x0 + hdr_info[2][0,1]*y0
	y = hdr_info[2][1,0]*x0 + hdr_info[2][1,1]*y0
	if hdr_info[3]=="DEC":
		a = x.copy()
		x = y.copy()
		y = a.copy()
		ra0 = hdr_info[0][1]
		dec0 = hdr_info[0][0]/raddeg
	else:
		ra0 = hdr_info[0][0]
		dec0 = hdr_info[0][1]/raddeg
	if hdr_info[5]=="TAN":
		r_theta = scipy.sqrt(x*x+y*y)/raddeg
		theta = arctan(1./r_theta)
		phi = arctan2(x,-1.*y)
	elif hdr_info[5]=="SIN":
		r_theta = scipy.sqrt(x*x+y*y)/raddeg
		theta = arccos(r_theta)
		phi = artan2(x,-1.*y)
	ra = ra0 + raddeg*arctan2(-1.*cos(theta)*sin(phi-pi),
				   sin(theta)*cos(dec0)-cos(theta)*sin(dec0)*cos(phi-pi))
	dec = raddeg*arcsin(sin(theta)*sin(dec0)+cos(theta)*cos(dec0)*cos(phi-pi))

	return ra,dec

def elaz2radec_lst(el, az, lst, lat = 38.43312) :
    """DO NOT USE THIS ROUTINE FOR ANTHING THAT NEEDS TO BE RIGHT.  IT DOES NOT
    CORRECT FOR PRECESSION.

    Calculates the Ra and Dec from elavation, aximuth, LST and Latitude.

    This function is vectorized with numpy so should be fast.  Standart numpy
    broadcasting should also work.

    All angles in degrees, lst in seconds. Latitude defaults to GBT.
    """

    # Convert everything to radians.
    el = sp.radians(el)
    az = sp.radians(az)
    lst = sp.array(lst, dtype = float)*2*sp.pi/86400
    lat = sp.radians(lat)
    # Calculate dec.
    dec = sp.arcsin(sp.sin(el)*sp.sin(lat) +
                    sp.cos(el)*sp.cos(lat)*sp.cos(az))
    # Calculate the hour angle
    ha = sp.arccos((sp.sin(el) - sp.sin(lat)*sp.sin(dec)) /
                   (sp.cos(lat)*sp.cos(dec)))
    ra = sp.degrees(lst - ha) % 360

    return ra, sp.degrees(dec)

def xyz2longlat(x,y,z):
    """ converts cartesian x,y,z coordinates into spherical longitude and latitude """
    r = sc.sqrt(x*x + y*y + z*z)
    long = sc.arctan2(y,x)
    d = sc.sqrt(x*x + y*y)
    lat = sc.arcsin(z/r)
    return long*deg, lat*deg, r

def FhklDWBA4(x,y,z,h,k,l=None,occ=None,alphai=0.2,alphaf=None,substrate=None,wavelength=1.0,e_par=0.,e_perp=1.0,gpu_name="CPU",use_fractionnal=True,language="OpenCL",cl_platform="",separate_paths=False):
  """
  Calculate the grazing-incidence X-ray scattered intensity taking into account
  4 scattering paths, for a nanostructure object located above a given substrate.
  The 5th path is the scattering from the substrate, assumed to be below the
  interface at z=0.
  
  x,y,z: coordinates of the atoms in fractionnal coordinates (relative to the 
         substrate unit cell)- if use_fractionnal==False, these should be given in Angstroems
  h,k,l: reciprocal space coordinates. If use_fractionnal==False, these should be given
         in inverse Angstroems (multiplied by 2pi - physicist 'k' convention, |k|=4pisin(theta)/lambda,
         i.e. these correspond to k_x,k_y,k_z).
  alphai, alphaf: incident and outgoing angles, in radians
  substrate: the substrate material, as a pynx.gid.Crystal object - this will be used
             to calculate the material refraction index.
  wavelength: in Angstroems
  e_par,e_perp: percentage of polarisation parallel and perpendicular to the incident plane
  use_fractionnal: if True (the default), then coordinates for atoms and reciprocal
                   space are given relatively to the unit cell, otherwise in Angstroems
                   and 2pi*inverse Angstroems.
  
  Note: Either l *OR* alphaf must be supplied - it is assumed that the lattice
  coordinates are such that the [001] direction is perpendicular to the surface.
  """
  nrj=W2E(wavelength)
  if use_fractionnal: 
    c=substrate.uc.parameters()[2]
    s_fact=1.0
  else: 
    c=2*pi
    s_fact=1/c
  if alphaf==None:
    # alphaf, computed from l: l.c* = (sin(alpha_f) + sin(alpha_i))/wavelength
    alphaf=scipy.arcsin(l/c*wavelength-scipy.sin(alphai))

  # Incident wave
  w=Wave(alphai,e_par,e_perp,nrj)
  dw=DistortedWave(None,substrate,w)
  # Reflected wave after the dot
  w1=Wave(alphaf,e_par,e_perp,nrj)
  dw1=DistortedWave(None,substrate,w1)

  # First path, direct diffraction
  l=c*scipy.sin(alphaf+alphai)/wavelength
  f1=gpu.Fhkl_thread(h*s_fact,k*s_fact,l*s_fact,x,y,z,occ=occ,gpu_name=gpu_name,language=language,cl_platform=cl_platform)[0]

  # Second path, reflection before
  l=c*scipy.sin(alphaf-alphai)/wavelength
  f2=gpu.Fhkl_thread(h*s_fact,k*s_fact,l*s_fact,x,y,z,occ=occ,gpu_name=gpu_name,language=language,cl_platform=cl_platform)[0]*dw.Riy

  # Third path, reflection after dot
  l=c*scipy.sin(-alphaf+alphai)/wavelength
  f3=gpu.Fhkl_thread(h*s_fact,k*s_fact,l*s_fact,x,y,z,occ=occ,gpu_name=gpu_name,language=language,cl_platform=cl_platform)[0]*dw1.Riy

  # Fourth path, reflection before and after dot
  l=c*scipy.sin(-alphaf-alphai)/wavelength
  f4=gpu.Fhkl_thread(h*s_fact,k*s_fact,l*s_fact,x,y,z,occ=occ,gpu_name=gpu_name,language=language,cl_platform=cl_platform)[0]*dw.Riy*dw1.Riy
  if separate_paths: return f1,f2,f3,f4
  return f1+f2+f3+f4

def sin_seq(A = 1.0, f = 1.0, T = 1.0, delta = 0.1):
    t = 0.0
    events = []
    # Note use of round() to compensate for floating-point arithmetic errors
    # that lead to inexact results.
    while t <= T:
        v = delta * floor(round(A*sin(2*pi*f*t) / delta,10))
        events += [(t, v)]
        tm = fmod(abs(t), 0.5/f)
        if tm < 0.25/f:
            vq = delta * (floor(round(A*sin(2*pi*f*tm) / delta,10)) + 1.0)
            dt = arcsin(vq/A)/(2*pi*f) - tm
        else:
            vq = delta * (floor(round(A*sin(2*pi*f*tm) / delta,10)) - 1.0)
            dt = (pi - arcsin(vq/A))/(2*pi*f) - tm             
        t += dt
    return array(events)

 def azimuth(self, date):
     dec = self.declination(date)
     ha = self.hour_angle(date)
     az = sp.arcsin(sp.cos(dec) * sp.sin(ha) / sp.cos(self.elevation(date)))
     if (sp.cos(ha) >= (sp.tan(dec) / sp.tan(self.lat))):
         return az
     else:
         return (sp.pi - az)

 def _arcsin_sqrt_transform(self, verbose=False):
     a = sp.array(self.values)
     if min(a) < 0 or max(a) > 1:
         log.debug('Some values are outside of range [0,1], hence skipping transformation!')
         return False
     else:
         vals = sp.arcsin(sp.sqrt(a))
     self._perform_transform(vals,"arcsin")
     return True

def sky2pix(header,ra,dec):
	hdr_info = parse_header(header)
	if scipy.isscalar(ra):
		ra /= raddeg
		dec /= raddeg
	else:
		ra = ra.astype(scipy.float64)/raddeg
		dec = dec.astype(scipy.float64)/raddeg
	if hdr_info[3]=="DEC":
		ra0 = hdr_info[0][1]/raddeg
		dec0 = hdr_info[0][0]/raddeg
	else:
		ra0 = hdr_info[0][0]/raddeg
		dec0 = hdr_info[0][1]/raddeg

	phi = pi + arctan2(-1*cos(dec)*sin(ra-ra0),sin(dec)*cos(dec0)-cos(dec)*sin(dec0)*cos(ra-ra0))
	argtheta = sin(dec)*sin(dec0)+cos(dec)*cos(dec0)*cos(ra-ra0)
	if scipy.isscalar(argtheta):
		theta = arcsin(argtheta)
	else:
		argtheta[argtheta>1.] = 1.
		theta = arcsin(argtheta)

	if hdr_info[5]=="TAN":
		r_theta = raddeg/tan(theta)
		x = r_theta*sin(phi)
		y = -1.*r_theta*cos(phi)
	elif hdr_info[5]=="SIN":
		r_theta = raddeg*cos(theta)
		x = r_theta*sin(phi)
		y = -1.*r_theta*cos(phi)
	if hdr_info[3]=="DEC":
		a = x.copy()
		x = y.copy()
		y = a.copy()
	inv = linalg.inv(hdr_info[2])
	x0 = inv[0,0]*x + inv[0,1]*y
	y0 = inv[1,0]*x + inv[1,1]*y

	x = x0+hdr_info[1][0]-1
	y = y0+hdr_info[1][1]-1

	return x,y

def snell(n_1, n_2, th_1):
    """
    Return angle theta in layer 2 with refractive index n_2, assuming
    it has angle th_1 in layer with refractive index n_1. Use Snell's law. Note
    that "angles" may be complex!!
    """
    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin!! (They
    # give different results e.g. for arcsin(2).)
    # Use real_if_close because e.g. arcsin(2 + 1e-17j) is very different from
    # arcsin(2) due to branch cut
    return sp.arcsin(np.real_if_close(n_1 * np.sin(th_1) / n_2))

def list_snell(n_list,th_0):
    """
    return list of angle theta in each layer based on angle th_0 in layer 0,
    using Snell's law. n_list is index of refraction of each layer. Note that
    "angles" may be complex!!
    """
    #Important that the arcsin here is scipy.arcsin, not numpy.arcsin!! (They
    #give different results e.g. for arcsin(2).)
    #Use real_if_close because e.g. arcsin(2 + 1e-17j) is very different from
    #arcsin(2) due to branch cut
    return sp.arcsin(np.real_if_close((n_list[:,0]*np.sin(th_0))[:,None] / n_list))

def cart2sphere(coordlist):
    r2d = 180.0/sp.pi
    d2r = sp.pi/180.0
    
    X_vec = coordlist[:,0]
    Y_vec = coordlist[:,1]
    Z_vec = coordlist[:,2]
    R_vec = sp.sqrt(X_vec**2+Y_vec**2+Z_vec**2)
    Az_vec = sp.arctan2(X_vec,Y_vec)*r2d
    El_vec = sp.arcsin(Z_vec/R_vec)*r2d
    sp_coords = sp.array([R_vec,Az_vec,El_vec]).transpose()
    return sp_coords

def lbToEq (l_deg, b_deg):   
    """ Converts galactic l,b in to Equatorial ra, dec; from Binney and Merrifield, p. 31;  
    l, b must be arrays of same shape"""
    l, b = (l_deg*rad), (b_deg*rad)
    # Conversion Code
    t = lCP - l
    dec = sc.arcsin(sc.sin(decGP)*sc.sin(b) + sc.cos(decGP)*sc.cos(b)*sc.cos(t) )
    r = sc.arctan2( (sc.cos(b)*sc.sin(t)),
                    ( (sc.cos(decGP)*sc.sin(b)) - (sc.sin(decGP)*sc.cos(b)*sc.cos(t)))  )
    if type(r) != type(arr):  r = sc.array([r])
    for i in range(len(r)):
        r[i] = angle_bounds((r[i] + raGP)*deg)
    return r, (dec*deg)

def EqTolb (ra_deg, dec_deg):
    """ Converts equatorial ra, dec, into galactic l,b;  from Binney and Merrifield, p. 31
    following the method of http://star-www.st-and.ac.uk/~spd3/Teaching/AS3013/programs/radec2lb.f90
    NOT QUITE - I use arctan2 method instead"""
    ra, dec = (ra_deg*rad), (dec_deg*rad)
    # Conversion Code
    r = (ra - raGP)
    b = sc.arcsin( sc.sin(decGP)*sc.sin(dec) + sc.cos(decGP)*sc.cos(dec)*sc.cos(r) )
    t = sc.arctan2((sc.cos(dec)*sc.sin(r)),
                   (sc.cos(decGP)*sc.sin(dec) - sc.sin(decGP)*sc.cos(dec)*sc.cos(r)) )
    l = (lCP - t)
    b, l = angle_bounds2((b*deg), (l*deg))
    return l, b