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