;+ ; NAME: ; trace_field ; PURPOSE: ; Trace field lines in a flow field (such as the B-field) ; CALLING SEQUENCE: ; trace_field,field,location,ptlist ; INPUTS: ; FIELD - an N+1-dimensional array. The first Ndimensions are the ; independent variables of the flow field; the last ; dimension should have ordinality N; each of the N N-cubes ; represents one dimensionof theflowfield. ; location - An N-vector that indicates where to start ; OUTPUTS: ; ptlist- a list of the visited points (see "every" keyword). ; An NxM array, where N is the dimension of the flow field and ; M is the number ofpoints plotted. ; KEYWORD INPUTS: ; /forward - trace forward (default unless /backward is set)) ; /backward - trace backward ; /both - shorthand for "/forward,/backward" ; every - set to the number of pixels to skip between outputs ; (as in "every 4 pixels") ; separation - if set, we stop if we come within "separation" ; voxels of a previously visited one. ; maxticks - how many steps we're allowed to take before ; quitting. Defaults to (10/step) * max dimension ; KEYWORDS FOR INPUT/OUTPUT ; out - An integer array the size of the input field, indicating ; which voxels have been visited. ; ptno - the number at which to start marking visited voxels. ; ; RESTRICTIONS: ; Only works for N=2 and N=3. ; ; AUTHOR: Craig DeForest ; HISTORY: Written 10-Dec-1997 ;- pro trace_field,field,location,ptlist,separation=separation,every=every,out=out,ptno=ptno,maxticks=maxticks,plot=plot,verbose=verbose,debug=debug,ranges=ranges,small=small,ticktime=ticktime,steprecord=steprecord,forward=forward,backward=backward,both=both debug=keyword_set(debug) if isvalid(verbose) then v = verbose else v = 0 if not keyword_set(backward) and not keyword_set(forward) then forward=1 if keyword_set(both) then begin forward=1 backward=1 end if(not isvalid(location)) then message, "Must specify a starting location" n = nlm(location) initial = call_function(field,0,location) if (nlm(initial) ne n) then message,"Function "+field+"("+fug(location)+") didn't return the proper number of dimensions (wanted "+fug(n)+"; got "+fug(nlm(initial))+")" if(not isvalid(ranges)) then message,"Must specify a 2xN array of ranges" rs = size(ranges) if((rs(0) ne 2) or (rs(1) ne 2) or (rs(2) ne n)) then $ message,"Ranges must be a 2xN array!" if not isvalid(separation) then separation=1 if not isvalid(steprecord) then steprecord=separation if not isvalid(ticktime) then ticktime = 100 if not isvalid(ptno) then ptno = 1 if not isvalid(maxticks) then maxticks = 5000 if not isvalid(location) then begin location = double(fs(1))/2 for i=1,n-1 do location=[location,double(fs(i+1))/2] end plot = keyword_set(plot) and (n eq 2) if not isvalid(small) then small = 1d-2 else small=double(small) if not isvalid(every) then every = 1 fw = ( keyword_set(both) or $ keyword_set(forward) or $ not keyword_set(backward) ) bw = ( keyword_set(both) or $ keyword_set(backward) ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;; Done with parameter setup -- configure ;;;;;;;;;; tracing internals ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; for bk = (bw eq 0),(fw ne 0) do begin ticks = 0 t0 = 0 lastpt = ptno if(bk eq 1) then h = -1d5 else h = 1d5 print,"bk = ",bk potential = 0d0 l0 = double(location) tooclose=0 ;; Mark and output our first point if not isvalid(ptlist) then begin ptlist = [[2d0,l0*0d0],[-1,l0*0d0],[ptno,l0]] end else begin ptlist = [[ptlist(*,0:ptlist(0))],[-1,l0*0d0],[ptno,l0]] ptlist(0) = ptlist(0)+2 end p0 = l0 l3 = l0 ptstart = ptlist(0) ptno=ptno+1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;; Done with setup -- start tracing. ;;;;;;;;;; We use the RK4 algorithm with variable stepsize, using ;;;;;;;;;; half-steps to determine when we need to adjust the ;;;;;;;;;; stepsize. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print,"Starting trace..." repeat begin l0 = l3 repeat begin f0 = call_function(field,0,l0) l1 = rk4(l0,f0,potential,h,field,/double) l2 = rk4(l0,f0,potential,h/2,field,/double) l3 = rk4(l2,call_function(field,0,l2),potential+h/2,h/2,field,/double) diff = abs((l3 - l1))/(abs(l3 - l0)>small) notsogood = diff gt small if(total(notsogood) ne 0) then begin h = h / 2.0 end end until total(notsogood) eq 0 potential = potential + h if( total(diff gt (small/10.0)) eq 0 ) then begin h = h * 1.1 end ;; Check bounds before proceeding outofbounds = (l3 lt ranges(0,*)) or (l3 gt ranges(1,*)) p0dist = sqrt(total((l3 - p0) * (l3 - p0))) ticks=ticks+1 if(v gt 1) then begin print,"tick #",fug(ticks),"; h=",fug(h),"; Location=",l3 print," diff=",diff print," l1 =",l1 print," l3 =",l3 end if(p0dist ge steprecord) then begin if(v eq 1) then print,"*** tick #",fug(ticks),"; h=",fug(h),"; Location=",l3 if(v gt 1) then print,"-----------step-----------" t0 = ticks p0 = l0 ptno = ptno + 1 if(ptno-lastpt ge every) then begin ptlist(0) = ptlist(0)+1 if(ptlist(0) ge nlm(ptlist(0,*))) then begin ptlist = [[ptlist],[ptlist]] ptlist( *, ptlist(0) : (nlm(ptlist(0,*))-1) )=0d0 end ptlist(*,ptlist(0)) = [ptno,l3] if(plot and ptlist(0)-ptstart gt 1) then plots,[ptlist(1,ptlist(0)-2:ptlist(0)-1)-ranges(0,0)],[ptlist(2,ptlist(0)-2:ptlist(0)-1)-ranges(0,1)],color=255,/device end if(ptlist(0)-ptstart gt 3) then begin d = (ptlist(1:n,ptstart:ptlist(0)-2)) for i=0,n-1 do d(i,*) = d(i,*)-l3(i) d=total(d*d,1) tooclose = (total(d lt separation*separation) gt 0) end end d2 = total((l3-l0)*(l3-l0)) end until (total(outofbounds) gt 0) or (ticks gt maxticks) or (tooclose) or (d2 lt small*small) if(total(outofbounds) gt 0) then print,"outofbounds = ",outofbounds if(ticks gt maxticks) then print,"ticks exceeded",maxticks if(tooclose) then print,"Too close to itself" if(d2 lt small) then print,"Got near a singularity" end end