def hertz_to_bark(cf):
    """
    ::

        Convert frequency in Hz to Bark band 
    """
    return 13 * np.arctan(0.00076 * cf) + 3.5 * np.arctan( (cf / 7500.)**2)

def calc_theta(x_predicted, y_predicted, x_true, y_true, x_ref, y_ref):
    """ Calculate the angle the predicted position and the true position, where the zero degree corresponds to the line joing the true halo position and the reference point given.
    Arguments:
        x_predicted, y_predicted: vector for predicted x- and y-positions (1 to 3 elements)
        x_true, y_true: vector for known x- and y-positions (1 to 3 elements)
        Note that the input of these are matched up so that the first elements of each
        vector are associated with one another
        x_ref, y_ref: scalars of the x,y coordinate of reference point
    Returns:
        Theta: A vector containing the angles of the predicted halo w.r.t the true halo
        with the vector joining the reference point and the halo as the zero line. 
    """

    num_halos=len(x_predicted)
    theta=np.zeros([num_halos+1],float) #Set up the array which will pass back the values
    psi = np.arctan( (y_true-y_ref)/(x_true-x_ref) ) # Angle at which the halo is at
                                                     #with respect to the reference poitn
    phi = np.arctan((y_predicted-y_true)/(x_predicted-x_true)) # Angle of the estimate
                                                               #wrt true halo centre

    #Before finding the angle with the zero line as the line joiing the halo and the reference
    #point I need to convert the angle produced by Python to an angle between 0 and 2pi
    phi =convert_to_360(phi, x_predicted-x_true,\
         y_predicted-y_true)
    psi = convert_to_360(psi, x_true-x_ref,\
                             y_true-y_ref)
    theta = phi-psi #The angle with the baseline as the line joing the ref and the halo

    
    theta[theta< 0.0]=theta[theta< 0.0]+2.0*mt.pi #If the angle of the true pos wrt the ref is
                                                  #greater than the angle of predicted pos
                                                  #and the true pos then add 2pi
    return theta

 def _fgreen3d(self, z, y, x):
     ''' Return the periodic integrated greens funcion on the 'original'
     domain
     Qiang, Lidia, Ryne,Limborg-Deprey, PRSTAB 10, 129901 (2007)
     Args:
         x,y,z: arrays, e.g. x, y, z = np.meshgrid(xx, yy, zz)
     '''
     abs_r = np.sqrt(x * x + y * y + z * z)
     inv_abs_r = 1./abs_r
     tmpfgreen =  (-(  +    z*z * np.arctan(x*y*inv_abs_r/z)
                   +   y*y * np.arctan(x*z*inv_abs_r/y)
                   +   x*x * np.arctan(y*z*inv_abs_r/x)
                )/2.
                 + y*z*np.log(x+abs_r)
                 + x*z*np.log(y+abs_r)
                 + x*y*np.log(z+abs_r))
     fgreen = np.zeros((2 * self.mesh.nz,
                        2 * self.mesh.ny,
                        2 * self.mesh.nx), dtype=np.complex128)
     # evaluate the indefinite integral per cell (int_a^b f = F(b) - F(a))
     fgreen[:self.mesh.nz, :self.mesh.ny, :self.mesh.nx] = (
              tmpfgreen[ 1:,  1:,  1:]
             -tmpfgreen[:-1,  1:,  1:]
             -tmpfgreen[ 1:, :-1,  1:]
             +tmpfgreen[:-1, :-1,  1:]
             -tmpfgreen[ 1:,  1:, :-1]
             +tmpfgreen[:-1,  1:, :-1]
             +tmpfgreen[ 1:, :-1, :-1]
             -tmpfgreen[:-1, :-1, :-1]
             ) * 1./self.mesh.volume_elem # divide by vol_elem to average!
     return fgreen

