;+
; NAME:
;	find_launchpoints
; PURPOSE:
;	Find the places (near dipole sources) from whence to launch 
;	field lines.
; CALLING SEQUENCE:
;	launchlist = find_launchpoints(dloc,m,n=n,r=r)
; INPUTS:
;	dloc -- the 3-coordinates ofthe dipoles -- a (3xn) array.
;	m    -- the dipoles' magnetic moment vectors
; KEYWORD INPUTS:
;	n - the number of field lines overall
;	r - if specified, this is the radius from each dipole from which
;		fieldlines should be launched.
; RETURNS:
;	a 3xN array of launchpoints that are close to the dipoles.  
;	In this case, N is only "about equal" to n, because of the 
;	integer nature of fieldlines.
;
; LIMITATIONS:
;	At the moment, only launches fieldlines on the outbound side of the
;	dipoles.  This should be fixed.
; AUTHOR: 
;	Craig DeForest
; HISTORY:
;	Written 16-Dec-97
;-
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