def lla2ecef(lla: Sequence[float], cst: ConstantsFile, lla_as_degrees: bool=False) -> Tuple[float, float, float]:
    """
    converts LLA (Latitude, Longitude, Altitude) coordinates
    to ECEF (Earth-Centre, Earth-First) XYZ coordinates.
    """

    lat, lon, alt = lla
    if lla_as_degrees:
        lat = radians(lat)
        lon = radians(lon)

    a = cst.semi_major_axis
    b = cst.semi_minor_axis
    # calc. ellipsoid flatness
    f = (a - b) / a
    # calc. eccentricity
    e = sqrt(f * (2 - f))

    # Calculate length of the normal to the ellipsoid
    N = a / sqrt(1 - (e * sin(lat)) ** 2)
    # Calculate ecef coordinates
    x = (N + alt) * cos(lat) * cos(lon)
    y = (N + alt) * cos(lat) * sin(lon)
    z = (N * (1 - e ** 2) + alt) * sin(lat)
    # Return the ecef coordinates
    return x, y, z

def _genBFEdgeZero(plasma, zeros, rcent, zcent):
    """ this will absolutely need to be rewritten"""

    theta = scipy.linspace(-scipy.pi,scipy.pi,zeros)
    cent = geometry.Point(geometry.Vecr([rcent,0,zcent]),plasma)
    zerobeam = []
    outline = []
    for i in xrange(len(plasma.norm.s)-1):
        outline += [geometry.Vecx([plasma.sagi.s[i],
                                   0,
                                   plasma.norm.s[i]])-cent]
        
    for i in xrange(zeros):
        temp2 = geometry.Vecr([scipy.cos(theta[i]),
                               0,
                               scipy.sin(theta[i])])
        s = 0
        for j in outline:
            temp4 = j*temp2
            if temp4 > s:
                s = temp4

        temp2.s = s
        zerobeam += [Ray(geometry.Point(cent+temp2,
                                        plasma),
                         geometry.Vecr([scipy.sin(theta[i]),
                                        0,
                                        -scipy.cos(theta[i])]))]

    return zerobeam

def cyl_to_rect_vec(vr,vt,vz,phi):
    """
    NAME:

       cyl_to_rect_vec

    PURPOSE:

       transform vectors from cylindrical to rectangular coordinate vectors

    INPUT:

       vr - radial velocity

       vt - tangential velocity

       vz - vertical velocity

       phi - azimuth

    OUTPUT:

       vx,vy,vz

    HISTORY:

       2011-02-24 - Written - Bovy (NYU)

    """
    vx= vr*sc.cos(phi)-vt*sc.sin(phi)
    vy= vr*sc.sin(phi)+vt*sc.cos(phi)
    return (vx,vy,vz)

 def test_correlate(self) :
     Data = self.blocks[0]
     Data.calc_freq()
     map = self.map
     gain = 3.45
     const = 2.14
     # Set all data = gain*(cos(time_ind)).
     Data.data[:,:,:,:] = gain*sp.cos(sp.arange(1,11)
                                 [:,sp.newaxis,sp.newaxis,sp.newaxis])
     # Explicitly set time mean to something known.
     Data.data -= ma.mean(Data.data, 0)
     Data.data += gain*const*Data.freq/800.0e6
     # Now the Map.
     map[:,:,:] = 0.0
     # Set 10 pixels to match cos part of data.
     map[:, range(10), range(10)] = (
                 sp.cos(sp.arange(1,11)[None, :]))
     map[:, range(10), range(10)] -= ma.mean(
         map[:, range(10), range(10)], 1)[:, None]
     # Give Map a mean to test things out. Should really have no effect.
     map[...] += 0.352*map.get_axis('freq')[:, None, None]/800.0e6
     # Rig the pointing to point to those 10 pixels.
     def rigged_pointing() :
         Data.ra = map.get_axis('ra')[range(10)]
         Data.dec = map.get_axis('dec')[range(10)]
     Data.calc_pointing = rigged_pointing
     solved_gains = smd.sub_map(Data, map, correlate=True)
     # Now data should be just be gain*const*f, within machine precision.
     Data.data /= gain*Data.freq/800.0e6
     self.assertTrue(sp.allclose(Data.data[:,:,:,:], const))
     self.assertTrue(sp.allclose(solved_gains, gain))

 def f(self, x):
     res = 0
     for i in range(self.xdim):
         Ai = sum(self.A[i] * sin(self.alphas) + self.B[i] * cos(self.alphas))
         Bix = sum(self.A[i] * sin(x) + self.B[i] * cos(x))
         res += (Ai - Bix) ** 2
     return res

