def write_corner(fhandler, r, n, last_id, m_fact=1.0):
    """Write the Ghost particles as well as the common mirror particle, for a
    corner.
    This method is only for convex bodies, and Ghost particles.
    
    Position arguments:
    r -- Mirroring position
    n -- Outward normal of the corner
    last_id -- Number of already written particles
    
    Keyword arguments:
    m_fact -- Mass multiplier factor (To eventually overlap ghost particles)
    
    Returned value:
    Number of written particles
    """
    if BC != 'GP':
        return 0
    n_box = 1
    write_particle(fhandler, r, n=n, mass=dr, imove=-2)
    dir = (math.copysign(1.0, n[0]), math.copysign(1.0, n[1]))
    
    dx = 0.5 * dr
    while dx < sep * h:
        x = r[0] + dir[0] * dx
        dy = 0.5 * dr
        while dy < sep * h:
            y = r[1] + dir[1] * dy
            write_particle(fhandler, (x, y), n=n, mass=m_fact*refd*dr**2,
                           imove=-1, associated=last_id)
            n_box += 1
            dy += dr
        dx += dr
    return n_box

    def assertFloatsAreIdentical(self, x, y):
        """assert that floats x and y are identical, in the sense that:
        (1) both x and y are nans, or
        (2) both x and y are infinities, with the same sign, or
        (3) both x and y are zeros, with the same sign, or
        (4) x and y are both finite and nonzero, and x == y

        """
        msg = 'floats {!r} and {!r} are not identical'

        if isnan(x) or isnan(y):
            if isnan(x) and isnan(y):
                return
        elif x == y:
            if x != 0.0:
                return
            # both zero; check that signs match
            elif copysign(1.0, x) == copysign(1.0, y):
                return
            else:
                msg += ': zeros have different signs'
                if test_support.due_to_ironpython_bug("http://ironpython.codeplex.com/workitem/28352"):
                    print msg.format(x, y)
                    return
        self.fail(msg.format(x, y))

    def apply_matching(self,error_type,matching):
        """ applies appropriate flips to the array to bring it back to the codespace using the given matching """

        channel=0 if error_type=="X" else 1
        
        flips=[]
        
        for pair in matching:
            
            [p0,p1]=pair[0]
            [q0,q1]=pair[1]

            if channel==1 and (p1==-1 or p1==self.size*2+1)and(q1==-1 or q1==self.size*2+1):
                flips+=[]
            elif channel==0 and (p0==-1 or p0==self.size*2+1)and(q0==-1 or q0==self.size*2+1):
                flips+=[]
            else:

                s0=int(math.copysign(1,q0-p0))
                s1=int(math.copysign(1,q1-p1))

                range0=range(1,abs(q0-p0),2)
                range1=range(1,abs(q1-p1),2)
                       
                for x in range1:
                    flips+=[[p0,p1+s1*x]]
                for y in range0:
                    flips+=[[p0+s0*y,q1]]
                
        for flip in flips:
            self.array[flip[0]][flip[1]][channel]*=-1

    def move_vector(self):
        #horizontal?
        if self.x_real == self.x_dest:
            self.x_move = 0
            #up by default
            self.y_move = 1
            #whoops, we need down
            if(self.y_dest < self.y_real):
                self.y_move = -1
                #negative safely unreald slope
                self.m = (-999,1)
            #positive safely unreal slope
            else:
                self.m = (999,1)

        #we have vertical movement, find slope
        else:
            self.m = ((self.y_dest - self.y_real), (self.x_dest - self.x_real))

        self.y_move = self.m[0]/(abs(self.m[0])+abs(self.m[1]))
        self.x_move = 1 - self.y_move

        #make sure the sign is correct
        self.y_move = math.copysign(self.y_move, self.m[0])
        self.x_move = math.copysign(self.x_move, self.m[1])

	def scrollEvent(self, dx=0, dy=0):
		"""
		Generate scroll events from parametters and displacement

		@param int dx		   delta movement from last call on x axis
		@param int dy		   delta movement from last call on y axis
		@param bool free		set to true for free ball move

		@return float		   absolute distance moved this tick

		"""
		# Compute mouse mouvement from interger part of d * scale
		self._scr_dx += dx * self._scr_xscale
		self._scr_dy += dy * self._scr_yscale
		_syn = False
		if int(self._scr_dx):
			self.relEvent(rel=Rels.REL_HWHEEL, val=int(copysign(1, self._scr_dx)))
			self._scr_dx -= int(self._scr_dx)
			_syn = True
		if int(self._scr_dy):
			self.relEvent(rel=Rels.REL_WHEEL,  val=int(copysign(1, self._scr_dy)))
			self._scr_dy -= int(self._scr_dy)
			_syn = True
		if _syn:
			self.synEvent()

 def readline(self):
     if self.readList is None:
         return ''
     n = 0
     timeDiff = self.lastTempAt - time.time()
     self.lastTempAt = time.time()
     if abs(self.temp - self.targetTemp) > 1:
         self.temp += math.copysign(timeDiff *
                                    10, self.targetTemp - self.temp)
         if self.temp < 0:
             self.temp = 0
     if abs(self.bedTemp - self.bedTargetTemp) > 1:
         self.bedTemp += math.copysign(timeDiff *
                                       10, self.bedTargetTemp - self.bedTemp)
         if self.bedTemp < 0:
             self.bedTemp = 0
     while len(self.readList) < 1:
         time.sleep(0.1)
         n += 1
         if n == 20:
             return ''
         if self.readList is None:
             return ''
     time.sleep(0.001)
     return self.readList.pop(0)

