def ellipse(emittance, beta, alpha, gamma):
     phi = linspace(0, 2*const.pi, 1e3)
     a = sqrt(emittance/2*(beta+gamma+sqrt((beta+gamma)**2-4)))
     b = sqrt(emittance/2*(beta+gamma-sqrt((beta+gamma)**2-4)))
     if alpha > 0:
         PHI = acos(+sqrt((beta-b/a)*emittance/(a**2-b**2)))
     else:
         PHI = acos(-sqrt((beta-b/a)*emittance/(a**2-b**2)))
     pos = a*cos(phi)*cos(PHI)+b*sin(phi)*sin(PHI)
     mom = -a*cos(phi)*sin(PHI)+b*sin(phi)*cos(PHI)
     return pos, mom

def getlatlon(x,y,z):
    
    #radEarth = 6378145. #[m]
    rmag = np.sqrt(x**2+y**2+z**2)
    longitude = np.arctan2(y,x)*180./np.pi
    latitude = 90. - np.acos(z/rmag)*180./np.pi
    return (latitude,longitude)

 def point2unit_spherical(point):
     ub = dot(self.axinv,(point-self.origin)) 
     u=zeros((self.DIM,))
     u[0] = sqrt(ub[0]*ub[0]+ub[1]*ub[1]+ub[2]*ub[2]) 
     u[1] = atan2(ub[1],ub[0])*o2pi+.5 
     u[2] = acos(ub[2]/u[0])*o2pi*2.0 
     return u

  def develop(self):
    """
    Develop the 3D shape into 2D
    """
    return
    ## Every flat generatrice is 2 2D points
    flatgen = np.zeros((n_gen, 2, 2))
    flatgen[0,0] = np.zeros(2)
    flatgen[0,1] = np.array([0, length(self.gen[0])])
    ## 1st generatrice is already in same plane as 0st
    u1 = self.support[1] - self.support[0]
    c = cos(self.gen[1], u1)
    s = cos(self.gen[1], u1)
    d1 = np.array(length(u1)*s + length(u1)*c)
    flatgen[1,0] = flatgen[0,0] + d1
    
    for i in np.range(1, n_gen):
        ## v is rotation axis
        v = self.gen[i]
        u1 = self.support[i] - self.support[i-1]
        u2 = self.support[i] + self.gen[i] - self.support[i-1] - self.gen[i-1] 
	## Pick the best triangle for reference plane
    	c1 = np.cross(u1, v)
        c2 = np.cross(v, u2)
        if length(c1) > length(c2):
           u = u1
        else:
           u = u2
        ## Compute angle between gen and v
        c = cos(self.gen[i+1], v)
        alpha = np.acos(c)
        ## alpha is preserved in rotation, such as length of gen
        
    return self

    def drawKochPyramid(self, level):
        angle = rad2deg(acos(1/3))

        if (level == 0):
            glBegin(GL_TRIANGLES)
            self.basicTriangle()
            glEnd()
        else:
            glPushMatrix()
            glScalef(1.0/3, 1.0/3, 1.0/3)

            self.drawKochPyramid(level - 1)

            glTranslatef(1.0, 0, 0)
            self.drawKochPyramid(level - 1)

            glTranslatef(1.0, 0, 0)
            self.drawKochPyramid(level - 1)

            glTranslatef(-0.5, sqrt(3)/2.0, 0)
            self.drawKochPyramid(level - 1)

            self.innerPyra(level, angle)

            glTranslatef(-1,0, 0)
            self.drawKochPyramid(level - 1)

            self.innerPyra(level, angle)

            glTranslatef(0.5, sqrt(3)/2.0, 0)
            self.drawKochPyramid(level - 1)

            self.innerPyra(level, angle)

            glPopMatrix()

def chessboard_binary_noise(hologram):
    """Dump power into a binary chessboard grating"""
    signs = np.zeros(hologram.shape, dtype=np.float)
    signs[0::2, 1::2] = 1
    signs[1::2, 0::2] = 1
    # signs is a chessboard grating
    return nearest(hologram) + signs * np.acos(normalised_amplitude(hologram))

def get_phi_theta(vec):
    """ returns a tupel of the phi and theta angles of the given vector
    """
    try:
        return (atan2(vec[1], vec[0]), acos(np.clip(vec[2] / length(vec), -1, 1))) * u.rad
    except ValueError:
        return (0, 0)

