;+ ; ; NAME: ilct.pro ; ; PURPOSE: Finds velocity field V_REC (in pixels/time step) ; that satisfies three equations --- ; 1) curl(U_REC B_z) = curl(U_LCT B_z), ; where U_* := V_perp - B_perp*(V_z/B_z) ; 2) dB_z/dt = div(V cross B)_perp ; 3) V dot B = 0 ; Method finds phi(x,y) and psi(x,y) such that ; (V_z B_perp - V_perp B_z) = - grad_perp phi + curl(psi z_hat) ; ; ; REQUIRED INPUTS: ; mag = (nx, ny, 4) fltarr of B-field values [Gauss], with ; time-centered field (B_x,B_y,B_z) in(*,*,0:2), ; and change in B_z (final - initial) in (*,*,3) ; vel = (nx, ny, 2) fltarr, with (V_x,V_y) [cm/s] from LCT in (*,*,0:1) ; dx = pixel size, in cm ; dt = time step, in sec ; ; KEWORDS: ; thr = threshold in |B_z| below which velocities are zeroed; ; default is 100 G; if (thr < 1), assumed fraction of max(|B_z|) ; mask = set by default; zeros out results in |B_z| < thr regions ; check = set to check that induction eqn. is satisfied. ; generates (nx,ny,3) fltarr "induct_check:" ; induct_check(*,*,0)= dB_z/dt from input ; induct_check(*,*,1)= -div(U_LCT B_z); LCT's expected dB_z/dt ; induct_check(*,*,2)= -div(U_ILCT B_z); ILCT's expected dB_z/dt ; ; OUTPUT: ; vel = (nx, ny, 3) array of (V_x,V_y,V_z), in cm/s ; phi = (optional) "inductive" potential scalar field found by ILCT ; psi = (optional) "electrostatic" potential scalar field found by ILCT ; ; HISTORY: Written by BTW, Jun 2003. ; Modified to handle regions where B_z = 0, BTW, Oct 2004. ; Modified again for "v.small" B_z values, BTW, Nov 2004. ; ; USAGE: ilct, vel, mag, dx, dt, thr=thr, check=1, induct_check=induct_check ; ;- ; Three subroutines follow; will compile with main procedure. ;============================================================ ;============================================================ function poisson, infield, dx, inspect=inspect, outspect=outspect ; This subroutine solves 2D Poission's equation: ; the Laplacian of the 2D scalar output ; field gives the 2D scalar input field ; Approach: for each possible value of (k_x, k_y), find output ; FFT spectrum by dividing input FFT spectrum by (k_x^2 + k_y^2). dims = size(infield, /dim) nx = dims(0) ; define dimensionality ny = dims(1) ; define arrays outfield = fltarr(nx,ny) ; this will be output array outspect = complexarr(nx,ny) ; this spectrum w/FFT gives outfield ; define frequency arrays as IDL does -- for even and odd cases if ((nx)mod(2) eq 0) then xfreqs = $ [0,(findgen(nx/2) + 1)/(nx),-reverse(findgen(nx/2-1) + 1)/(nx)] $ else xfreqs = $ [0,(findgen(nx/2.-.5) + 1)/(nx),-reverse(findgen(nx/2.-.5) + 1)/(nx)] if ((ny)mod(2) eq 0) then yfreqs = $ [0,(findgen(ny/2) + 1)/(ny),-reverse(findgen(ny/2-1) + 1)/(ny)] $ else yfreqs = $ [0,(findgen(ny/2.-.5) + 1)/(ny),-reverse(findgen(ny/2.-.5) + 1)/(ny)] inspect = fft(infield, -1) for i = 0,nx-1 do begin for j = 0,ny-1 do begin ; laplacian of 0-modes is zero, so start w/non-zero modes. if ((i ne 0) or (j ne 0)) then $ outspect(i,j) = -(2*!pi)^(-2)/ $ ((xfreqs(i))^2+(yfreqs(j))^2) $ *inspect(i,j) endfor endfor outfield = real_part(fft(outspect,1)) if not(keyword_set(dx)) then dx = 1. outfield = outfield*dx^2 return, outfield end ;============================================================ ;============================================================ function xderiv, infield, dx, inspect=inspect, outspect=outspect ; This subroutine uses FFT's to compute the ; the 2D x-derivative of the 2-d input array. ; Approach: take FFT, multiply by ikx, transform back if not(keyword_set(dx)) then dx = 1. dims = size(infield, /dim) nx = dims(0) ; define dimensionality ny = dims(1) ; define arrays outfield = fltarr(nx,ny) ; this will be output array outspect = complexarr(nx,ny) ; this spectrum w/FFT gives outfield ; define frequency arrays as IDL does -- for even and odd cases if ((nx)mod(2) eq 0) then xfreqs = $ [0,(findgen(nx/2) + 1)/(nx),-reverse(findgen(nx/2-1) + 1)/(nx)] $ else xfreqs = $ [0,(findgen(nx/2.-.5) + 1)/(nx),-reverse(findgen(nx/2.-.5) + 1)/(nx)] inspect = fft(infield, -1) for i = 0,nx-1 do begin for j = 0,ny-1 do begin ; derivative of 0-mode is zero, so start w/non-zero modes. if ((i ne 0) or (j ne 0)) then $ outspect(i,j) = (2*!pi)* $ ( xfreqs(i) *complex(0.,1.)) $ *inspect(i,j) endfor endfor outfield = real_part(fft(outspect,1)) outfield = outfield/dx return, outfield end ;============================================================ ;============================================================ function yderiv, infield, dy, inspect=inspect, outspect=outspect ; This subroutine uses FFT's to compute the ; the 2D y-derivative of the 2D input array. ; Approach: take FFT, multiply by iky, transform back. if not(keyword_set(dy)) then dy = 1. dims = size(infield, /dim) nx = dims(0) ; define dimensionality ny = dims(1) ; define arrays outfield = fltarr(nx,ny) ; this will be output array outspect = complexarr(nx,ny) ; this spectrum w/FFT gives outfield ; define frequency arrays as IDL does -- for even and odd cases if ((ny)mod(2) eq 0) then yfreqs = $ [0,(findgen(ny/2) + 1)/(ny),-reverse(findgen(ny/2-1) + 1)/(ny)] $ else yfreqs = $ [0,(findgen(ny/2.-.5) + 1)/(ny),-reverse(findgen(ny/2.-.5) + 1)/(ny)] inspect = fft(infield, -1) for i = 0,nx-1 do begin for j = 0,ny-1 do begin ; derivative of 0-mode is zero, so start w/non-zero modes. if ((i ne 0) or (j ne 0)) then $ outspect(i,j) = (2*!pi)* $ ( yfreqs(j) *complex(0.,1.)) $ *inspect(i,j) endfor endfor outfield = real_part(fft(outspect,1)) outfield = outfield/dy return, outfield end ;============================================================ ;============================================================ ; MAIN PROCEDURE ;============================================================ ;============================================================ pro ilct, vel, mag, dx, dt, thr=thr, phi=phi, psi=psi, $ check=check, induct_check=induct_check, mask=mask, $ verbose=verbose nparms=N_Params() if (nparms lt 4) then begin Message, 'Usage: ilct, vel, mag, dx, dt', /info Message, 'Purpose: Compute velocities from LCT and normal', /info Message, ' component of induction equation.', /info Message, 'Input: vel = (*,*,2) fltarr of LCT results, cm/s', /info Message, ' mag = (*,*,4) fltarr of B-field comps', /info Message, ' in (*,*,0:2) and change in B in (*,*,3)', /info Message, ' dx = pixel size, in cm', /info Message, ' dt = time interval, in seconds',/info Message, 'Output: vel = (*,*,3) fltarr of velocities, cm/s', /info Message, 'Require: idl version 5.4 or more recent', /info Message, 'Written: B.Welsch, June 2003', /info return endif if (n_elements(mask) eq 0) then mask = 0 if (n_elements(hanning) eq 0) then hanning = 0 nx = n_elements(mag(*,0,0)) ny = n_elements(mag(0,*,0)) u = vel(*,*,0) ; Make inputs easier to reference! v = vel(*,*,1) Bx = mag(*,*,0) ; Time-centered B field components. By = mag(*,*,1) Bz = mag(*,*,2) dBzdt = mag(*,*,3)/float(dt) ; Change in B_z over dt, final - initial ; FIND PHI ;========== fftphi = poisson(dBzdt,dx) ;jacobi, dBzdt, phi = phi, tol=tol, dx=dx, /verbose ;jacphi = phi phi = fftphi dphidx = xderiv(phi, dx) ; and grad phi dphidy = yderiv(phi, dx) ddphiddx = xderiv(dphidx, dx) ; and Laplacian phi ddphiddy = yderiv(dphidy, dx) lapphi = ddphiddx + ddphiddy cm = check_math() if (cm ne 0) and (cm ne 32) then print,temporary(fail_here) ; FIND PSI ;========== uxBz = u*Bz uyBz = v*Bz curlulBz = (xderiv(uyBz,dx) - yderiv(uxBz,dx)) fftpsi = poisson(-curlulBz, dx) ;jacobi, -curlulBz, phi = psi, tol=tol, dx=dx, /verbose ;jacpsi = psi psi = fftpsi dpsidx = xderiv(psi, dx) ; and grad psi dpsidy = yderiv(psi, dx) ddpsiddx = xderiv(dpsidx, dx) ; and Laplacian psi ddpsiddy = yderiv(dpsidy, dx) lappsi = ddpsiddx + ddpsiddy cm = check_math() if (cm ne 0) and (cm ne 32) then print,temporary(fail_here) ; WHERE |B_z| > nzthr defines non-zero data ;=========================================== if not(keyword_set(thr)) then thr = 100 if (thr le 1) then nzthr = thr*max(abs(Bz)) else nzthr = thr hiBz = where( abs(Bz) gt nzthr, n_hiBz) loBz = where( abs(Bz) le nzthr, n_loBz) zeros = where(abs(Bz) lt 1e-4, n_zeros) if (n_zeros ne 0) then begin Bzmiss = 1e+10*max(abs(Bz)) Bz(zeros) = Bzmiss print,'% ILCT: Zeros present, adjusting.' endif else print,'% ILCT: No zeros present.' vx = fltarr(nx,ny) vy = fltarr(nx,ny) vz = fltarr(nx,ny) Fx = -dphidx + dpsidy Fy = -dphidy - dpsidx magb2 = (Bx^2 + By^2 + Bz^2) if (n_hiBz eq 0) then begin print,"% ILCT: No magnetic data above threshold! Stopping..." print, fail_variable endif else begin vz = -(Bx*Fx + By*Fy)/magb2 vx = ( Fx + Bx*vz)/Bz vy = ( Fy + By*vz)/Bz endelse ; IF KEYWORD_SET(check) THEN BEGIN CHECKING... ;====================================================== if keyword_set(check) then begin induct_check = fltarr(nx,ny,4) ; Construct U_REC B_z = V_perp * B_z - B_perp * V_z xsoln = (vx*Bz - vz*Bx) ysoln = (vy*Bz - vz*By) ; Div. of this solution vector should = - dB_z/dt divsoln = (xderiv(xsoln,dx) + $ yderiv(ysoln,dx)) ; Div. of scalar fields' derivatives should = - dB_z/dt divsclrs = (xderiv(Fx,dx) + $ yderiv(Fy,dx)) ; How does LCT alone do? Find div (V_LCT * B_z). smp = 5 ; helps if smoothing's applied divulctBz = (xderiv(smooth(Bz*u,smp), dx) + $ yderiv(smooth(Bz*v,smp), dx)) if keyword_set(mask) then begin dBzdt(loBz) = 0. divulctBz(loBz) = 0. divsoln(loBz) = 0. divsclrs(loBz) = 0. endif ; Generate output array induct_check(*,*,0) = dBzdt induct_check(*,*,1) = -divulctBz induct_check(*,*,2) = -divsoln induct_check(*,*,3) = -divsclrs endif if keyword_set(mask) then begin vx(loBz) = 0. vy(loBz) = 0. vz(loBz) = 0. endif vel = [[[vx]],[[vy]],[[vz]]] end