def CalcXY2GPSParam_2p(x1,x2,g1,g2,K=[0,0]):
	# Kx = dLng/dx; Ky = dlat/dy;
	# In China:
	# Kx = (133.4-1.2*lat)*1e3
	# Ky = (110.2+0.002*lat)*1e3

	X1 = array(x1)
	Y1 = array(g1)
	X2 = array(x2)
	Y2 = array(g2)
	detX = X2-X1
	detY = Y2-Y1
	lat = Y1[1]
	if K[0] == 0:
		Kx = (133.4-1.2*lat)*1e3
		Ky = (110.2+0.002*lat)*1e3
		K = array([Kx,Ky])
	else:
		Kx = K[0]
		Ky = K[1]
	detKY = detY*K

	alpha =  myArctan(detX[0],detX[1]) - myArctan(detKY[0],detKY[1])
	A = array([[sp.cos(alpha),sp.sin(alpha)],[-sp.sin(alpha),sp.cos(alpha)]])
	X01 = X1 - dot(linalg.inv(A),Y1*K) 
	X02 = X2 - dot(linalg.inv(A),Y2*K)
	X0 = (X01+X02) /2

	return A,X0,K

def sparse_orth(d):
    """ Constructs a sparse orthogonal matrix.
    
    The method is described in:
    Gi-Sang Cheon et al., Constructions for the sparsest orthogonal matrices,
    Bull. Korean Math. Soc 36 (1999) No.1 pp.199-129
    """
    from scipy.sparse import eye
    from scipy import r_, pi, sin, cos

    if d % 2 == 0:
        seq = r_[0:d:2, 1:d - 1:2]
    else:
        seq = r_[0:d - 1:2, 1:d:2]
    Q = eye(d, d).tocsc()
    for i in seq:
        theta = random() * 2 * pi
        flip = (random() - 0.5) > 0;
        Qi = eye(d, d).tocsc()
        Qi[i, i] = cos(theta)
        Qi[(i + 1), i] = sin(theta)
        if flip > 0:
            Qi[i, (i + 1)] = -sin(theta)
            Qi[(i + 1), (i + 1)] = cos(theta)
        else:
            Qi[i, (i + 1)] = sin(theta)
            Qi[(i + 1), (i + 1)] = -cos(theta)
        Q = Q * Qi;
    return Q

def evaluateSphericalVariation (phi,theta,cntPhi,cntTheta):
  global conf,cons
  success = False

  alpha0     =sp.zeros([dim['alpha']],complex)
  alpha0[conf['id0']]=sp.cos(phi)*sp.sin(theta)+0j
  alpha0[conf['id1']]=sp.sin(phi)*sp.sin(theta)+0j
  alpha0[conf['id2']]=sp.cos(theta)+0j