def _compute_static_prob(tri, com):
    """
    For an object with the given center of mass, compute
    the probability that the given tri would be the first to hit the
    ground if the object were dropped with a pose chosen uniformly at random.

    Parameters
    ----------
    tri: (3,3) float, the vertices of a triangle
    cm:  (3,) float, the center of mass of the object

    Returns
    -------
    prob: float, the probability in [0,1] for the given triangle
    """
    sv = [(v - com) / np.linalg.norm(v - com) for v in tri]

    # Use L'Huilier's Formula to compute spherical area
    a = np.arccos(min(1, max(-1, np.dot(sv[0], sv[1]))))
    b = np.arccos(min(1, max(-1, np.dot(sv[1], sv[2]))))
    c = np.arccos(min(1, max(-1, np.dot(sv[2], sv[0]))))
    s = (a + b + c) / 2.0

    # Prevents weirdness with arctan
    try:
        return 1.0 / np.pi * np.arctan(np.sqrt(np.tan(s / 2) * np.tan(
            (s - a) / 2) * np.tan((s - b) / 2) * np.tan((s - c) / 2)))
    except BaseException:
        s = s + 1e-8
        return 1.0 / np.pi * np.arctan(np.sqrt(np.tan(s / 2) * np.tan(
            (s - a) / 2) * np.tan((s - b) / 2) * np.tan((s - c) / 2)))

def expected(scheme, angle_degrees):
    angle = angle_degrees * np.pi / 180.0

    cohesion = 10
    friction_degrees = 20
    tip_smoother = 4
    mean = -10
    friction = friction_degrees * np.pi / 180.0

    if (scheme == "native"):
        coh = cohesion
        fric = friction
    elif (scheme == "outer_tip"):
        coh = 2 * np.sqrt(3) * cohesion * np.cos(friction) / (3.0 - np.sin(friction))
        fric = np.arctan(2 * np.sin(friction) / np.sqrt(3) / (3.0 - np.sin(friction)))
    elif (scheme == "inner_tip"):
        coh = 2 * np.sqrt(3) * cohesion * np.cos(friction) / (3.0 + np.sin(friction))
        fric = np.arctan(2 * np.sin(friction) / np.sqrt(3) / (3.0 + np.sin(friction)))
    elif (scheme == "lode_zero"):
        coh = cohesion * np.cos(friction)
        fric = np.arctan(np.sin(friction) / 3.0)
    elif (scheme == "inner_edge"):
        coh = 3 * cohesion * np.cos(friction) / np.sqrt(9.0 + 3.0 * np.power(np.sin(friction), 2))
        fric = np.arctan(np.sin(friction) / np.sqrt(9.0 + 3.0 * np.power(np.sin(friction), 2)))

    bar = np.sqrt(np.power(coh - mean * 3.0 * np.tan(fric), 2) - np.power(tip_smoother, 2))
    x = bar * np.cos(angle)
    y = bar * np.sin(angle)
    return (x, y)

    def deflections(self, xin, yin):
        from numpy import arctanh, arctan, arctan2, log, sin, cos

        #x,y = self.align_coords(xin,yin)
        x = xin - self.x
        y = yin - self.y
        b, rs = self.b, self.rs
        X = b / rs
        if X < 1.:
            amp = X**2 / (8 * arctanh(
                ((1 - X) / (1 + X))**0.5) / (1 - X**2)**0.5 + 4 * log(X / 2.))
        elif X == 1:
            amp = 0.25 / (1. + log(0.5))
        else:
            amp = X**2 / (8 * arctan(
                ((X - 1) / (1 + X))**0.5) / (X**2 - 1)**0.5 + 4 * log(X / 2.))

        r2 = (x**2 + y**2) / rs**2
        r = r2**0.5
        F = r * 0.
        F[r < 1.] = arctanh((1 - r2[r < 1.])**0.5) / (1 - r2[r < 1.])**0.5
        F[r == 1.] = 1.
        F[r > 1.] = arctan((r2[r > 1.] - 1.)**0.5) / (r2[r > 1.] - 1)**0.5

        dr = 4 * amp * rs * (log(r / 2) + F) / r
        A = arctan2(y, x)
        return dr * cos(A), dr * sin(A)