def squashMatching(size,error_type,matching):

   channel=0 if error_type=="X" else 1

   flip_array=[[1]*(2*size+1) for _ in range(2*size+1)]

   for [(pt,p0,p1),(qt,q0,q1)] in matching:


      flips = []
    
      if (p0 in [-1,size*2+1] and q0 in [-1,size*2+1]) or (p1 in [-1,size*2+1] and q1 in [-1,size*2+1]):
         flips+=[]

      else:

         s0=int(math.copysign(1,q0-p0))
         s1=int(math.copysign(1,q1-p1))
         
         range0=range(1,abs(q0-p0),2)
         range1=range(1,abs(q1-p1),2)         
         for x in range1:
            flips+=[[p0,p1+s1*x]]
         for y in range0:
            flips+=[[p0+s0*y,q1]]

         
      for flip in flips:
         flip_array[flip[0]][flip[1]]*=-1



   return flip_array

    def update(self, txn):
        if self.sid != txn.sid:
            raise Exception('updating position with txn for a '
                            'different sid')

        total_shares = self.amount + txn.amount

        if total_shares == 0:
            self.cost_basis = 0.0
        else:
            prev_direction = copysign(1, self.amount)
            txn_direction = copysign(1, txn.amount)

            if prev_direction != txn_direction:
                # we're covering a short or closing a position
                if abs(txn.amount) > abs(self.amount):
                    # we've closed the position and gone short
                    # or covered the short position and gone long
                    self.cost_basis = txn.price
            else:
                prev_cost = self.cost_basis * self.amount
                txn_cost = txn.amount * txn.price
                total_cost = prev_cost + txn_cost
                self.cost_basis = total_cost / total_shares

            # Update the last sale price if txn is
            # best data we have so far
            if self.last_sale_date is None or txn.dt > self.last_sale_date:
                self.last_sale_price = txn.price
                self.last_sale_date = txn.dt

        self.amount = total_shares

    def stepThatWay(self):
        """
        Not used.
        """

        fX = self.goal[0] - self.pos.means[0]
        fY = self.goal[1] - self.pos.means[1]

        fM = math.hypot(fX, fY)
        fD = math.atan2(fY, fX)

        # Convert to robot frame
        fDR = fD - self.pos.means[2]
        if math.fabs(fDR) > math.pi:
            fDR -= math.copysign(2 * math.pi, fDR)

        if math.fabs(fDR) > math.pi/10.0:
            provX = 0.0
        else:
            provX = min(self.maxSpeed, self.kSpeed * fM)

        provZ = math.fabs(self.kTurn * fDR)
        provZ = min(self.maxTurn, provZ)
        provZ = math.copysign(provZ, fDR)

        twist = Twist()
        twist.linear.x = provX
        twist.angular.z = provZ

        return twist

    def intersection_point(center_x, center_y, source_x, source_y, \
            width, height=None):
        """Determines where the rhombus centered at (center_x, center_y)
        intersects with a line drawn from (source_x, source_y) to
        (center_x, center_y).

        @see: ShapeDrawer.intersection_point"""
        height = height or width

        if height == 0 and width == 0:
            return center_x, center_y

        delta_x, delta_y = source_x - center_x, source_y - center_y

        # Treat edge case when delta_x = 0
        if delta_x == 0:
            if delta_y == 0:
                return center_x, center_y
            else:
                return center_x, center_y + copysign(height / 2, delta_y)

        width = copysign(width, delta_x)
        height = copysign(height, delta_y)

        f = height / (height + width * delta_y / delta_x)
        return center_x + f * width / 2, center_y + (1-f) * height / 2

    def update_vel(self, delta_time, desired_vel):
        v = float(desired_vel.v)
        omega = float(desired_vel.omega)

        sigma = max(v/self.max_vel, omega/self.max_rot_vel, 1.0)
        if sigma != 1.0:
            if v/self.max_vel > omega/self.max_rot_vel:
                v = copysign(self.max_vel, v)
                omega /= sigma
            else: # omega/self.max_rot_vel > v/self.max_vel
                v /= sigma
                omega = copysign(self.max_rot_vel, omega)

        # these are the velocities wheels would get
        # if they didn't have to accelerate smoothly
        target_lvel = v - 0.5 * self.width * omega
        target_rvel = v + 0.5 * self.width * omega
        if abs(self.lvel - target_lvel) < delta_time * config.BOT_ACCEL_CAP:
            self.lvel = target_lvel
        else:
            self.lvel += copysign(config.BOT_ACCEL_CAP, target_lvel - self.lvel) * delta_time
        if abs(self.rvel - target_rvel) < delta_time * config.BOT_ACCEL_CAP:
            self.rvel = target_rvel
        else:
            self.rvel += copysign(config.BOT_ACCEL_CAP, target_rvel - self.rvel) * delta_time
        # cap velocity
        v = max(abs(self.lvel), abs(self.rvel))
        if v > self.max_vel:
            self.lvel *= self.max_vel / v
            self.rvel *= self.max_vel / v