#  if (sp.absolute(alpha0[conf['id0']]) <= 1e-10): 
#    alpha0[conf['id0']]=0.0+0j
#  if (sp.absolute(alpha0[conf['id1']]) <= 1e-10): 
#    alpha0[conf['id1']]=0.0+0j
#  if (sp.absolute(alpha0[conf['id2']]) <= 1e-10): 
#    alpha0[conf['id2']]=0.0+0j
#  
  # normalize coefficients for alpha -> defines net-power
  alpha0[:]=alpha0[:]/sp.linalg.norm(alpha0[:])*cons['alpha_norm']
  __,res   = MemoryPulseFunctional.evaluateFunctional(alpha0,1.0+0j)
 
  myRes    = sp.zeros([conf['entries']+3])
  myRes[0] = alpha0[conf['id0']].real
  myRes[1] = alpha0[conf['id1']].real
  myRes[2] = alpha0[conf['id2']].real
  myRes[3:]= res
 
  print "### spherical map: phi/pi={0:5.3f}, theta/pi={1:5.3f}, fun={2:f}".format(phi/sp.pi,theta/sp.pi,myRes[conf['funval']])

  myRes[conf['funval']] = min(conf['cutoff'],res[conf['funval']])

  return myRes,cntPhi,cntTheta

 def __init__(self,id,matName,orientation,source=0.0):
     self.id = id
     self.matName = matName
     self.orientation = orientation
     self.source = source
     self.T = array([[cos(orientation),-sin(orientation)],
                     [sin(orientation), cos(orientation)]])

    def form_point_set(self, histo, point_set):
        (slices, numbins) = histo.shape
        phases = numpy.arange(numbins)
        phases = phases * (360. / numbins)
        phases += phases[1] / 2.
        phi_step = phases[0]
        
        for time in xrange(slices):
            z = float(time)
            for bin in xrange(numbins):
                r = histo[time,bin]
                theta = phi_step * (bin+1)
                theta *= (scipy.pi / 180.)
                x = r*scipy.cos(theta)
                y = r*scipy.sin(theta)
                point_set.InsertNextPoint(x, y, z)

            for bin in xrange(numbins):
                curbin = bin
                lastbin = bin-1
                if lastbin < 0:
                    lastbin = numbins-1

                r = (histo[time,bin] -  histo[time,lastbin]) / 2.
                theta = curbin * 360. / numbins
                x = r*scipy.cos(theta)
                y = r*scipy.sin(theta)
                point_set.InsertNextPoint(x, y, z)

def binary_ephem(P, T, e, a, i, O_node, o_peri, t):
   # Grados a radianes
   d2rad = pi/180.
   rad2d = 180./pi
   i = i*d2rad
   O_node = (O_node*d2rad)%(2*pi)
   o_peri = (o_peri*d2rad)%(2*pi)
 
   # Anomalia media
   M = ((2.0*pi)/P)*(t - T)  # radianes
    
   if M >2*pi: M = M - 2*pi 
   M=M%(2*pi)

   # Anomalia excentrica (1ra aproximacion)
   E0 = M  + e*sin(M) + (e**2/M) * sin(2.0*M)

   for itera in range(15):
      M0 = E0 - e*sin(E0)
      E0 = E0 + (M-M0)/(1-e*cos(E0))

   true_anom = 2.0*arctan(sqrt((1+e)/(1-e))*tan(E0/2.0))
 
   #radius = (a*(1-e**2))/(1+e*cos(true_anom))
   radius = a*(1-e*cos(E0))

   theta = arctan( tan(true_anom + o_peri)*cos(i) ) + O_node
   rho = radius * (cos(true_anom + o_peri)/cos(theta - O_node))
   
   # revuelve rho ("), theta (grad), Anomalia excentrica (grad), Anomalia verdadera (grad)
   return rho, (theta*rad2d)%360. #, E0*rad2d, M*rad2d, true_anom*rad2d

 def rotate(self, angle, mask=None):
     """Rotate the grids (arena centered)
     
     Grids to be rotated can be optionally specified by bool/index array
     *mask*, otherwise population is rotated. Specified *angle* can be a
     scalar value to be applied to the population or a population- or
     mask-sized array depending on whether *mask* is specified.
     """
     rot2D = lambda psi: [[cos(psi), sin(psi)], [-sin(psi),  cos(psi)]]
     if mask is not None and type(mask) is np.ndarray:
         if mask.dtype.kind == 'b':
             mask = mask.nonzero()[0]
         if type(angle) is np.ndarray and angle.size == mask.size:
             for i,ix in enumerate(mask):
                 self._phi[ix] = np.dot(self._phi[ix], rot2D(angle[i]))
         elif type(angle) in (int, float, np.float64):
             angle = float(angle)
             self._phi[mask] = np.dot(self._phi[mask], rot2D(angle))
         else:
             raise TypeError, 'angle must be mask-sized array or float'
         self._psi[mask] = np.fmod(self._psi[mask]+angle, 2*pi)
     elif mask is None:
         if type(angle) is np.ndarray and angle.size == self.num_maps:
             for i in xrange(self.num_maps):
                 self._phi[i] = np.dot(self._phi[i], rot2D(angle[i]))
         elif type(angle) in (int, float, np.float64):
             angle = float(angle)
             self._phi = np.dot(self._phi, rot2D(angle))
         else:
             raise TypeError, 'angle must be num_maps array or float'
         self._psi = np.fmod(self._psi+angle, 2*pi)
     else:
         raise TypeError, 'mask must be bool/index array'

 def __init__(self,alphai,eparall,eperp,nrj):
   """
       Incident wave above a surface.
       Coordinates:
         - z is perpendicular to the surface, >0 going UP (different from H Dosch's convention)
         - x is the projection of the wavevector on the surface
         - y is parallel to the surface
       
       alphai: incident angle, with respect to the surface
       eparallel: component of the electric field parallel to the incident plane (vertical plane)
       eperp: component of the electric field perpendicular to the incident plane (along y)
       nrj: values of the energy of the incident wave, in eV
       
       alphai *or* nrj can be arrays, but not together
   """
   self.alphai=alphai
   self.eparall=eparall
   self.eperp=eperp
   self.ex=scipy.sin(alphai)*eparall
   self.ey=eperp
   self.ez=scipy.cos(alphai)*eparall
   self.kx= 2*pi/W2E(nrj)*scipy.cos(alphai)
   self.ky= 2*pi/W2E(nrj)*0
   self.kz=-2*pi/W2E(nrj)*scipy.sin(alphai)
   self.nrj=nrj

