function find_launchpoints,dloc,m,n=n,r=r,verbose=v
common seed,seed
v = (keyword_set(v))
if not isvalid(n) then n = 100

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; np - number of points.
np = nlm(n(0))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; If necessary, calculate a 'small' radius.
if not isvalid(r) then if np eq 1 then r=1 else begin
    if(v) then print,"Calculating radius..."
    pdist = fltarr(3,(n*(n+1))/2)
    idx = 0
    for i=0,np-2 do begin
        pdist(*,idx:idx+np-i-1) = dloc(*,i:np-2) - dloc(*,i)
        idx = idx + np - 1 - i
    end
    pdist = total(pdist*pdist,1)
    r = min(pdist(where(pdist gt 0)))/20.0
    if(v) then print,"r=",r
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Find out how many field lines to draw per unit of flux...
mmag = sqrt(total(m*m,1))
f_n =  n / total(mmag)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Loop over the dipoles and find out where to put 'em for each one...
; The dipole field is concentrated near the pole -- the field strength
; varies around the dipole.  So we have to stochastically launch them 
; with the right distribution.  The right distribution in colatitude 
; about the dipole is sqrt((1 + 3 (cos(theta))^2))*sin(theta).
; That function has a maximum of about 1.155.  So we divide it by 1.2
; and then use it as a rejection criterion for randomly chosen latitudes.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
out = dblarr(4,2*n)
out(0)=1
out(0,1) = -1

if(f_n * max(mmag) lt 0.5) then message,"Hey -- not enough lines for those dipoles!"

for pno=0,np-1 do begin
    for i=0,fix(mmag(pno)*f_n + 0.5) do begin
        ;;; Construct a basis for the little coordinate system
        ;;; around the point
        x = normal(cross([1,0,0],m(*,pno)))
            if(max(x*x) eq 0) then x = normal(cross([0,1,0],m(*,pno)))
        y = normal(cross(m(*,pno),x))
        z = normal(m(*,pno))
        
        ;;;
        ;;;Pick a longitude
        repeat begin
            phi = randomu(seed)*!pi/2
        end until randomu(seed) lt (sqrt((1 + 3 *(cos(phi))^2))*sin(phi))/1.2
        theta = randomu(seed)*!pi*2
        
        dx = cos(phi)*cos(theta)
        dy = cos(phi)*sin(theta)
        dz = sin(phi)
        
        offset = ([[x],[y],[z]]#[dx,dy,dz])*r   

        out(*,(out(0)))=[pno,dloc(*,pno)+offset]
        out(0)=out(0)+1
    end
end

return,out
end