def angles_from_matrix(rot_matrix):
    if rot_matrix.shape == (2, 2):
        theta = np.atan2(rot_matrix[1, 0], rot_matrix[0, 0])
        return theta,
    elif rot_matrix.shape == (3, 3):
        if rot_matrix[2, 2] == 1.:  # cannot use last row and column
            theta = 0.
            # upper-left block is 2d rotation for phi + psi, so one needs
            # to be fixed
            psi = 0.
            phi = np.atan2(rot_matrix[1, 0], rot_matrix[0, 0])
            if phi < 0:
                phi += 2 * np.pi  # in [0, 2pi)
        else:
            phi = np.atan2(rot_matrix[0, 2], -rot_matrix[1, 2])
            psi = np.atan2(rot_matrix[2, 0], rot_matrix[2, 1])
            theta = np.acos(rot_matrix[2, 2])

            if phi < 0. or psi < 0.:
                phi += np.pi
                psi += np.pi
                theta = -theta

        return phi, theta, psi

    else:
        raise ValueError('shape of `rot_matrix` must be (2, 2) or (3, 3), '
                         'got {}'.format(rot_matrix.shape))

    def beta(self,m1,d,g=1.4,i=0):
        p=-(m1*m1+2.)/m1/m1-g*np.sin(d)*np.sin(d)
        q=(2.*m1*m1+1.)/ np.pow(m1,4.)+((g+1.)*(g+1.)/4.+
                                          (g-1.)/m1/m1)*np.sin(d)*np.sin(d)
        r=-np.cos(d)*np.cos(d)/np.pow(m1,4.)

        a=(3.*q-p*p)/3.
        b=(2.*p*p*p-9.*p*q+27.*r)/27.

        test=b*b/4.+a*a*a/27.

        if (test>0.0):
            return -1.0
        elif (test==0.0):
          x1=np.sqrt(-a/3.)
          x2=x1
          x3=2.*x1
          if(b>0.0):
            x1*=-1.
            x2*=-1.
            x3*=-1.

        if(test<0.0):
          phi=np.acos(np.sqrt(-27.*b*b/4./a/a/a))
          x1=2.*np.sqrt(-a/3.)*np.cos(phi/3.)
          x2=2.*np.sqrt(-a/3.)*np.cos(phi/3.+np.pi*2./3.)
          x3=2.*np.sqrt(-a/3.)*np.cos(phi/3.+np.pi*4./3.)
          if(b>0.0):
            x1*=-1.
            x2*=-1.
            x3*=-1.

        s1=x1-p/3.
        s2=x2-p/3.
        s3=x3-p/3.

        if(s1<s2 and s1<s3):
          t1=s2
          t2=s3
        elif(s2<s1 and s2<s3):
          t1=s1
          t2=s3
        else:
          t1=s1
          t2=s2

        b1=np.asin(np.sqrt(t1))
        b2=np.asin(np.sqrt(t2))

        betas=b1
        betaw=b2
        if(b2>b1):
          betas=b2
          betaw=b1

        if(i==0):
            return betaw
        if(i==1):
            return betas

 def angle_to(self, other):
     """computes the angle between two vectors
         cos theta = (n * m) / (n.length * m.length)
     """
     arc = self.dot(other) / self.length / other.length
     if abs(arc - 1) <= 1e-6 or abs(arc + 1) <= 1e-6:
         arc = 1
     return np.acos(arc)

  def evaluate(self, xOrig):
    """
      Calculates and returns RV curve according to current model parameters.
      
      .. note:: The units of the model RV curve are **stellar-radii per second**.
      
      Parameters
      ----------
      xOrig : array
          The time stamps at which to calculate the model RV curve.
          Note that the orbit period and central transit time are used
          to convert time into "true anomaly".
    """
    x = self.trueAnomaly(xOrig)
    Xp = self.Xp(x)
    Zp = self.Zp(x)
    rho = self.rho(Xp, Zp)
    etap = self.etap(Xp, Zp)
    zeta = self.zeta(etap)
    x0 = self.x0(etap)
    xc = self.xc(zeta, x0)

    # dphase is the phase difference between the primary transit and the time points
    # It is used to exclude the secondary transit from the calculations 
    dphase = numpy.abs((xOrig-self["T0"])/self["P"])
    dphase = numpy.minimum( dphase-numpy.floor(dphase), numpy.abs(dphase-numpy.floor(dphase)-1))

    y = numpy.zeros(len(x))
    indi = numpy.where(numpy.logical_and(rho < (1.0-self["gamma"]), dphase < 0.25))[0]

    y[indi] = Xp[indi]*self["Omega"]*sin(self["Is"])* self["gamma"]**2 * \
        (1.0 - self["epsilon"]*(1.0 - self.W2(rho[indi]))) / \
        (1.0 - self["gamma"]**2 - self["epsilon"]*(1./3. - self["gamma"]**2*(1.0-self.W1(rho[indi]))))
    
    indi = numpy.where(numpy.logical_and( \
                numpy.logical_and(rho >= 1.-self["gamma"], rho < 1.0+self["gamma"]), dphase < 0.25))[0]
    z0 = self.z0(etap, indi)
    
    y[indi] = (Xp[indi]*self["Omega"]*sin(self["Is"])*( \
        (1.0-self["epsilon"]) * (-z0[indi]*zeta[indi] + self["gamma"]**2*acos(zeta[indi]/self["gamma"])) + \
        (self["epsilon"]/(1.0+etap[indi]))*self.W4(x0[indi], zeta[indi], xc[indi], etap[indi]))) / \
        (pi*(1.-1.0/3.0*self["epsilon"]) - (1.0-self["epsilon"]) * (asin(z0[indi])-(1.+etap[indi])*z0[indi] + \
        self["gamma"]**2*acos(zeta[indi]/self["gamma"])) - self["epsilon"]*self.W3(x0[indi], zeta[indi], xc[indi], etap[indi]))

    return y

