R_MAX = 50.0				# limit the radius
PHI_0 = 1.e6				# limit the potential
DS = 0.01				# target step size
D_PHI = 1.e3				# limit the step

def compute_field_line(xi, yi,		# starting point
                       q, xq, yq,	# charges
                       direction):	# +1 = forward step, -1 = backward

    x = [xi]
    y = [yi]
    
    # Loop until the line is too far from the origin
    # or too close to a charge.
    
    while xi*xi+yi*yi < R_MAX*R_MAX \
            and abs(potential(xi, yi, q, xq, yq)) < PHI_0:

        # Compute the field.
        Ex,Ey = field(xi, yi, q, xq, yq)
        E = sqrt(Ex*Ex + Ey*Ey)

        # Choose the step length and limit the change in potential.
        
        ds = DS
        if ds > D_PHI/E: ds = D_PHI/E

        # Take a step in the direction (or opposite) of the field.
        
        xi += direction*ds*Ex/E
        yi += direction*ds*Ey/E

        x.append(xi)
        y.append(yi)

    return x, y