def rotMatrixfromXYZ(station,mode='LBA'):
    """Return a rotation matrix which will rotate a station to (0,0,1)"""
    loc=station.antField.location[mode]
    longRotMat=rotationMatrix(0.,0.,-1.*n.arctan(loc[1]/loc[0]))
    loc0=n.dot(longRotMat,loc)
    latRotMat=rotationMatrix(0.,n.arctan(loc0[0,2]/loc0[0,0]),0.)
    return n.dot(latRotMat,longRotMat)

    def remap(self):
        """
        warp the linear pixel coordinates to a spherical corrected representation.

        Function is called when the monitor object is initialized and populate
        the `deg_coord_x` and `deg_coord_y` attributes.
        """

        resolution = [0, 0]
        resolution[0] = self.resolution[0] / self.downsample_rate
        resolution[1] = self.resolution[1] / self.downsample_rate

        map_coord_x, map_coord_y = np.meshgrid(range(resolution[1]),
                                               range(resolution[0]))

        new_map_x = np.zeros(resolution, dtype=np.float32)
        new_map_y = np.zeros(resolution, dtype=np.float32)

        for j in range(resolution[1]):
            new_map_x[:, j] = ((180.0 / np.pi) *
                               np.arctan(self.lin_coord_x[0, j] / self.dis))
            dis2 = np.sqrt(np.square(self.dis) +
                           np.square(self.lin_coord_x[0, j]))

            for i in range(resolution[0]):
                new_map_y[i, j] = ((180.0 / np.pi) *
                                   np.arctan(self.lin_coord_y[i, 0] / dis2))

        self.deg_coord_x = new_map_x + self.center_coordinates[1]
        self.deg_coord_y = new_map_y + self.center_coordinates[0]

def evalExactφ(sqp,sqpnext,sqpdesired):

	a=sqp/sqpdesired
	b=sqpnext/sqpdesired
	if(abs(a+1) > float("1e-12") and (a*a + b*b - 1)>0):

		φ1=2*np.arctan(     ( b -  (a*a + b*b - 1)**0.5 )/(a+1)    )
		φ2=2*np.arctan(     ( b +  (a*a + b*b - 1)**0.5 )/(a+1)    )

		if(np.isnan(φ1)):
			print("Ok nan found, locating butter chicken")
			print( (a*a + b*b -1) )
			print( (b + (a*a + b*b -1)**0.5 ) )

		φ1=np.arcsin(np.sin(φ1))
		φ2=np.arcsin(np.sin(φ2))
		if(abs(φ1)<abs(φ2)):
			print(φ1,φ2)
			return φ1
		else:
			print(φ2,φ1)
			return φ2
	else:
		print("Glaba, something went globular :P")
		return 0

def extend_and_norm_feature(element, feature, command, num_elems):
    feature['element'] = element

    feature['tagname_edit'] = str_util.get_mean_distance_of_words(command, [feature['tagname']])

    # alt. str_util.get_normed_dist_for_words or str_util.get_mean_distance_of_words
    distance = str_util.new_dist
    feature['text_words_edit'] = distance(command, feature['text_words'])
    feature['sibling_text_words_edit'] = distance(command, feature['sibling_text_words'])

    feature['n_children'] = 1 - float(feature['n_children']) / num_elems

    w,h = feature['width'], feature['height']
    mx, my, sx, sy = 54.611, 25.206, 43.973, 6.467

    prior = 0.02
    likelihood_button, marginal = likelihood_and_marginal(w, h)

    if marginal != 0:
        feature['button_model'] = ((likelihood_button * prior) / marginal) - .147
    else:
        feature['button_model'] = 0

    # relative x and y can  be more than 1 because things can be beyond the edge of the window
    # so nudge things to be between -1 and 1
    feature['relative_x'] = np.arctan(1 * (feature['relative_x'] + 0.5)) / (np.pi / 2)
    feature['relative_y'] = np.arctan(1 * feature['relative_y']) / (np.pi / 2)

    # Feature ideas:
        # new on page?
        # position (relative to last action?)
        # color
        # relative tab index
        # text specificity
    return feature