def SunPositionLocalNoon ( DoY, longitude, latitude ):
    """
    Original code provided by: Jose Luis Gomez-Dans
    https://gist.github.com/733741

    Modified by: Gerardo Lopez-Saldana 

    Calculates the position of the sun given a position at local solar noon.
    Basically, all you need to know is here: 
    http://answers.google.com/answers/threadview/id/782886.html
    """
    import time
    from numpy import sin
    from numpy import cos
    from numpy import degrees
    from numpy import radians
    from numpy import arccos as acos
    from numpy import deg2rad

    #from math import sin, cos, degrees, radians, acos

    #DoY = int ( time.strftime( "%j", time.strptime( date, "%Y-%m-%d" ) ) )
    latitude = deg2rad ( latitude )
    #longitude = deg2rad ( longitude )

    # Calculate Local Solar Time LST
    n = numpy.round(DoY - 2451545 - 0.0009 - (longitude / 360.))
    J = 2451545. + 0.0009 + (longitude / 360.) + n
    J = (J - JulianDay) * 60

    #EoT
    M = (2*numpy.pi*JulianDay)/365.242
    EoT = -7.655 * numpy.sin(M) + 9.873*numpy.sin(2*M + 3.588)

    TimeZone = 0
    LST = (720-4*longitude-EoT+TimeZone*60)/1440

    longitude = deg2rad ( longitude )

    #LST_hours = int(LST*24.0)
    #LST_mins = int((LST*24.0 - LST_hours) * 60.)
    #LST_secs = (((LST*24.0 - LST_hours) * 60.) - LST_mins) * 60.
    #print "LST:", str(LST_hours), str(LST_mins), str(LST_secs)

    #( hh, mm )  = hour.split(":")
    #h = float(hh)
    #h = h + float(mm)/60.

    # To emulate MODIS products, set fixed LST = 12.00
    LST[:,:] = 12.0

    ##Now we can calculate the Sun Zenith Angle (SZA):
    h = (12.0 - (LST)) / 12.0 * numpy.pi
    delta = -23.45 * (numpy.pi/180.0) * cos (2 * numpy.pi/365.0 * (DoY+10))
    SZA = degrees( acos(sin(latitude) * sin(delta) + cos(latitude) * cos(delta) * cos(h)) )

    return (SZA, LST*24)

def acos(*args, **kw):
    arg0 = args[0]
    if isinstance(arg0, (int, float, long)):
        return _math.acos(*args,**kw)
    elif isinstance(arg0, complex):
        return _cmath.acos(*args,**kw)
    elif isinstance(arg0, _sympy.Basic):
        return _sympy.acos(*args,**kw)
    else:
        return _numpy.acos(*args,**kw)

