def _ang2vert_eq(nside,theta,phi,z):
		a= 4./3./nside
		b= 8./3./sc.pi
		deltaZ= a
		deltaPhi= a/b
		out= []
		out.append([sc.arccos(z+deltaZ/2),phi])
		out.append([theta,phi-deltaPhi/2.])
		out.append([sc.arccos(z-deltaZ/2),phi])
		out.append([theta,phi+deltaPhi/2.])
		return sc.array(out)

	def _xsys2ang(xs,ys):
		if sc.fabs(ys) <= sc.pi/4.:
			return [sc.arccos(8./3./sc.pi*ys),xs]
		else:
			xt= (xs % (sc.pi/2.))
			fabsys= sc.fabs(ys)
			theta= sc.arccos((1.-1./3.*(2.-4.*fabsys/sc.pi)**2.)*ys/fabsys)
			if fabsys == sc.pi/2.:
				phi= xs-fabsys+sc.pi/4.
			else:
				phi= xs-(fabsys-sc.pi/4.)/(fabsys-sc.pi/2.)*(xt-sc.pi/4.) 
			return [theta % (sc.pi+0.0000000001),phi % (2.*sc.pi)] #Hack

def W_Anom(planet,t=0,t0=0):
    """Wahre Anomalie"""
    E = ex_Anom(planet,t,t0)
    eps = Planet[planet]["numEx"] #Zur besseren Lesbarkeit
    
    if (E >= 0 and E <= sc.pi):
        return sc.arccos((sc.cos(E) - eps) / (1 - eps*sc.cos(E)))
        
    elif (E >= sc.pi and E <= 2*sc.pi):
        return 2*sc.pi - sc.arccos((sc.cos(E) - eps) / (1 - eps*sc.cos(E)))
            
    else: 
        return 0 #Außerhalb der Intervalle nicht definiert

def rext_calc(df, lat=float):
    """
    Function to calculate extraterrestrial radiation output in J/m2/day
    Ref:http://www.fao.org/docrep/x0490e/x0490e07.htm

    :param df: dataframe with datetime index
    :param lat: latitude (negative for Southern hemisphere)
    :return: Rext (J/m2)
    """
    # set solar constant [MJ m^-2 min^-1]
    s = 0.08166
    #convert latitude [degrees] to radians
    latrad = lat*math.pi / 180.0
    #have to add in function for calculating single value here
    # extract date, month, year from index
    date = pd.DatetimeIndex(df.index).day
    month = pd.DatetimeIndex(df.index).month
    year = pd.DatetimeIndex(df.index).year
    doy = met.date2doy(dd=date, mm=month, yyyy=year)  # create day of year(1-366) acc to date
    l = sp.size(doy)
    if l < 2:
        dt = 0.409 * math.sin(2 * math.pi / 365 * doy - 1.39)
        ws = sp.arccos(-math.tan(latrad) * math.tan(dt))
        j = 2 * math.pi / 365.25 * doy
        dr = 1.0 + 0.03344 * math.cos(j - 0.048869)
        rext = s * 1440 / math.pi * dr * (ws * math.sin(latrad) * math.sin(dt) + math.sin(ws) * math.cos(latrad) * math.cos(dt))
    #Create dummy output arrays sp refers to scipy
    else:
        rext = sp.zeros(l)
        dt = sp.zeros(l)
        ws = sp.zeros(l)
        j = sp.zeros(l)
        dr = sp.zeros(l)
        #calculate Rext
        for i in range(0, l):
            #Calculate solar decimation dt(d in FAO) [rad]
            dt[i] = 0.409 * math.sin(2 * math.pi / 365 * doy[i] - 1.39)
            #calculate sunset hour angle [rad]
            ws[i] = sp.arccos(-math.tan(latrad) * math.tan(dt[i]))
            # calculate day angle j [radians]
            j[i] = 2 * math.pi / 365.25 * doy[i]
            # calculate relative distance to sun
            dr[i] = 1.0 + 0.03344 * math.cos(j[i] - 0.048869)
            #calculate Rext dt = d(FAO) and latrad = j(FAO)
            rext[i] = (s * 1440.0 / math.pi) * dr[i] * (ws[i] * math.sin(latrad) * math.sin(dt[i]) + math.sin(ws[i])* math.cos(latrad) * math.cos(dt[i]))

    rext = sp.array(rext) * 1000000
    return rext

