def arc_by_radian(x, y, height, radian_range, thickness, gaussian_width):
    """
    Radial arc with Gaussian fall-off after the solid ring-shaped
    region with the given thickness, with shape specified by the
    (start,end) radian_range.
    """

    # Create a circular ring (copied from the ring function)
    radius = height/2.0
    half_thickness = thickness/2.0

    distance_from_origin = sqrt(x**2+y**2)
    distance_outside_outer_disk = distance_from_origin - radius - half_thickness
    distance_inside_inner_disk = radius - half_thickness - distance_from_origin

    ring = 1.0-bitwise_xor(greater_equal(distance_inside_inner_disk,0.0),greater_equal(distance_outside_outer_disk,0.0))

    sigmasq = gaussian_width*gaussian_width

    if sigmasq==0.0:
        inner_falloff = x*0.0
        outer_falloff = x*0.0
    else:
        with float_error_ignore():
            inner_falloff = exp(divide(-distance_inside_inner_disk*distance_inside_inner_disk, 2.0*sigmasq))
            outer_falloff = exp(divide(-distance_outside_outer_disk*distance_outside_outer_disk, 2.0*sigmasq))
            
    output_ring = maximum(inner_falloff,maximum(outer_falloff,ring))

    # Calculate radians (in 4 phases) and cut according to the set range)

    # RZHACKALERT:
    # Function float_error_ignore() cannot catch the exception when
    # both dividend and divisor are 0.0, and when only divisor is 0.0
    # it returns 'Inf' rather than 0.0. In x, y and
    # distance_from_origin, only one point in distance_from_origin can
    # be 0.0 (circle center) and in this point x and y must be 0.0 as
    # well. So here is a hack to avoid the 'invalid value encountered
    # in divide' error by turning 0.0 to 1e-5 in distance_from_origin.
    distance_from_origin += where(distance_from_origin == 0.0, 1e-5, 0)

    with float_error_ignore():
        sines = divide(y, distance_from_origin)
        cosines = divide(x, distance_from_origin)
        arcsines = arcsin(sines)

    phase_1 = where(logical_and(sines >= 0, cosines >= 0), 2*pi-arcsines, 0)
    phase_2 = where(logical_and(sines >= 0, cosines <  0), pi+arcsines,   0)
    phase_3 = where(logical_and(sines <  0, cosines <  0), pi+arcsines,   0)
    phase_4 = where(logical_and(sines <  0, cosines >= 0), -arcsines,     0)
    arcsines = phase_1 + phase_2 + phase_3 + phase_4

    if radian_range[0] <= radian_range[1]:
        return where(logical_and(arcsines >= radian_range[0], arcsines <= radian_range[1]),
                     output_ring, 0.0)
    else:
        return where(logical_or(arcsines >= radian_range[0], arcsines <= radian_range[1]),
                     output_ring, 0.0)

    def FugaP(self,Z,A,B):
        """ Fugacity Coefficient of Pure Substances"""
##        print Z-B
        L = ( 1/( 2*sqrt(2) ) ) * log( ( Z +B*(1+sqrt(2) ) ) / ( Z +B*(1-sqrt(2) ) ) )
        LogFug = Z-1 - log(Z-B) - A/B*L
        Fug = exp( LogFug )
        return Fug

    def function(self,params):
        """Archemidean spiral function."""

        aspect_ratio = params['aspect_ratio']
        x = self.pattern_x/aspect_ratio
        y = self.pattern_y
        thickness = params['thickness']
        gaussian_width = params['smoothing']
        size = params['size']

        half_thickness = thickness/2.0
        spacing = size*2*pi

        distance_from_origin = sqrt(x**2+y**2)
        distance_from_spiral_middle = fmod(spacing + distance_from_origin - size*arctan2(y,x),spacing)

        distance_from_spiral_middle = minimum(distance_from_spiral_middle,spacing - distance_from_spiral_middle)
        distance_from_spiral = distance_from_spiral_middle - half_thickness

        spiral = 1.0 - greater_equal(distance_from_spiral,0.0)

        sigmasq = gaussian_width*gaussian_width

        with float_error_ignore():
            falloff = exp(divide(-distance_from_spiral*distance_from_spiral, 2.0*sigmasq))

        return maximum(falloff, spiral)

def erf(x):
    """
    Approximation to the erf-function with fractional error
    everywhere less than 1.2e-7

    @param x: value
    @type  x: float

    @return: value
    @rtype: float
    """
    if x > 10.: return 1.
    if x < -10.: return -1.

    z = abs(x)
    t = 1. / (1. + 0.5 * z)

    r = t * N.exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + \
                                                               t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * \
                                                                                                         (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * \
                                                                                                                                               0.17087277)))))))))

    if x >= 0.:
        return 1. - r
    else:
        return r - 1.