def getPossibleCoordinates(fromPos, cvSize, a=None):
    a = DEFLECTIONBORDERLENGTH/2. if a is None else a
    counter = 0
    while(True):
        counter+=1
        if(counter==1000):
            print("Warning: needed 1000 attempts to find new coordinates for trajectory")
        newPos = [None, None]
        alpha = uniform(0, 2*np.pi)
        translationLength = uniform(TRANSLATIONLENGTH[0],TRANSLATIONLENGTH[1])
        
        if(a==0):
            x = math.floor(fromPos[0]+np.cos(alpha)*translationLength)
            y = math.floor(fromPos[1]-np.sin(alpha)*translationLength)
            return [x,y]
        else:
            if(fromPos[0] < cvSize[0]/2.):
                alpha1 = np.arctan(fromPos[0]/a)
                if(alpha < np.pi/2. + alpha1 or alpha > np.pi*3/2. - alpha1):
                    newPos[0] = math.floor(fromPos[0]+np.cos(alpha)*translationLength)
            else:
                alpha1 = np.arctan((cvSize[0]-fromPos[0])/a)
                if(alpha > np.pi/2. - alpha1 and alpha < np.pi*3/2. + alpha1):
                    newPos[0] = math.floor(fromPos[0]+np.cos(alpha)*translationLength)            
            if(fromPos[1]<cvSize[1]/2.):
                alpha2 = np.arctan(fromPos[1]/a)
                if(alpha < alpha2 or alpha > np.pi-alpha2):
                    newPos[1] = math.floor(fromPos[1]-np.sin(alpha)*translationLength)
            else:
                alpha2 = np.arctan((cvSize[1]-fromPos[1])/a) 
                if(alpha < np.pi+alpha2 or alpha > 2*np.pi - alpha2):
                    newPos[1] = math.floor(fromPos[1]-np.sin(alpha)*translationLength) 
            
            if(newPos[0] is not None and newPos[1] is not None and keepMiddlepointOnCanvas(cvSize, newPos)):
                return newPos

def kep_eqtnP(del_t, p, mu=Earth.mu):
    """Parabolic solution to Kepler's Equation (Algorithm 3)
    
    A trigonometric approach to solving Kepler's Equation for 
    parabolic orbits. For reference, see Algorithm 3 in Vallado (Fourth 
    Edition), Section 2.2 (pg 69).
    
    Parameters
    ----------
    del_t: double
        Change in time
    p: double
        Semi-parameter
    mu: double, optional, default = Earth.mu
        Gravitational parameter; defaults to Earth
        
    Returns
    -------
    B: double
        Parabolic Anomaly (radians)
    """
    
    p3 = p**3
    n_p = 2.0*np.sqrt(mu/p3)
    s = 0.5*np.arctan(2.0/(3.0*n_p*del_t))
    w = np.arctan((np.tan(s))**(1/3.0))
    B = 2.0/np.tan(2.0*w)
    return B

def calcMultipliers(meanM, stdM, Cfactors, phiShifts, coeffs, intercept, T):
    
    K = int(len(coeffs)/2)
    C = np.zeros(K)
    phi = np.zeros(K)
    for k in range(K):
        C[k] = np.sqrt(coeffs[2*k]**2 + coeffs[2*k+1]**2)
        if coeffs[2*k] > 0:
            phi[k] = np.arctan(coeffs[2*k+1]/coeffs[2*k])
        elif coeffs[2*k] < 0:
            phi[k] = np.arctan(coeffs[2*k+1]/coeffs[2*k]) + math.pi
        else:
            phi[k] = math.pi/2
        
    y1 = np.zeros(T)
    y2 = np.zeros(T)
    for t in range(T):
        y1[t] = y1[t] + intercept
        y2[t] = y2[t] + intercept
        for k in range(K):
            y1[t] = y1[t] + C[k]*np.cos((2*math.pi*(k+1)*(t+1)/T)-phi[k])
            y2[t] = y2[t] + C[k]*Cfactors[k]*np.cos((2*math.pi*(k+1)*(t+1)/T)-(phi[k]-phiShifts[k]))
    
    meanMultipliers = meanM * y2 / y1    
    stdMultipliers = stdM * np.ones(len(meanMultipliers))
    
    return meanMultipliers, stdMultipliers

def _filter_small_slopes(hgt, dx, min_slope=0):
    """Masks out slopes with NaN until the slope if all valid points is at 
    least min_slope (in degrees).
    """

    min_slope = np.deg2rad(min_slope)
    slope = np.arctan(-np.gradient(hgt, dx))  # beware the minus sign
    # slope at the end always OK
    slope[-1] = min_slope

    # Find the locs where it doesn't work and expand till we got everything
    slope_mask = np.where(slope >= min_slope, slope, np.NaN)
    r, nr = label(~np.isfinite(slope_mask))
    for objs in find_objects(r):
        obj = objs[0]
        i = 0
        while True:
            i += 1
            i0 = objs[0].start-i
            if i0 < 0:
                break
            ngap =  obj.stop - i0 - 1
            nhgt = hgt[[i0, obj.stop]]
            current_slope = np.arctan(-np.gradient(nhgt, ngap * dx))
            if i0 <= 0 or current_slope[0] >= min_slope:
                break
        slope_mask[i0:obj.stop] = np.NaN
    out = hgt.copy()
    out[~np.isfinite(slope_mask)] = np.NaN
    return out

