Simple streamline generators

streamlines.py

def fieldlines(self,v1,v2,seeds,dS=1):
    """ Simple Euler fieldline integrator
    v1,v2 - y and x vectors
    seeds - array of coordinate (array indexes) pairs
    dS [optional] - step size
    """
    from scipy import ndimage
    
    field = []
    #Get extent of domain
    max1 = v1.shape[0]
    max2 = v2.shape[1]
    min2 = 0
    min1 = 0
    #integrate for each fieldline
    for seed in seeds:
        c1,c2 = seed
        out1 = [c1] #first point is seed
        out2 = [c2]
        cnt = 0
        
        while (c1 <= max1 and c1 >= min1) and (c2 <= max2 and c2 >= min2):
            #Interpolate the vector field to the coord
            coords = np.array([[c1],[c2]])
            v1i = ndimage.map_coordinates(v1,coords)[0]
            v2i = ndimage.map_coordinates(v2,coords)[0]
            vi = np.sqrt(v1i**2 + v2i**2)
            
            d1 = ( v1i * dS )/ vi
            d2 = ( v2i * dS )/ vi
            c1 -= d1 
            c2 -= d2
            
            out1.append(c1)
            out2.append(c2)
            
            cnt += 1
            if cnt > 500: # Maximum iteration limit
                print "limit"
                break
        out = np.zeros([len(out1),len(out2)])
        out[0] = out1
        out[1] = out2
        field.append(out)
    
    return np.array(field)