def add_geometric_edge_properties(G):
    """ Adds angle to cortical surface (in degrees), cortical depth, volume and
    cross section to each edge in the graph.
    INPUT: G:  Vascular graph in iGraph format.
    OUTPUT: None - the vascular graph G is modified in place.
    """
           
    depth = []
    angle = []
    crossSection = [] 
    volume = []

    ez = array([0,0,1])
    
    for edge in G.es:         
        
        a = G.vs[edge.source]['r']
        b = G.vs[edge.target]['r']
        v = a-b    
        depth.append((a[2]+b[2])/2.0)
        
        theta=arccos(dot(v,ez)/norm(v))/2/pi*360
        if theta > 90:
            theta = 180-theta
        angle.append(theta)
    
        crossSection.append(np.pi * edge['diameter']**2. / 4.)
        volume.append(crossSection[-1] * edge['length'])
    
    
    G.es['depth'] = depth
    G.es['angle'] = angle
    G.es['volume'] = volume
    G.es['crossSection'] = crossSection 

def create_map(parameters,thetas,phis,vrs):
    
    
    pix = hp.ang2pix(parameters.nside,thetas,phis)
    number_of_pixels = hp.nside2npix(parameters.nside)
    vrmap = sp.ones(number_of_pixels)*parameters.badval
    
#    pdb.set_trace()
    vrs = sp.array(vrs)
    vrs_mean_of_repeated_pixels = copy.copy(vrs)
    for p in set(pix):
        vrs_mean_of_repeated_pixels[pix == p] = sp.mean(vrs[pix == p])
    
    vrmap[pix] = vrs_mean_of_repeated_pixels

#    pdb.set_trace()    
    theta_max = sp.arccos(1-2*parameters.skyfraction)
    pix_all = sp.array(range(number_of_pixels))
    pix_unseen = pix_all[hp.pix2ang(parameters.nside,pix_all)[0]>theta_max]
    vrmap[pix_unseen] = parameters.unseen

    empty_pixels = number_of_pixels-len(set(pix))
    print "The number of empty pixels inside the survey is", empty_pixels
    print "The corresponds to a fraction of", empty_pixels/number_of_pixels
    print "The number of pixels outside survey is", len(pix_unseen)
#    pdb.set_trace()
    return vrmap

def overlap(x, xdown, y, ydown, z, zdown, r, rdown, alpha):
    overlap_fraction = np.zeros(np.size(x))
    for i in range(0, np.size(x)):
        #define dx as the upstream x coordinate - the downstream x coordinate then rotate according to wind direction
        dx = xdown - x[i]
        #define dy as the upstream y coordinate - the downstream y coordinate then rotate according to wind direction
        dy = abs(ydown - y[i])
        dz = abs(zdown - z[i])
        d = sp.sqrt(dy**2.+dz**2.)
        R = r[i]+dx*alpha #The radius of the wake depending how far it is from the turbine
        A = rdown**2*pi #The area of the turbine
        if dx > 0:
            #if d <= R-rdown:
            #    overlap_fraction[i] = 1 #if the turbine is completely in the wake, overlap is 1, or 100%
            if d == 0:
                print "Area of turbine: ", A
                print "Area of wake: ", pi*R**2
                if A <= pi*R**2:
                    overlap_fraction[i] = 1.
                else: 
                    overlap_fraction[i] = pi*R**2/A
            elif d >= R+rdown:
                overlap_fraction[i] = 0 #if none of it touches the wake, the overlap is 0
            else:
                #if part is in and part is out of the wake, the overlap fraction is defied by the overlap area/rotor area
                overlap_area = rdown**2.*sp.arccos((d**2.+rdown**2.-R**2.)/(2.0*d*rdown))+R**2.*sp.arccos((d**2.+R**2.-rdown**2.)/(2.0*d*R))-0.5*sp.sqrt((-d+rdown+R)*(d+rdown-R)*(d-rdown+R)*(d+rdown+R))
                overlap_fraction[i] = overlap_area/A
        else:
            overlap_fraction[i] = 0 #turbines cannot be affected by any wakes that start downstream from them

    # print overlap_fraction
    print overlap_fraction
    return overlap_fraction #retrun the n x n matrix of how each turbine is affected by all of the others