def drag3(event):
    global x, y, dpx
    dpx=max(abs(start3[0]-event.x)/sizex,abs(start3[1]-event.y)/sizey)
    x=start3[0]+copysign(dpx*sizex, event.x-start3[0])
    y=start3[1]+copysign(dpx*sizey, event.y-start3[1])
    canv.coords(rect,start3[0],start3[1],x,y)
    canv.itemconfig(rect,state='normal')

 def compare(self, flower1, flower2):
     if (flower1.height == flower2.height):
         return 0 # never happen
     if ( self.block(flower1,flower2) ):
         return int(math.copysign(1, flower1.height - flower2.height))
     else:
         return int(math.copysign(1, flower2.height - flower1.height))

    def move_arc(self, radius, theta):
        turn_speed = 0.0
        initial_distance = self.total_distance
        initial_heading = self.heading
        arc_length = radius * theta # some calculation
        while(self.total_distance - initial_distance < abs(arc_length)):
            turn_speed = 0.0
            # calculate desired angle based on distance travelled
            target_heading = (self.total_distance - initial_distance)/radius
            target_heading = math.copysign(target_heading, arc_length)
            # offset by starting angle
            target_heading += initial_heading
            # ensure that 0 <= target_heading <= 2*pi
            if target_heading < 0:
                target_heading += 2*math.pi
            if target_heading > 2*math.pi:
                target_heading -= 2*math.pi

            error = target_heading - self.heading
            if error < -math.pi:
                error += 2*math.pi
            if error > math.pi:
                error -= 2*math.pi
            # correct angle
            if(abs(error) > 0.008):
                if abs(error) < 0.05:
                    turn_speed = math.copysign(0.3, error)  
                else:
                    turn_speed = error*5
            self.set_speeds(0.2, turn_speed)
            print "relative distance: ", self.total_distance - initial_distance
            print "error: ", error
        self.stop_all_motion()
        return 1

    def GetG4XYZ(self):
        """Returns angles (theta,phi,psi) for:
        an x-rotation by theta followed by a y-rotation by phi followed by a z-rotation by psi.
        Passive rotations are considered here, as expected for a GEANT4 rotation. That is equivalent to 
        an Active rotation by -psi,-phi,-theta in the ZYX order.
        Note that this is the same algorithm as GetXYZ(), for the INVERSE matrix.
        This works in all quadrants
        Note that if you construct a rotation matrix using three angles, you may not get back the same angles from this function
        (i.e. the mapping between xyz-angle-triples and SO(3) is not 1-1)"""

        phi= -math.asin(self.item(0,2))

        if math.fabs(1.-math.fabs(self.item(0,2))) >  0.0000001:
            psi  = math.atan2(self.item(0,1),self.item(0,0))
            theta= math.atan2(self.item(1,2),self.item(2,2))
        else:
            psi = 0
            if math.fabs(self.item(1,1)) > 0.0000001 and math.fabs(self.item(2,0))>0.0000001:
                if math.copysign(1,self.item(1,1)) == math.copysign(1,self.item(2,0)):
                    theta=math.atan2(self.item(1,0),self.item(1,1))
                else:
                    theta=math.atan2(-self.item(1,0),self.item(1,1))
            else:
                if not math.copysign(1,self.item(1,0)) == math.copysign(1,self.item(2,1)):
                    theta=math.atan2(self.item(1,0),self.item(1,1))
                else:
                    theta=math.atan2(-self.item(1,0),self.item(1,1))
                
        return(theta,phi,psi)

