Download
;+
;NAME:
; bff
;PURPOSE:
; Computes the linear force-free field (constant alpha)
;CATEGORY:
;CALLING SEQUENCE:
; magout = bff(magin)
;INPUTS:
; magin = magnetic field structure (from bfits.pro)
;OPTIONAL INPUT PARAMETERS:
;KEYWORD PARAMETERS
; alpha = The value of alpha (0 is the default (potential)).
; Be sure to get the sign right!
; nz = Number of equally spaced grid points in the z direction
; (default = 30)
; zscale = Sets the z-scale (1.0 = same scale as the x,y axes before the
; heliographic transformation). Default= 1.0.
; 0.1 is probably the minimum you would want to use.
; /quiet = If set, don't display all the plots
; binup = Decreases grid spacing by given (integer) factor by
; interpolation of the original magnetogram.
; filter = If present, applies digital filter to the FFT (not too useful)
; /no_heliographic = if set, put the output back to the image coordinates.
; Still returns Bx, By, and BZ, though.
; /rzdep = when /no_heliographic is set, normally only z=0 is returned.
; If /rzdep is set the full z dependence is returned.
; /image = put the output back to the image coordinates and return the
; longitudinal and transverse field as Bx, By, Bz
; /doxform = force bff to do heliographic transformation
; /finite_energy = if set, restricts the soutions to those with finite
; energy by turning off the oscillatory Fourier
; components.
; surround = Surround the magnetogram with zeros to make the FFTs less
; sensitive to replication. This will slow down the program
; quite a bit, but, in principle, makes the result more
; accurate. Padded with zeroes such that
; nx,ny -> (surround+1)*(nx,ny): /surround doubles the sizes.
; The zeroes are stripped off on output; the padded zeroes
; are only used internally. xysize below is similar but the
; zeroes are retained on output.
; xyexpand = A 2-element array which gives the x-y expansion of the input.
; xysize works like surround, except that the padded zeroes are
; left in the output. This is useful for plotting field lines
; which leave the input box. default = [1.,1.]. If both
; xyexpand and surround are defined, xyexpand gives the
; expansion factors, but the zeroes are stripped as in
; surround.
; /get_alpha = Compute the best value of alpha (returned in alpha keyword)
;OUTPUTS:
; magout = new magnetic field structure with the fields bxp,byp,bzp defined
; These are 3D arrays in x,y,z giving the x,y,z components of the
; force free field.
;COMMON BLOCKS:
;SIDE EFFECTS:
;RESTRICTIONS:
; The magnetic field input must already have the 180 degree ambiguity
; resolved.
;PROCEDURE:
; Uses the FFT equations from Gary, ApJ, 69, 323, 1989. By default the
; solution is NOT restricted to finite energy. This can be overcome with
; the /finite_energy keyword.
;MODIFICATION HISTORY:
;
; Written by T. Metcalf 1/31/92
; 30-Jul-93 TRM Fixed a bug in the computation of the FFT
; frequencies
; 10-Nov-03 TRM Fixed bug in calculation of cff.
;-
function bff,magin,nz=nz,alpha=alpha,quiet=quiet,zscale=zscale,binup=binup, $
filter=filter,no_heliographic=no_heliographic,doxform=doxform, $
finite_energy=finite_energy,surround=surround, $
get_alpha=get_alpha,xyexpand=xysize,image=image,rzdep=rzdep
if n_elements(xysize) NE 2 then begin ; Sets the input expansion
xysize = [1.,1.]
if n_elements(surround) EQ 1 then xysize = xysize*((surround+1)>1)
endif
xysize = xysize>[1.,1.]
if n_elements(zscale) LE 0 then begin
if tag_index(magin,'zscale',/quiet) GE 0 then zscale=magin.zscale $
else zscale = 1.0
endif
if tag_index(magin,'zscale',/quiet) GE 0 then begin
if zscale NE magin.zscale then print,'WARNING: bff: Inconsistent zscale.'
endif
if n_elements(alpha) EQ 0 then alpha = 0.0
if tag_index(magin,'POINT') LT 0 then begin
print,'ERROR: bff: Structure field "POINT" is missing from mag structure'
return,0
endif
; Do the transformation to heliographic coordinates if necessary
xf = set_xform(magin.point.b0,magin.point.p,magin.point.lat,magin.point.cmd)
if tag_index(magin,'BX') LT 0 OR $
tag_index(magin,'BY') LT 0 OR $
tag_index(magin,'BZ') LT 0 OR $
tag_index(magin,'B_PXD') LT 0 OR $
tag_index(magin,'B_PYD') LT 0 OR keyword_set(doxform) then begin
if tag_index(magin,'B_LONG') GE 0 AND $
tag_index(magin,'B_TRANS') GE 0 AND $
tag_index(magin,'B_AZIM') GE 0 AND $
tag_index(magin,'POINT') GE 0 then begin
print,'INFO: bff: Computing heliographic transformation'
nx = n_elements(magin.b_long(*,0))
ny = n_elements(magin.b_long(0,*))
pix_size=double(magin.point.pix_size)
if n_elements(pix_size) EQ 1 then pix_size = [pix_size(0),pix_size(0)]
mpix_size = min(pix_size)
x = (lindgen(nx,ny) MOD nx)*pix_size(0)/mpix_size
y = (long(lindgen(nx,ny) / nx))*pix_size(1)/mpix_size
; xform assumes x is West
xform, -magin.b_trans*sin(magin.b_azim*!dtor), $
magin.b_trans*cos(magin.b_azim*!dtor), $
magin.b_long, $
bdx,bdy,bdz, $
transform=xf.image2helio
xform,x,y,pxd,pyd,transform=xf.zimage2helio ; Gets helio pixel coords
if tag_index(magin,'BX',/quiet) LT 0 then $
magin=add2str(magin,['bx'],bdx,/quiet) $
else magin.bx=bdx
if tag_index(magin,'BY',/quiet) LT 0 then $
magin=add2str(magin,['by'],bdy,/quiet) $
else magin.by=bdy
if tag_index(magin,'BZ',/quiet) LT 0 then $
magin=add2str(magin,['bz'],bdz,/quiet) $
else magin.bz=bdz
if tag_index(magin,'B_PXD',/quiet) LT 0 then $
magin=add2str(magin,['b_pxd'],pxd,/quiet) $
else magin.b_pxd=pxd
if tag_index(magin,'B_PYD',/quiet) LT 0 then $
magin=add2str(magin,['b_pyd'],pyd,/quiet) $
else magin.b_pyd=pyd
endif $
else begin
print,'ERROR: bff: Not enough information in the input structure'
return,0
endelse
endif
; Decrease grid spacing by interpolation. Also protects the input structure
; from change
if n_elements(binup) GT 0 then begin
mag = binup_b(magin,binup,alpha)
nx = n_elements(magin.bz(*,0))
ny = n_elements(magin.bz(0,*))
nnx = fix(nx * binup)
nny = fix(ny * binup)
binupx = float(nnx)/float(nx)
binupy = float(nny)/float(ny)
endif $
else begin
mag = magin
binup = 1.
binupx = 1.
binupy = 1.
endelse
; Put the magnetogram on a uniform grid in heliographic coordinates
; Uses the Bx, By, Bz variables.
nx = n_elements(mag.bz(*,0))
ny = n_elements(mag.bz(0,*))
x = (lindgen(nx,ny) MOD nx)
y = (long(lindgen(nx,ny) / nx))
; p,q will be the forward poly_2d transformation and pinv,qinv will be for
; the inverse transform which puts the grid back to the original at the
; end of the computations.
bpxd = mag.b_pxd*binupx ; Binup is 1 unless grid spacing has been changed
bpyd = mag.b_pyd*binupy
if tag_index(mag,'J_PXD',/quiet) GE 0 and $
tag_index(mag,'J_PYD',/quiet) GE 0 then begin
jpxd = mag.j_pxd*binupx ; Binup is 1 unless grid spacing has been changed
jpyd = mag.j_pyd*binupy
endif
polywarp,x,y,bpxd-min(bpxd),bpyd-min(bpyd),3,p,q
polywarp,bpxd-min(bpxd),bpyd-min(bpyd),x,y,3,pinv,qinv
if n_elements(jpxd) GT 0 and n_elements(jpyd) GT 0 then $
polywarp,x,y,jpxd-min(jpxd),jpyd-min(jpyd),3,pj,qj
; x2=fltarr(4) & y2=fltarr(4)
; temp = xf.zimage2helio # [0,0] & x2(0) = temp(0) & y2(0)=temp(1)
; temp = xf.zimage2helio # [nx-1,0] & x2(1) = temp(0) & y2(1)=temp(1)
; temp = xf.zimage2helio # [0,ny-1] & x2(2) = temp(0) & y2(2)=temp(1)
; temp = xf.zimage2helio # [nx-1,ny-1] & x2(3) = temp(0) & y2(3)=temp(1)
; xfzi = invert(xf.zimage2helio)
; pinv = [-min(x2),xf.zimage2helio(0,1),xf.zimage2helio(0,0),0]
; qinv = [-min(y2),xf.zimage2helio(1,1),xf.zimage2helio(1,0),0]
; p = [min(x2),xfzi(0,1),xfzi(0,0),0]
; q = [min(y2),xfzi(1,1),xfzi(1,0),0]
nx = (fix(max(bpxd) - min(bpxd)) + 1) ; Make sure that pixel scale
ny = (fix(max(bpyd) - min(bpyd)) + 1) ; does not change.
bx2 = poly_2d(mag.bx,p,q,1,nx,ny,missing=0.0) ; Field on heliographic grid
by2 = poly_2d(mag.by,p,q,1,nx,ny,missing=0.0)
bz2 = poly_2d(mag.bz,p,q,1,nx,ny,missing=0.0)
; Print some information about the transformation
;pq2rss,p,q,erot,exscl,eyscl,exshft,eyshft,enrss,quiet=quiet
nxs = long(nx*xysize(0))
nys = long(ny*xysize(1))
temp = fltarr(nxs,nys) & temp(*)=0.0
istart = ([nxs,nys]-[nx,ny])/2L
iend = istart+[nx,ny]-1L
temp(istart(0):iend(0),istart(1):iend(1)) = bx2 & bx2 = temp
temp(istart(0):iend(0),istart(1):iend(1)) = by2 & by2 = temp
temp(istart(0):iend(0),istart(1):iend(1)) = bz2 & bz2 = temp
if tag_index(magin,'JZ',/quiet) GE 0 and $
tag_index(magin,'J_PXD',/quiet) GE 0 and NOT keyword_set(binup) then $
jz2 = poly_2d(mag.jz,pj,qj,1,nx,ny,missing=0.0)
if tag_index(magin,'J_TOTAL',/quiet) GE 0 and $
tag_index(magin,'J_PXD',/quiet) GE 0 and NOT keyword_set(binup) then $
jtot2 = poly_2d(mag.j_total,pj,qj,1,nx,ny,missing=0.0)
if NOT keyword_set(quiet) then begin
mag1 = add2str(mag1,['bx','by','bz','point'], $
bx2(istart(0):iend(0),istart(1):iend(1)), $
by2(istart(0):iend(0),istart(1):iend(1)), $
bz2(istart(0):iend(0),istart(1):iend(1)), $
mag.point,/quiet)
mag1.point.pix_size(*) = min(magin.point.pix_size)
vmgram,mag1,/xyz,/arrows,/inc,len=2,/helio, $
title='bff: Input Field'
endif
if n_elements(nz) EQ 0 then nz = 30
; These will hold the force-free field
bxp = fltarr(nxs,nys,nz)
byp = fltarr(nxs,nys,nz)
bzp = fltarr(nxs,nys,nz)
u = fltarr(nxs,nys) ; The Fourier variables
v = fltarr(nxs,nys)
; Compute the Fourier frequencies (cycles per pixel)
t = FINDGEN(nxs)/FLOAT(nxs) & nn=nxs/2 & t(nn+1:nxs-1) = t(nn+1:nxs-1) - 1.0
FOR j=0,nys-1 DO u(*,j) = t
t = FINDGEN(nys)/FLOAT(nys) & nn=nys/2 & t(nn+1:nys-1) = t(nn+1:nys-1) - 1.0
FOR i=0,nxs-1 DO v(i,*) = t
den = 2.*!pi*(u^2+v^2)
den(0,0)=1.
overden = 1.0/den
;overden(0,0) = 0.
fz = fft(bz2,-1)
fx = fft(bx2,-1)
fy = fft(by2,-1)
message,/info,'fz(0,0) = '+strcompress(string(fz(0,0)))+ $
' (net flux): forced to zero'
fz(0,0) = complex(0.0,0.0) ; The condition that the net flux is zero
;fx(0,0) = complex(0.0,0.0)
;fy(0,0) = complex(0.0,0.0)
; Compute the "best" value of alpha. The standard deviation gives an
; indication of how constant alpha really is.
;junk = check_math(1,1) ;Turn math errors off avoiding over and underflows
calpha = 2*!pi*complex(0.,1.)*(u(1:*)*fy(1:*)-v(1:*)*fx(1:*))/fz(1:*) ; Corrected formula from Gary
;junk = check_math(1,0) ; Turn math errors back on
calpha = float(calpha)
scalpha = stdev(calpha,mcalpha)
calpha = median(calpha) ;;mcalpha
if keyword_set(get_alpha) then begin
alpha = calpha
message,/info,strcompress('Using computed alpha of '+string(alpha)+ $
' +/- '+string(scalpha))
endif $
else if n_elements(alpha) EQ 1 then message,/info, strcompress( $
'Using alpha = '+string(alpha)+' ; Computed alpha is '+string(calpha)+ $
' +/- '+string(scalpha))
k = 4.*!pi*!pi*(u^2+v^2)-alpha^2
good = where(k GE 0.,countg)
bad = where(k LT 0.,countb)
if keyword_set(finite_energy) then begin
if countb GT 0 then k(bad) = -k(bad)
endif $
else begin
k = complex(k,fltarr(n_elements(k)))
;junk = check_math(1,1) ;Turn math errors off avoiding over and underflows
cff = complex(0.,1.)*2.0*!pi/sqrt(-k(1:*))
; TRM 2003 Nov 10, fixed bug in indexing of cff
cff = -cff * (u(1:*)*fx(1:*)+v(1:*)*fy(1:*))/fz(1:*) ; Corrected expression from Gary's paper
;cff = -cff(1:*) * (u(1:*)*fx(1:*)+v(1:*)*fy(1:*))/fz(1:*) ; Corrected expression from Gary's paper
;junk = check_math(1,0) ; Turn math errors back on
;cff = cff * (u*fx-v*fy)/fx ; Direct from Gary's paper (WRONG)
; get the "best" value of the constant cff by averaging all values
cff = float(cff) ; take it to be a real constant (it should be)
if countb GT 0 then $
cff = median(cff(bad)) $ ;;total(cff(bad))/n_elements(cff(bad)) $
;;cff = total(cff(bad))/n_elements(cff(bad)) $
else cff = 0.
a1 = (1.0 + complex(0.0,1.0)*cff)/2.0 ; Corrected expressions from Gary
a2 = (1.0 - complex(0.0,1.0)*cff)/2.0
endelse
k = SQRT(k)
if keyword_set(filter) then begin
noise = 0.50 * max(rebin(abs(fz),1,1))
filt = abs(fz)^2 / (abs(fz)^2+noise^2)
endif $
else filt = 1.0
if keyword_set(finite_energy) then $
print,'Force Free: Bad count = ',countb, ' of ',countb+countg
if countg GT 0 OR NOT keyword_set(finite_energy) then begin
junk = check_math(1,1) ;Turn math errors off avoiding over and underflows
FOR iz=0,Nz-1 DO BEGIN
f=complexarr(nxs,nys,/nozero)
z = zscale*float(iz)
f = -complex(0.,1.)*(k*u-alpha*v)*fz*exp(-k*z)*overden
if (countb GT 0) then begin
if (keyword_set(finite_energy)) then $
f(bad) = complex(0.) $
else begin
f2 = -complex(0.,1.)*(-k*u-alpha*v)*fz*exp(k*z)*overden
f(bad) = a1(bad)*f(bad) + a2(bad)*f2(bad)
endelse
endif
bxp(*,*,iz) = FFT(f*filt,1)
f = -complex(0.,1.)*(k*v+alpha*u)*fz*exp(-k*z)*overden
if (countb GT 0) then begin
if (keyword_set(finite_energy)) then $
f(bad) = complex(0.) $
else begin
f2 = -complex(0.,1.)*(-k*v+alpha*u)*fz*exp(k*z)*overden
f(bad) = a1(bad)*f(bad) + a2(bad)*f2(bad)
endelse
endif
byp(*,*,iz) = FFT(f*filt,1)
f = fz*exp(-k*z)
if (countb GT 0) then begin
if (keyword_set(finite_energy)) then begin
if (iz NE 0) then f(bad) = complex(0.)
endif $
else begin
f2 = fz*exp(k*z)
f(bad) = a1(bad)*f(bad) + a2(bad)*f2(bad)
endelse
endif
bzp(*,*,iz) = FFT(f*filt,1)
if NOT keyword_set(quiet) then $
message,/info,'Completed height ' $
+ strcompress(string(iz),/remove_all) + ': ' $
+ strcompress(string(z),/remove_all) + ' pixels (' $
+ strcompress(string(long(722.*z*min(mag.point.pix_size))), $
/remove_all) $
+' km)'
ENDFOR
junk = check_math(1,0)
if keyword_set(surround) then begin
; If surround is set, un-surrounds
bxp = bxp(istart(0):iend(0),istart(1):iend(1),0:nz-1)
byp = byp(istart(0):iend(0),istart(1):iend(1),0:nz-1)
bzp = bzp(istart(0):iend(0),istart(1):iend(1),0:nz-1)
endif
; Now the field has been computed and I need to put it back onto the
; grid it started on
if NOT keyword_set(quiet) then begin
mag3=add2str(mag3,['bx','by','bz'], $
bxp(*,*,0),byp(*,*,0),bzp(*,*,0),/quiet)
vmgram,mag3,/xyz,/arrows,/inc,len=2,/helio, $
title='bff: Output Field (heliographic)'
endif
nx = n_elements(mag.bz(*,0))
ny = n_elements(mag.bz(0,*))
if ((NOT keyword_set(no_heliographic)) OR keyword_set(rzdep)) and $
NOT keyword_set(image) then begin
bxpigrid = fltarr(nx,ny,nz)
bypigrid = fltarr(nx,ny,nz)
bzpigrid = fltarr(nx,ny,nz)
for iz=0,nz-1 do begin
bxpigrid(*,*,iz) = poly_2d(bxp(*,*,iz),pinv,qinv,1,nx,ny,missing=0.)
bypigrid(*,*,iz) = poly_2d(byp(*,*,iz),pinv,qinv,1,nx,ny,missing=0.)
bzpigrid(*,*,iz) = poly_2d(bzp(*,*,iz),pinv,qinv,1,nx,ny,missing=0.)
endfor
endif $
else begin
bxpigrid = fltarr(nx,ny,1)
bypigrid = fltarr(nx,ny,1)
bzpigrid = fltarr(nx,ny,1)
bxpigrid(*,*,0) = poly_2d(bxp(*,*,0),pinv,qinv,1,nx,ny,missing=0.)
bypigrid(*,*,0) = poly_2d(byp(*,*,0),pinv,qinv,1,nx,ny,missing=0.)
bzpigrid(*,*,0) = poly_2d(bzp(*,*,0),pinv,qinv,1,nx,ny,missing=0.)
endelse
if keyword_set(no_heliographic) then begin
bxp2 = bxpigrid
byp2 = bypigrid
bzp2 = bzpigrid
endif $
else if keyword_set(image) then begin
bxp2 = fltarr(nx,ny,nz)
byp2 = fltarr(nx,ny,nz)
bzp2 = fltarr(nx,ny,nz)
xform,bxpigrid,bypigrid,bzpigrid,bxp2,byp2,bzp2,trans=xf.helio2image
endif $
else begin
bxp2 = bxp
byp2 = byp
bzp2 = bzp
endelse
b2point = magin.point
b2point.pix_size(*) = min(magin.point.pix_size)/float(binup)
; If the mgrams are left in helio coords, p should be set to zero so
; that the field vectors plot in the correct orientation.
if NOT keyword_set(no_heliographic) then begin
b2point = add2str(b2point,'P_TRUE',b2point.p)
b2point.p = 0.0
endif
; Get the image plane fields to put in the structure
xform,bxpigrid(*,*,0),bypigrid(*,*,0),bzpigrid(*,*,0),bxpi,bypi,bzpi, $
trans=xf.helio2image
btrans = sqrt(bxpi^2+bypi^2)
check = check_math(0,1)
bazim = !radeg*atan(-bxpi,bypi)
check = check_math(1,0)
;x = lindgen(nx,ny) MOD nx
;y = long(lindgen(nx,ny) / nx)
; Add the fields to the mag structure
mag2=add2str(mag2,['bxp','byp','bzp','zscale','alpha','bx','by','bz', $
'b_long','b_trans','b_azim','b_pxd','b_pyd',$
'pix_size','point','bx_obs','by_obs','bz_obs'], $
bxp2,byp2,bzp2,zscale,alpha,bxp2(*,*,0),byp2(*,*,0), $
bzp2(*,*,0),bzpi,btrans,bazim, $
mag.b_pxd,mag.b_pyd,min(b2point.pix_size), $
b2point,mag.bx,mag.by,mag.bz,/quiet)
if tag_index(magin,'I_CONT',/quiet) GE 0 then $
mag2 = add2str(mag2,['i_cont'],magin.i_cont,/quiet)
if tag_index(magin,'J_PXD',/quiet) GE 0 and $
tag_index(magin,'J_PYD',/quiet) then $
mag2 = add2str(mag2,['j_pxd','j_pyd'],magin.j_pxd,magin.j_pyd,/quiet)
if n_elements(jz2) GT 0 then mag2 = add2str(mag2,['jz'],jz2,/quiet)
if n_elements(jtot2) GT 0 then mag2 = add2str(mag2,['j_total'],jtot2,/quiet)
if NOT keyword_set(quiet) and $
(keyword_set(no_heliographic) OR keyword_set(image)) then begin
vmgram,mag2,/xyz,/arrows,/inc,len=2, $
title='bff: Output Field (image)'
endif
endif $
else begin
; All the fourier components have been destroyed
print,'ERROR: No good k values left - alpha is too large.'
stop
endelse
return,mag2
end