def Rz_to_uv(R,z,delta=1.):
    """
    NAME:

       Rz_to_uv

    PURPOSE:

       calculate prolate confocal u and v coordinates from R,z, and delta

    INPUT:

       R - radius

       z - height

       delta= focus

    OUTPUT:

       (u,v)

    HISTORY:

       2012-11-27 - Written - Bovy (IAS)

    """
    coshu, cosv= Rz_to_coshucosv(R,z,delta)
    u= sc.arccosh(coshu)
    v= sc.arccos(cosv)
    return (u,v)

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 __new__(self, nbinaries=1e6):
     arr = sp.ones(nbinaries, dtype=[(name, 'f8') for name in ['period', 'mass_ratio', 'eccentricity', 'phase', 'theta', 'inclination']])
     arr['eccentricity'] = 0.
     arr['phase'] = sp.rand(nbinaries)
     arr['theta'] = sp.rand(nbinaries) * 2 * sp.pi
     arr['inclination'] = sp.arccos(sp.rand(nbinaries) * 2. - 1.)
     return arr.view(OrbitalParameters)

def resolve_tri(A,B,a,b,up=True):
    AB=A-B
    c=l.norm(AB)
    aa=s.arctan2(AB[1],AB[0])
    bb=s.arccos((b**2+c**2-a**2)/(2*b*c))
    if up: return B+b*s.array((s.cos(aa+bb),s.sin(aa+bb)))
    else: return B+b*s.array((s.cos(aa-bb),s.sin(aa-bb)))

def dotties():
    data = fi.read_data("../sgrnorth_paper/planefit_results.txt", ",")
    planes = []
    for i in range(6):
        planes.append(data[i,:])
    #for plane in planes:  print plane
    dots = sc.zeros((6,6))
    for i in range(6):
        for j in range(6):
            dots[i,j] = (sc.arccos(planes[i][0]*planes[j][0] + planes[i][1]*planes[j][1]
                                   + planes[i][2]*planes[j][2]))*deg
    print dots
    print
    errors = sc.zeros((6,6))
    for i in range(6):
        for j in range(6):
            dotty = planes[i][0]*planes[j][0]+planes[i][1]*planes[j][1]+planes[i][2]*planes[j][2]
            factor = -1.0 / sc.sqrt(1 - dotty*dotty)
            print factor
            errors[i,j] = factor*(planes[j][0]*planes[i][4] + planes[i][0]*planes[j][4] +
                                  planes[j][1]*planes[i][5] + planes[i][1]*planes[j][5] +
                                  planes[j][2]*planes[i][6] + planes[i][2]*planes[j][6])*deg
            #if i==3 and j==5:  print dotty
    print errors
    Sgr = [14.655645800014774, 2.2704309216189817, -5.8873772610912614]
    print "# - dist to Sgr Core:"
    for plane in planes:
        print (Sgr[0]*plane[0] + Sgr[1]*plane[1] + Sgr[2]*plane[2] + plane[3]), \
        (Sgr[0]*plane[4] + Sgr[1]*plane[5] + Sgr[2]*plane[6] + plane[7])

def besselFourierKernel(m, zero, rho):
    """ Function kernel for the bessel Fourier method inversion
    
    Uses the mathematical formulation laid out in L. Wang and R. Granetz,
    Review of Scientific Instruments, 62, p.842, 1991. This generates
    the necessary weighting for a given chord in bessel/fourier space
    for a given tangency radius rho. There should be little reason to
    use this function unless modified and generating a new inversion
    scheme utilizing this kernel.

    Args:
        m: geometry Object with reference origin

        zero: geometry Object with reference origin
        
        rho: normalized tangency radius

    Returns:
        numpy array: Vector points from pt1 to pt2.
    
    """


    # I SHOULD TRY AND VECTORIZE THIS AS MUCH AS POSSIBLE
    jprime = (scipy.special.jn(m+1, zero) - scipy.special.jn(m-1, zero))
    return jprime*scipy.integrate.quad(_beam.bessel_fourier_kernel,
                                       0,
                                       scipy.arccos(rho),
                                       args = (m, zero, rho))[0]