def VectorMath(ExtPts):
    aDx = ExtPts[1].X-ExtPts[0].X
    aDy = ExtPts[1].Y-ExtPts[0].Y

    bDx = ExtPts[-2].X-ExtPts[0].X
    bDy = ExtPts[-2].Y-ExtPts[0].Y

    VectorA=[aDx,aDy]
    VectorB=[bDx,bDy]
    VectorX=[1,0]
    VectorY=[0,1]

    A_Len = np.linalg.norm(VectorA)
    B_Len = np.linalg.norm(VectorB)

    dpAB=dotProduct(VectorA, VectorB)
    AngAB=math.acos(dpAB/A_Len/B_Len)
    cpAB=(A_Len*B_Len)*math.sin(AngAB)

    dpAX=dotProduct(VectorA, VectorX)
    AngAx=math.acos(dpAX/A_Len/1.0)*math.copysign(1,aDy)

    dpBX=dotProduct(VectorB, VectorX)
    AngBx=math.acos(dpBX/B_Len/1.0)*math.copysign(1,bDy)

    eParts=[A_Len,B_Len, AngAB, AngAx, AngBx,cpAB]

    return(eParts)

    def GetXYZ(self):
        """Returns angles (theta,phi,psi) for:
        an x-rotation by theta followed by a y-rotation by phi followed by a z-rotation by psi.
        Active rotations are considered here.
        This works in all quadrants
        Note that if you construct a rotation matrix using three angles, you may not get back the same angles from this function
        (i.e. the mapping between xyz-angle-triples and SO(3) is not 1-1)"""

        phi=math.asin(-self.item(2,0))
        if math.fabs(1.-math.fabs(self.item(2,0)))>.0000001: 
# if the phi rotation is not by pi/2 or -pi/2 (two solutions, return one of them)
            psi=math.atan2(self.item(1,0),self.item(0,0))
            theta=math.atan2(self.item(2,1),self.item(2,2))
        else: 
# if the phi rotation was in fact by pi/2 or -pi/2 (infinite solutions- pick 0 for psi and solve for theta)
            psi=0
            if math.fabs(self.item(1,1))>.0000001 and math.fabs(self.item(0,2))>.0000001: #as long as these matrix elements aren't 0
                if math.copysign(1,self.item(1,1))==math.copysign(1,self.item(0,2)): 
# relative sign of these matrix elements determines theta's quadrant
                    theta=math.atan2(self.item(0,1),self.item(1,1))
                else:
                    theta=math.atan2(-self.item(0,1),self.item(1,1))
            else:
# if one of them was zero then use these matrix elements to test for theta's quadrant instead
                if not math.copysign(1,self.item(0,1))==math.copysign(1,self.item(1,2)): 