def concentricrings(x, y, white_thickness, gaussian_width, spacing):
    """
    Concetric rings with the solid ring-shaped region, then Gaussian fall-off at the edges.
    """

    # To have zero value in middle point this pattern calculates zero-value rings instead of
    # the one-value ones. But to be consistent with the rest of functions the parameters
    # are connected to one-value rings - like half_thickness is now recalculated for zero-value ring:
    half_thickness = ((spacing-white_thickness)/2.0)*greater_equal(spacing-white_thickness,0.0)

    distance_from_origin = sqrt(x**2+y**2)

    distance_from_ring_middle = fmod(distance_from_origin,spacing)
    distance_from_ring_middle = minimum(distance_from_ring_middle,spacing - distance_from_ring_middle)

    distance_from_ring = distance_from_ring_middle - half_thickness

    ring = 0.0 + greater_equal(distance_from_ring,0.0)

    sigmasq = gaussian_width*gaussian_width

    with float_error_ignore():
        falloff = exp(divide(-distance_from_ring*distance_from_ring, 2.0*sigmasq))

    return maximum(falloff,ring)

def sigmoid(axis, slope):
    """
    Sigmoid dividing axis into a positive and negative half,
    with a smoothly sloping transition between them (controlled by the slope).

    At default rotation, axis refers to the vertical (y) axis.
    """
    with float_error_ignore():
        return (2.0 / (1.0 + exp(-2.0*slope*axis))) - 1.0

def smooth_rectangle(x, y, rec_w, rec_h, gaussian_width_x, gaussian_width_y):
    """
    Rectangle with a solid central region, then Gaussian fall-off at the edges.
    """

    gaussian_x_coord = abs(x)-rec_w/2.0
    gaussian_y_coord = abs(y)-rec_h/2.0
        
    box_x=less(gaussian_x_coord,0.0)
    box_y=less(gaussian_y_coord,0.0)
    sigmasq_x=gaussian_width_x*gaussian_width_x
    sigmasq_y=gaussian_width_y*gaussian_width_y

    with float_error_ignore():
        falloff_x=x*0.0 if sigmasq_x==0.0 else \
            exp(divide(-gaussian_x_coord*gaussian_x_coord,2*sigmasq_x))
        falloff_y=y*0.0 if sigmasq_y==0.0 else \
            exp(divide(-gaussian_y_coord*gaussian_y_coord,2*sigmasq_y))

    return minimum(maximum(box_x,falloff_x), maximum(box_y,falloff_y))

def logMean( alpha, beta ):
    """
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: standarddev of log-transformed distribution
    @type  beta: float

    @return: mean of the original lognormal distribution
    @rtype: float
    """
    return N.exp( alpha + (beta**2)/2. )

def logMedian( alpha, beta=None ):
    """
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: not needed
    @type  beta: float

    @return: median of the original lognormal distribution
    @rtype: float
    """
    return N.exp( alpha )

def logSigma( alpha, beta ):
    """
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: standarddev of log-transformed distribution
    @type  beta: float

    @return: 'standard deviation' of the original lognormal distribution
    @rtype: float
    """
    return logMean( alpha, beta ) * N.sqrt( N.exp(beta**2) - 1.)

def exponential(x, y, xscale, yscale):
    """
    Two-dimensional oriented exponential decay pattern.
    """
    if xscale==0.0 or yscale==0.0:
        return x*0.0
    
    with float_error_ignore():
        x_w = divide(x,xscale)
        y_h = divide(y,yscale)
        return exp(-sqrt(x_w*x_w+y_h*y_h))

def hyperbola(x, y, thickness, gaussian_width, axis):
    """
    Two conjugate hyperbolas with Gaussian fall-off which share the same asymptotes.
    abs(x^2/a^2 - y^2/b^2) = 1 
    As a = b = axis, these hyperbolas are rectangular.
    """

    difference = absolute(x**2 - y**2)
    hyperbola = 1.0 - bitwise_xor(greater_equal(axis**2,difference),greater_equal(difference,(axis + thickness)**2))

    distance_inside_hyperbola = sqrt(difference) - axis
    distance_outside_hyperbola = sqrt(difference) - axis - thickness

    sigmasq = gaussian_width*gaussian_width

    with float_error_ignore():
        inner_falloff = exp(divide(-distance_inside_hyperbola*distance_inside_hyperbola, 2.0*sigmasq))
        outer_falloff = exp(divide(-distance_outside_hyperbola*distance_outside_hyperbola, 2.0*sigmasq))
    
    return maximum(hyperbola,maximum(inner_falloff,outer_falloff))

def gabor(x, y, xsigma, ysigma, frequency, phase):
    """
    Gabor pattern (sine grating multiplied by a circular Gaussian).
    """
    if xsigma==0.0 or ysigma==0.0:
        return x*0.0
    
    with float_error_ignore():
        x_w = divide(x,xsigma)
        y_h = divide(y,ysigma)
        p = exp(-0.5*x_w*x_w + -0.5*y_h*y_h)
    return p * 0.5*cos(2*pi*frequency*y + phase)