def young(phase,
          sigma_sg,
          sigma_sl,
          surface_tension='pore.surface_tension',
          **kwargs):
    r'''
    Calculate contact angle using Young's equation
    
    Notes
    -----
    Young's equation is: sigma_lg*Cos(theta)=sigma_sg - sigma_sl
    where
    sigma_lg is the liquid-gas surface tension [N/m]
    sigma_sg is the solid-gas surface tension [N/m]
    sigma_sl is the solid-liquid interfacial tension [J/m^2]
    theta is the Young contact angle [rad]
           
    '''
    if surface_tension.split('.')[0] == 'pore':
        sigma = phase[surface_tension]
        sigma = phase.interpolate_data(data=sigma)
    else:
        sigma = phase[surface_tension]
    theta = sp.arccos((sigma_sg - sigma_sl)/phase[surface_tension])
    theta = sp.rad2deg(theta)
    return theta

def pos2Ray(pos, tokamak, angle=None, eps=1e-6):
    r"""Take in GENIE pos vectors and convert it into TRIPPy rays
    
    Args:
        pos: 4 element tuple or 4x scipy-array
            Each pos is assembled into points of (R1,Z1,RT,phi)

        tokamak: 
            Tokamak object in which the pos vectors are defined.
            
    Returns:
        Ray: Ray object or typle of ray objects.
        
    """

    r1 = scipy.array(pos[0])
    z1 = scipy.array(pos[1])
    rt = scipy.array(pos[2])
    phi = scipy.array(pos[3])

    zt = z1 - scipy.tan(phi)*scipy.sqrt(r1**2 - rt**2)
    angle2  = scipy.arccos(rt/r1)

    if angle is None:
        angle = scipy.zeros(r1.shape)

    pt1 = geometry.Point((r1,angle,z1),tokamak)
    pt2 = geometry.Point((rt,angle+angle2,zt),tokamak)

    output = Ray(pt1,pt2)
    output.norm.s[-1] = eps
    tokamak.trace(output)
    return output

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 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 errorAngle(actual,estimate):        
     error = actual*estimate.conjugate()
     error.normalise()
     # Cast error to float so that we can deal with fixed point numbers
     # throws an error otherwise
     eAngle = degrees(2*arccos(float(error.w)))
     
     eAngle = abs(eAngle - 360) if eAngle > 180 else eAngle
     return eAngle

def guess_initial_parameters(mean, phi0, omega):
    mean_max = mean.max()
    mean_min = mean.min()
    offset = (mean_max + mean_min)/2.
    amplitude = (mean_max - mean_min)/2.
    if phi0 is None or omega is None:
        y = mean-offset
        y2 = savgol_filter(y, 11, 1, mode="nearest")
        if phi0 is None:
            cos_phi0 = clip(y2[0]/amplitude, -1., 1.)
            if y2[1] > y2[0]:
                phi0 = -arccos(cos_phi0)
            else:
                phi0 = arccos(cos_phi0)
        if omega is None:
            zero_crossings = scipy.where(scipy.diff(scipy.sign(y2)))[0]
            omega = pi/scipy.average(scipy.diff(zero_crossings))
    return offset, amplitude, phi0, omega

def aoa_diff_rad(aoa_a,aoa_b):
    """
    Calculate the difference between two angles of arrival, in radians.
    """
    # Prepend radius 1 and convert to cartesian.
    cart_a = sph2cart( (1,) + tuple(aoa_a) )
    cart_b = sph2cart( (1,) + tuple(aoa_b) )
    # Property of dot product between vectors.
    return sp.arccos(sp.dot(cart_a,cart_b))