def ned2ecef(lat, lon, alt, n, e, d):
    X0, Y0, Z0 = coord.geodetic2ecef(lat, lon, alt)
    lat, lon = radians(lat), radians(lon)
    
    pitch = math.pi/2 + lat
    yaw = -lon 
    
    my = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (cos(pitch), 0, -sin(pitch),
         0,1,0,
         sin(pitch), 0, cos(pitch)))
    
    mz = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (cos(yaw), sin(yaw),0,
         -sin(yaw),cos(yaw),0,
         0,0,1))
    
    mr = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (-cos(lon)*sin(lat), -sin(lon), -cos(lat) * cos(lon), 
         -sin(lat)*sin(lon), cos(lon), -sin(lon)*cos(lat),
         cos(lat), 0, -sin(lat)))
    
    geo = mat('[%f; %f; %f]' % (X0, Y0, Z0))
    ned = mat('[%f; %f; %f]' % (n, e, d))
    res = mr*ned + geo
    return res[0], res[1], res[2]  

 def setUp(self):
     # Make a positive definite noise matrix, clean map, and dirty_map.
     self.nra = 10
     self.ndec = 5
     self.nf = 20
     self.shape = (self.nf, self.nra, self.ndec)
     self.size = self.nra * self.ndec * self.nf
     # Clean map.
     clean_map = sp.empty(self.shape, dtype=float)
     clean_map = al.make_vect(clean_map, axis_names=('freq', 'ra', 'dec'))
     clean_map[...] = sp.sin(sp.arange(self.nf))[:,None,None]
     clean_map *= sp.cos(sp.arange(self.nra))[:,None]
     clean_map *= sp.cos(sp.arange(self.ndec))
     # Noise inverse matrix.
     noise_inv = sp.empty(self.shape * 2, dtype=float)
     noise_inv = al.make_mat(noise_inv, axis_names=('freq', 'ra', 'dec')*2,
                             row_axes=(0, 1, 2), col_axes=(3, 4, 5))
     rand_mat = rand.randn(*((self.size,) * 2))
     information_factor = 1.e6  # K**-2
     rand_mat = sp.dot(rand_mat, rand_mat.transpose()) * information_factor
     noise_inv.flat[...] = rand_mat.flat
     # Dirty map.
     dirty_map = al.partial_dot(noise_inv, clean_map)
     # Store in self.
     self.clean_map = clean_map
     self.noise_inv = noise_inv
     self.dirty_map = dirty_map