def calcAngle(accel):
	angle = np.array([
					np.arctan(accel[0]/np.sqrt(accel[1]**2+accel[2]**2))+1.5707963267948966,
					np.arctan(accel[1]/np.sqrt(accel[0]**2+accel[2]**2))+1.5707963267948966,
					np.arctan(accel[2]/np.sqrt(accel[1]**2+accel[0]**2))+1.5707963267948966,
					])
	return angle

def test_vector(Simulator, nl):
    """A network that represents sin(t), cos(t), arctan(t)."""
    N = 40

    m = nengo.Network(label='test_vector', seed=123)
    with m:
        m.config[nengo.Ensemble].neuron_type = nl()
        input = nengo.Node(
            output=lambda t: [np.sin(t), np.cos(t), np.arctan(t)])
        A = nengo.Ensemble(N * 3, 3, radius=2)
        nengo.Connection(input, A)
        in_p = nengo.Probe(input, 'output')
        A_p = nengo.Probe(A, 'decoded_output', synapse=0.02)

    sim = Simulator(m)
    sim.run(5)

    with Plotter(Simulator, nl) as plt:
        t = sim.trange()
        plt.plot(t, sim.data[in_p], label='Input')
        plt.plot(t, sim.data[A_p], label='Neuron approximation, pstc=0.02')
        plt.legend(loc='best', prop={'size': 10})
        plt.savefig('test_ensemble.test_vector.pdf')
        plt.close()

    target = np.vstack((np.sin(np.arange(5000) / 1000.),
                        np.cos(np.arange(5000) / 1000.),
                        np.arctan(np.arange(5000) / 1000.))).T
    logger.debug("In RMSE: %f", npext.rmse(target, sim.data[in_p]))
    assert npext.rmse(target, sim.data[in_p]) < 0.01
    assert npext.rmse(target, sim.data[A_p]) < 0.1

def main():
    from model import Model
    #print(rotate([[2,0,0],[0,1,0],[0,0,1]], 90, 'y'))
    #print(calc_rot_array_from_hkl(19,-24,28))
    modelfile = sys.argv[1]
    m = Model(modelfile)
    #rot_arr = calc_rot_array_from_hkl(41,60,-6)
    rot_arr = calculate_rotation_array(30, 60, 90)
    rotate(m, rot_arr)
    #m.write_real_xyz('temp.real.xyz')
    return

    # Below is a (the?) rotation matrix of Pei's t1 that gives some planes. Oriented for a specific plane ~.
    #rot_arr = [ -0.977103, -0.123352, -0.173361, -0.130450, 0.990997, 0.030118, 0.168085, 0.052043, -0.984398 ]
    #rot(m,rot_arr)

    # Angles in radians
    # Note that these are semi difficult to figure out from the vesta rotation matrix,
    # partly because there are negative angles, so you may need to do 2pi - angle you found.
    #t1 = np.pi*2 - 0.0371505
    #t2 = 0.162790
    #t3 = 0
    #rot_arr = calc_rot_array(m,t1,t2,t3)
    #rot(m,rot_arr)

    kx = -0.76094085
    ky = 0.028182994
    kz = -0.648208872
    t2 = np.arctan(-ky/kx)
    t3 = 0.0
    t1 = np.arctan( (kx*np.cos(t2)-ky*cos(t2))/kz )
    t1 = 0.0
    print(t1,t2,t3)
    rot_arr = calc_rot_array(t1,t2,t3)
    rot(m,rot_arr)

    def exact(self, t):
        # Valid for linear s(u)
        k, b, m, A_F, w_F, I, V = self.k, self.b, self.m, \
                                  self.A_F, self.w_F, self.I, self.V
        b_crit = 2*np.sqrt(k*m)
        w_e = np.sqrt(k/m)
        zeta = b/b_crit
        zeta1p = zeta + np.sqrt(zeta**2 - 1)
        zeta1m = zeta - np.sqrt(zeta**2 - 1)
        zeta2 = np.sqrt(zeta**2 - 1)

        if zeta > 1:
            # No oscillations
            sol1 = (V + w_e*zeta1p*I)/(2*w_e*zeta2)*np.exp(-w_e*zeta1m*t)
            sol2 = (V + w_e*zeta1m*I)/(2*w_e*zeta2)*np.exp(-w_e*zeta1p*t)
            u_h = sol1 - sol2
        elif zeta == 1:
            u_h = (I + (V + w_e*I)*t)*np.exp(-w_e*t)
        else:
            # Oscillations
            A = np.sqrt(I**2 + ((V + zeta*w_e*I)/(w_e*zeta2))**2)
            phi = np.arctan((V + zeta*w_e*I)/(I*w_e*zeta2))
            u_h = A*np.exp(-zeta*w_e*t)*np.cos(zeta2*w_e*t - phi)

        # Excitation: F=F_0*cos(w_F*t)
        # F_0 and w_F must be extracted from F......?
        phi_0 = np.arctan(b*w_F/(k - m*w_F**2))
        A_0 = A_F/np.sqrt((k - m*w_F**2)**2 + (b*w_F)**2)
        u_p = A_0*np.cos(omega*t - phi_0)

        # Test: all special cases...
        return u_h + u_p