def angle2vecs(vec1, vec2):
    """angle between two vectors"""
    # vector a * vector b = |a|*|b|* cos(angle between vector a and vector b)
    dot = np.dot(vec1, vec2)
    vec1_modulus = np.sqrt(np.multiply(vec1, vec1).sum())
    vec2_modulus = np.sqrt(np.multiply(vec2, vec2).sum())
    if (vec1_modulus * vec2_modulus) == 0:
        cos_angle = 1
    else: cos_angle = dot / (vec1_modulus * vec2_modulus)
    return math.degrees(acos(cos_angle))

    def drawFaces(self, level):
        glPushMatrix()
        glRotatef(180, 0, 0, 1) # let these be the upper face of pyramid
        glRotatef(acos(1.0/3.0)*(180.0/pi), 1,0,0) # apply the dihedral angle
        self.drawKochPyramid(level - 1)
        glPopMatrix()

        glPushMatrix()
        glTranslate(-1,0,0)
        glRotatef(300, 0, 0, 1) # let these be the upper face of pyramid
        glRotatef(acos(1.0/3.0)*(180.0/pi), 1,0,0) # apply the dihedral angle
        self.drawKochPyramid(level - 1)
        glPopMatrix()

        glPushMatrix()
        glTranslate(-1/2.0, -sqrt(3)/2.0, 0)
        glRotatef(420, 0, 0, 1) # let these be the lower right face of pyramid
        glRotatef(acos(1.0/3.0)*(180.0/pi), 1,0,0) # apply the dihedral angle
        self.drawKochPyramid(level - 1)
        glPopMatrix()

def angle(v1,v2):
    """
    calculates the angle between two vectors.
    v1 and v2 are numpy.array objects.
    returns a float containing the angle in radians.
    """
    length_product = norm(v1)*norm(v2)
    if length_product == 0:
        raise AngleGeometryError(\
        "Cannot calculate angle for vectors with length zero")
    angle = acos(scalar(v1,v2)/length_product)
    return angle

def angle(v1, v2):
    """ computes the angle between two vectors
        assuming carthesian coordinates

    Parameters
    ----------
    vec1 : numpy array
    vec2 : numpy array

    Returns
    -------
    the angle between vec1 and vec2 as a dimensioned astropy quantity
    """
    return acos(np.clip(v1.dot(v2) / (length(v1) * length(v2)), -1.0, 1.0))

def quaternion_to_axisangle(quat):
    data = quat._data
    sqr_length = sum(np.square(data[1:]))
    axis = np.zeros(3)
    if sqr_length > 0:
        angle = 2.0 * np.acos(data[0])
        inv_length = 1.0 / np.sqrt(sqr_length)
        axis = inv_length * data[1:]
    else:
        angle = 0
        axis[0] = 1
        axis[1] = 0
        axis[2] = 0
    return axis, angle

def aitoff(lon, lat):
    """
    Make Aitoff map projection.

    Take traditional longitude and latitude in radians and return a
    tuple (x, y).

    Notice that traditionally longitude is in [-pi:pi] from the meridian,
    and latitude is in [-pi/2:pi/2] from the equator. So, for example, if
    you would like to make a galactic map projection centered on the galactic
    center, before passing galactic longitude l to the function you should
    first do:
    l = l if l <= numpy.pi else l - 2 * numpy.pi

    Keyword arguments:
    lon -- Traditional longitude in radians, in range [-pi:pi]
    lat -- Traditional latitude in radians, in range [-pi/2:pi/2]
    """

    def sinc(x):
        # a quick unnormalized sinc function, with discontinuity removed
        if not x:
            return 0
        else:
            return numpy.sin(x) / x

    x = numpy.zeros_like(lon)
    y = numpy.zeros_like(lat)

    # check if the input values are in the range
    if lon > numpy.pi or lon < -numpy.pi or lat > numpy.pi / 2 or \
            lat < -numpy.pi / 2:
        print('Aitoff: Input longitude and latitude out of range.\n')
        print('           lon: [-pi,pi]; lat: [-pi/2,pi/2].\n')
        return None

    # take care of the sigularity at (0, 0), otherwise division by zero may
    # happen
    if lon == 0 and lat == 0:
        return 0.0, 0.0

    alpha = numpy.acos(numpy.cos(lat) * numpy.cos(lon / 2.0))

    # the sinc function used here is the unnormalized sinc function
    x = 2.0 * numpy.cos(lat) * numpy.sin(lon / 2.0) / sinc(alpha)

    y = numpy.sin(lat) / sinc(alpha)

    return x, y

 def generateAngle(self, nnew, emitted_energy, incident_energy):
     """
     Use Moller approximation, or delegate calculation to rsoopic (return
     result from call to rsoopic function h2crosssections.generateAngle)
     """
     if self.useMollerApproximation:
         theta = np.empty((nnew))
         for i in range(nnew):
             costheta = emitted_energy[i] * (incident_energy[i] + 2. * self.emassEV)
             costheta /= incident_energy[i] * (emitted_energy[i] + 2. * self.emassEV)
             costheta = np.sqrt(costheta)
             theta[i] = np.acos(costheta)
         return theta
     else:
         return h2crosssections.generateAngle(nnew, emitted_energy, incident_energy)