def gaussian(x, y, xsigma, ysigma):
    """
    Two-dimensional oriented Gaussian pattern (i.e., 2D version of a
    bell curve, like a normal distribution but not necessarily summing
    to 1.0).
    """
    if xsigma==0.0 or ysigma==0.0:
        return x*0.0

    with float_error_ignore():
        x_w = divide(x,xsigma)
        y_h = divide(y,ysigma)
        return exp(-0.5*x_w*x_w + -0.5*y_h*y_h)

    def P(self,T,m):
        P_j= array( exp( m["HAR_A"]+m["HAR_B"] / T  + m["HAR_C"]*log(T) ) )
        P = []
        T1 = log(T)
        T2 = power(T,2)
##        print m["HAR_D"]
        for j in range(len(m["HAR_A"])):
            P_r = P_j[j]
            i= 1
##            print j,P_r,T2
            while i<=20:
                P_i = m["HAR_A"][j]+m["HAR_B"][j]/ T  + m["HAR_C"][j]*T1 + m["HAR_D"][j]*P_r/T2
##                print P_i
                P_i = exp(P_i)
##                print P_i*0.13332236
                i +=1
                if abs(P_i-P_r)<=1:
                    P.append(P_i)
                    break
                P_r = P_i
                
        return array(P) * self.Factor

def ring(x, y, height, thickness, gaussian_width):
    """
    Circular ring (annulus) with Gaussian fall-off after the solid ring-shaped region.
    """
    radius = height/2.0
    half_thickness = thickness/2.0

    distance_from_origin = sqrt(x**2+y**2)
    distance_outside_outer_disk = distance_from_origin - radius - half_thickness
    distance_inside_inner_disk = radius - half_thickness - distance_from_origin

    ring = 1.0-bitwise_xor(greater_equal(distance_inside_inner_disk,0.0),greater_equal(distance_outside_outer_disk,0.0))

    sigmasq = gaussian_width*gaussian_width

    if sigmasq==0.0:
        inner_falloff = x*0.0
        outer_falloff = x*0.0
    else:
        with float_error_ignore():
            inner_falloff = exp(divide(-distance_inside_inner_disk*distance_inside_inner_disk, 2.0*sigmasq))
            outer_falloff = exp(divide(-distance_outside_outer_disk*distance_outside_outer_disk, 2.0*sigmasq))

    return maximum(inner_falloff,maximum(outer_falloff,ring))

def log_gaussian(x, y, x_sigma, y_sigma, mu):
    """
    Two-dimensional oriented Log Gaussian pattern (i.e., 2D version of a
    bell curve with an independent, movable peak). Much like a normal
    distribution, but not necessarily placing the peak above the center,
    and not necessarily summing to 1.0).
    """
    if x_sigma==0.0 or y_sigma==0.0:
        return x * 0.0

    with float_error_ignore():
        x_w = divide(log(x)-mu, x_sigma*x_sigma)
        y_h = divide(log(y)-mu, y_sigma*y_sigma)

        return exp(-0.5*x_w*x_w + -0.5*y_h*y_h)

def line(y, thickness, gaussian_width):
    """
    Infinite-length line with a solid central region, then Gaussian fall-off at the edges.
    """
    distance_from_line = abs(y)
    gaussian_y_coord = distance_from_line - thickness/2.0
    sigmasq = gaussian_width*gaussian_width

    if sigmasq==0.0:
        falloff = x*0.0
    else:
        with float_error_ignore():
            falloff = exp(divide(-gaussian_y_coord*gaussian_y_coord,2*sigmasq))

    return where(gaussian_y_coord<=0, 1.0, falloff)

    def von_mises( self, pars, x ):
        """
        Compute a simplified von Mises function.

        Original formulation in Richard von Mises, "Wahrscheinlichkeitsrechnung
        und ihre Anwendungen in der Statistik und theoretischen Physik", 1931,
        Deuticke, Leipzig; see also Mardia, K.V. and Jupp, P.E., " Directional
        Statistics", 1999, J. Wiley, p.36;
        http://en.wikipedia.org/wiki/Von_Mises_distribution
        The two differences are that this function is a continuous probability
        distribution on a semi-circle, while von Mises is on the full circle,
        and that the normalization factor, which is the inverse of the modified
        Bessel function of first kind and 0 degree in the original, is here a fit parameter.
        """
        a, k, t = pars
        return a * exp( k * ( cos( 2 * ( x - t ) ) - 1 ) )

def disk(x, y, height, gaussian_width):
    """
    Circular disk with Gaussian fall-off after the solid central region.
    """
    disk_radius = height/2.0

    distance_from_origin = sqrt(x**2+y**2)
    distance_outside_disk = distance_from_origin - disk_radius
    sigmasq = gaussian_width*gaussian_width

    if sigmasq==0.0:
        falloff = x*0.0
    else:
        with float_error_ignore():
            falloff = exp(divide(-distance_outside_disk*distance_outside_disk,
                                  2*sigmasq))

    return where(distance_outside_disk<=0,1.0,falloff)