# relative sign of these matrix elements determines theta's quadrant
                    theta=math.atan2(self.item(0,1),self.item(1,1))
                else:
                    theta=math.atan2(-self.item(0,1),self.item(1,1))
                
        return (theta,phi,psi)

    def update(self, dt):
        # call parent update (Animation)
        super().update(dt)

        # Apply gravity if grounded
        if not self.grounded:
            self.velocity_y += self.gravity*dt

        # Set up vars
        c_x = self.get_x(); c_y = self.get_y()
        vx = math.floor(self.velocity_x)
        vy = math.floor(self.velocity_y)

        # Collision handling

        # Collision x
        while not self.can_move(c_x, c_y, c_x + vx, c_y, 0) and not vx == 0:
            vx -= math.copysign(1, vx)

        # Move x
        self.set_position(c_x + vx, c_y)
        c_x = self.get_x()

        # Collision y
        while not self.can_move(c_x, c_y, c_x, c_y + vy, 1) and not vy == 0:
            vy -= math.copysign(1, vy)
            self.velocity_y = 0

        # Check if grounded
        self.grounded = not self.can_move(c_x, c_y, c_x, c_y + 1, 1)

        # Move y
        self.set_position(c_x, c_y + vy)

def almostEqualF(a, b, rel_err=2e-15, abs_err = 5e-323):
    """Determine whether floating-point values a and b are equal to within
    a (small) rounding error.  The default values for rel_err and
    abs_err are chosen to be suitable for platforms where a float is
    represented by an IEEE 754 double.  They allow an error of between
    9 and 19 ulps."""

    # special values testing
    if math.isnan(a):
        return math.isnan(b)
    if math.isinf(a):
        return a == b

    # if both a and b are zero, check whether they have the same sign
    # (in theory there are examples where it would be legitimate for a
    # and b to have opposite signs; in practice these hardly ever
    # occur).
    if not a and not b:
        return math.copysign(1., a) == math.copysign(1., b)

    # if a-b overflows, or b is infinite, return False.  Again, in
    # theory there are examples where a is within a few ulps of the
    # max representable float, and then b could legitimately be
    # infinite.  In practice these examples are rare.
    try:
        absolute_error = abs(b-a)
    except OverflowError:
        return False
    else:
        return absolute_error <= max(abs_err, rel_err * abs(a))

def rowJacobian(x, y, drop_data, tolerances):
    [xP, yP, RP, BP, wP] = drop_data.params
    # x, y = xy # extract the data points
    # s_0 = drop_data.s_0
    if ((x - xP) * cos(wP) - (y - yP) * sin(wP)) < 0:
        s_0 = drop_data.s_left
    else:
        s_0 = drop_data.s_right
    # if x < xP:
    #     s_0 = s_left
    # else:
    #     s_0 = s_right



    xs, ys, dx_dBs, dy_dBs, e_r, e_z, s_i = minimum_arclength(x, y, s_0, drop_data, tolerances) # functions at s* 
    if ((x - xP) * cos(wP) - (y - yP) * sin(wP)) < 0:
        drop_data.s_left = s_i
    else:
        drop_data.s_right = s_i
    # e_i = sqrt(e_r**2 + e_z**2)              # actual residual
    e_i = math.copysign(sqrt(e_r**2 + e_z**2), e_r)              # actual residual
    sgnx = math.copysign(1, ((x - xP) * cos(wP) - (y - yP) * sin(wP))) # calculates the sign for ddi_dX0
    ddi_dxP = -( e_r * sgnx * cos(wP) + e_z * sin(wP) ) / e_i             # derivative w.r.t. X_0 (x at apex)
    ddi_dyP = -(-e_r * sgnx * sin(wP) + e_z * cos(wP) ) / e_i                    # derivative w.r.t. Y_0 (y at apex)
    ddi_dRP = -( e_r * xs + e_z * ys) / e_i  # derivative w.r.t. RP (apex radius)
    ddi_dBP = - RP * (e_r * dx_dBs + e_z * dy_dBs) / e_i   # derivative w.r.t. Bo  (Bond number)
    ddi_dwP = (e_r * sgnx * (- (x - xP) * sin(wP) - (y - yP) * cos(wP)) + e_z * ( (x - xP) * cos(wP) - (y - yP) * sin(wP))) / e_i
    return [[ ddi_dxP, ddi_dyP, ddi_dRP, ddi_dBP, ddi_dwP], e_i]