def test_time_integration():

    dt = 0.01
    tmax = 1

    def check_time_integration(res, y0, exact):
        ts, ys = time_integrate(res, y0, dt, tmax)

        exacts = [exact(t) for t in ts]

        print(ts, ys, exacts)

        utils.assert_list_almost_equal(ys, exacts, 1e-3)

    tests = [
        (lambda t, y, dy: y - dy, 1.0, lambda t: array(exp(t))),
        (lambda t, y, dy: y + dy, 1.0, lambda t: array(exp(-t))),
        (
            lambda t, y, dy: array([-0.1 * sin(t), y[1]]) - dy,
            array([0.1 * cos(0.0), exp(0.0)]),
            lambda t: array([0.1 * cos(t), exp(t)]),
        ),
    ]

    for r, y0, exact in tests:
        yield check_time_integration, r, y0, exact

def lat_lon_2_vertex(lat,lon):
	"""
	A routine to return the location of a detector's vertex in 3D
	Cartesean Coordinates given the latitude and longitude of the
	detector. This routine approximates teh Earth as a sphere.
	"""
	return (metric.R_earth*cos(lon)*cos(lat), metric.R_earth*sin(lon)*cos(lat), metric.R_earth*sin(lat))

def radec_to_lb_single(ra,dec,T,degree=False):
    """
    NAME:
       radec_to_lb_single
    PURPOSE:
       transform from equatorial coordinates to Galactic coordinates for a single pair of ra,dec
    INPUT:
       ra - right ascension
       dec - declination
       T - epoch dependent transformation matrix (dictionary)
       degree - (Bool) if True, ra and dec are given in degree and l and b will be as well
    OUTPUT:
       l,b
    HISTORY:
       2009-11-12 - Written - Bovy (NYU)
    """
    T=T['T']
    if degree:
        thisra= ra/180.*sc.pi
        thisdec= dec/180.*sc.pi
    else:
        thisra= ra
        thisdec= dec
    XYZ=sc.array([sc.cos(thisdec)*sc.cos(thisra),sc.cos(thisdec)*sc.sin(thisra),sc.sin(thisdec)])
    galXYZ= sc.dot(T,XYZ)
    b= m.asin(galXYZ[2])
    l= m.atan(galXYZ[1]/galXYZ[0])
    if galXYZ[0]/sc.cos(b) < 0.:
        l+= sc.pi
    if l < 0.:
        l+= 2.*sc.pi
    if degree:
        return (l/sc.pi*180.,b/sc.pi*180.)
    else:
        return (l,b)

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 xyzfield(gcoefs, hcoefs, phi, theta, rparam=1.0, order=13):
    # no usar esta función, es peor en rendimiento
    x, y, z = 0, 0, 0

    legendre, dlegendre = scipy.special.lpmn(order + 1, order + 1, scipy.cos(theta))

    for l in range(1, order + 1):
        for m in range(0, l + 1):
            deltax = (
                rparam ** (l + 2)
                * (gcoefs[m, l] * scipy.cos(m * phi) + hcoefs[m, l] * scipy.sin(m * phi))
                * dlegendre[m, l]
                * (-scipy.sin(theta))
            )
            deltay = (
                rparam ** (l + 2)
                * (gcoefs[m, l] * scipy.sin(m * phi) - hcoefs[m, l] * scipy.cos(m * phi))
                * m
                * legendre[m, l]
                / (scipy.sin(theta))
            )
            deltaz = (
                rparam ** (l + 2)
                * (l + 1)
                * (gcoefs[m, l] * scipy.cos(m * phi) + hcoefs[m, l] * scipy.sin(m * phi))
                * legendre[m, l]
            )

            x += deltax
            y += deltay
            z += deltaz

    return (x, y, z)