This function calculates the intersection points for two circles defined by (x1, y1, r1) and (x2, y2, r2). It returns the coordinates of the two intersection points, which can also be used to calculate the intersection line for the circles.

Intersection Line for two circles

def intersection_line(x1, y1, r1, x2, y2 r2):
    """Calculates intersection line for two circles defined by (x1, y1, r1) and (x2, y2, r2)."""
    # Formulas obtained from
    # http://www.wolframalpha.com/input/?i=solve+x^2+%2B+y^2+%3D+pow(r,2)+
    # and+%28x-a%29^2+%2B+%28y-b%29^2+%3D+pow(q,2)+for+x%2Cy

    # Create transformed coordinates for calculations
    a = x2 - x1
    b = y2 - y1
    (r,q) = (r1, r2)

    # Check if circles are intersecting, but not contained within each other
    distance = sqrt(a*a + b*b)
    if (distance + r1) <= r2 or (distance + r2) <= r1:
        print 'One circle contained in another, intersection points not defined'
        return None
    if distance > (r1 + r2):
        print 'Circles do not intersect, intersection points not defined'
        return None
    
    if b == 0:
        # Calculate solutions for simple case
        sol_x1 = (pow(a,2)-pow(q,2)+pow(r,2))/(2*a)
        sol_x2 = sol_x1

        sol_y1 = sqrt(pow(r,2)-pow(pow(a,2)-pow(q,2)+pow(r,2),2)/(4*pow(a,2)))
        sol_y2 = -sol_y1
    else:
        # Calculate solutions for x coordinate
        const1 = pow(a,3)
        const2 = sqrt(-pow(b,2)*(pow(a,4)+2*pow(a,2)*pow(b,2)-2*pow(a,2)*pow(q,2)\
            -2*pow(a,2)*pow(r,2)+pow(b,4)-2*pow(b,2)*pow(q,2)-2*pow(b,2)*pow(r,2)\
            +pow(q,4)-2*pow(q,2)*pow(r,2)+pow(r,4)))
        const3 = a*pow(b,2)-a*pow(q,2)+a*pow(r,2)
        const4 = 2*(pow(a,2)+pow(b,2))
        
        sol_x1 = (const1 - const2 + const3)/const4
        sol_x2 = (const1 + const2 + const3)/const4

        # Calculate solutions for y coordinate
        const5 = pow(a,2)*pow(b,2)
        const6 = a*sqrt(-pow(b,2)*(pow(a,4)+2*pow(a,2)*pow(b,2)-2*pow(a,2)\
            *pow(q,2)-2*pow(a,2)*pow(r,2)+pow(b,4)-2*pow(b,2)*pow(q,2)-2*pow(b,2)\
            *pow(r,2)+pow(q,4)-2*pow(q,2)*pow(r,2)+pow(r,4)))
        const7 = pow(b,4)-pow(b,2)*pow(q,2)+pow(b,2)*pow(r,2)
        const8 = 2*b*(pow(a,2)+pow(b,2))

        sol_y1 = (const5 + const6 + const7)/const8
        sol_y2 = (const5 - const6 + const7)/const8

    # Create coordinate tuples for the two solutions
    solution1 = (sol_x1 + x1, sol_y1 + y1)
    solution2 = (sol_x2 + x1, sol_y2 + y1)

    return [solution1, solution2]