def place_label(x,y,label,indice=None,cotan=False,color='k'):
    """ Routine qui se débrouille pour mettre un label semi-transparent au 
    niveau de la courbe données par ses coordonnées x et y. Si on sait que le 
    label sera presque vertical avec possibilité de dépasser 90°, on peut 
    utiliser cotan=True pour corriger (considération purement esthétique). 
    'indice' correspond à la position dans les tableaux x et y où devra 
    s'afficher le label demandé. """
    print(x[0],y[0],label) # un peu de feedback pour savoir ce qu'on calcule
    N = len(x)//2          # Emplacement par défaut
    if indice: N=indice    # sauf si l'utilisateur impose la valeur
    xi,xf = plt.xlim()     # Les limites en x du graphe
    yi,yf = plt.ylim()     # Pareil en y
    Xsize = xf - xi        # La largeur
    # Pour la hauteur et la pente, cela dépend si les ordonnées sont en repère 
    # logarithmique ou non.
    if Plogscale:
        Ysize = np.log10(yf) - np.log10(yi)
        a = (np.log10(y[N+1])-np.log10(y[N-1]))/(x[N+1]-x[N-1]) * Xsize/Ysize
    else:
        Ysize = yf - yi
        a = (y[N+1]-y[N-1])/(x[N+1]-x[N-1]) * Xsize/Ysize
    bbox = plt.gca().get_window_extent() # Récupération de la taille de la figure
    a *= bbox.height / bbox.width        # Correction de la pente avec la taille 
    rot = np.degrees(np.arctan(a))       # Calcul de l'angle de rotation
    if cotan:                            # Si on dépasse la verticale
        rot = 90 - np.degrees(np.arctan(1/a))
    t = plt.text(x[N],y[N],label,        # On met le texte au bon endroit
    ha='center',va='center',color=color,rotation = str(rot)) # Avec la bonne rotation
    # On se débrouille pour que la "boîte" d'écriture soit semi-transparente
    #t.set_bbox(dict(facecolor='w',edgecolor='None',alpha=0.8))
    t.set_bbox(dict(boxstyle="round",facecolor='w',edgecolor='None',alpha=0.85))

    def kappa(self, r):
        from numpy import arctanh, arctan, arctan2, log, sin, cos, pi, logspace
        x = self.b / self.rs

        if x < 1.:
            norm = x**2 / (4 * arctanh(
                ((1 - x) / (1 + x))**0.5) / (1 - x**2)**0.5 + 2 * log(x / 2))
        elif x == 1.:
            norm = 1. / (2. + 2 * log(0.5))
        else:
            norm = x**2 / (4 * arctan(
                ((x - 1) / (x + 1))**0.5) / (x**2 - 1)**0.5 + 2 * log(x / 2))

        x = r / self.rs
        A = x * 0.
        C = x < 1.
        X = x[C].copy()
        A[C] = (1. - 2 * arctanh(
            ((1. - X) / (1. + X))**0.5) / (1 - X**2)**0.5) / (X**2 - 1.)
        C = x == 1.
        A[C] = 1. / 3
        C = x > 1.
        X = x[C].copy()
        A[C] = (1. - 2 * arctan(
            ((X - 1.) / (1. + X))**0.5) / (X**2 - 1.)**0.5) / (X**2 - 1.)
        return norm * A