pro a_jz_usage
;+
;
;	procedure:  a_jz
;
;	purpose:  compute zenith current Jz & d(Bz)/dz from a_* files
;
;	routines:  a_jz_usage  vec_deriv  asp_deriv  a_jz
;
;	author:  paul@ncar  4/98
;
;	credits:  A. Skumanich, M. Semel; papers & code examples
;
;	reference:
;		Semel, M.; A. Skumanich: An ambiguity-free determination
;		of Jz in solar active regions, Astromomy and Astrophysics.,
;		331, 383-391, 1998
;
;		Skumanich, A.; M. Semel: Field azimuth disambiguation using
;		ambiguity-free currents, Sol. Phys., 164, 291-302, 1996
;
;==============================================================================
;
;       Check number of parameters.
;
if 1 then begin
	print
	print, "usage:	a_jz, jzl, jzt [, dir]"
	print
	print, "	Compute zenith current & d(Bz)/dz from a_* files."
	print
	print, "	Arguments"
	print, "		dir	- a_* file directory path"
	print, "			  (def: current working directory)."
	print, "		jzl	- output z current, line of sight"
	print, "			  component, G/Mm (a_* vector format)."
	print, "		jzt	- output z current, plane of sky"
	print, "			  component, G/Mm (a_* vector format)."
	print, "	Keywords"
	print, "		dbzl	- output d(Bz)/dz, line of sight"
	print, "			  component, G/Mm (a_* vector format)."
	print, "		dbzt	- output d(Bz)/dz, plane of sky"
	print, "			  component, G/Mm (a_* vector format)."
	print, "		nofill	- set to disable fill factor applied"
	print, "			  to field (def: apply fill factor)."
	print
	print, "		bt_cutoff  - use most continuous derivatives"
	print, "			  where surface tranverse field"
	print, "			  strength is less than bt_cutoff"
	print, "			  (def: 10. Gauss)."
	print
	print, "		continuous  - for all points use locally the"
	print, "			  most continuous azimuth resolution"
	print, "			  to take derivatives."
	print, "			  (def: use rotated current formula)."
	print
	print, "		a_str	- input a_image.pro or azam structure"
	print, "			  to return 2D arrays"
	print, "			  (def: return a_* file format)."
	print
	return
endif
;-
end
;------------------------------------------------------------------------------
;
; Derivation of rotated ambiguity-free current formula:
;
; Let the rotated frame have coordinates (x',y').  The transverse to line
; of sight field along the surface in the rotated frame is (bx',by').
; Ambiguous field is (-bx',-by').   Note derivatives d(bx'bx')/dy' and
; d(by'by')/dx' below are independent of azimuth resolution, called
; ambiguity-free.  Note the ambiguous current is -j; j is said to have
; parity sign.  Goal is to ratate to avoid division by small bx' or by'.
; Later bx'=by' is set, optimum case.
;
; j = d(by')/dx'-d(bx')/dy'
;
; j = ( 1/by' d(by'by')/dx' - 1/bx' d(bx'bx')/dy' )/2
;
; j = ( bx' d(by'by')/dx' - by' d(bx'bx')/dy' )/(2*bx'by')
;
; x' = xC-yS   x =  x'C+y'S   S=sin(a)
; y' = xS+yC   y = -x'S+y'C   C=cos(a)
;
; bx'= bxC-byS
; by'= bxS+byC
;
; j = ( (bxC-byS) d/dx' (bxbxSS+2bxbySC+bybyCC)
;     - (bxS+byC) d/dy' (bxbxCC-2bxbySC+bybySS) )
;     / 
;     ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; j = ( (bxC-byS) ( d/dx (bxbxSS+2bxbySC+bybyCC) C
;                 - d/dy (bxbxSS+2bxbySC+bybyCC) S )
;     - (bxS+byC) ( d/dx (bxbxCC-2bxbySC+bybySS) S
;                 + d/dy (bxbxCC-2bxbySC+bybySS) C ) )
;     / 
;     ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; j = 
; ( bx ( d/dx ( ( CCSS-CCSS)bxbx + 2( CCCS+CSSS)bxby + ( CCCC-SSSS)byby )
;      + d/dy ( (-CSSS-CCCS)bxbx + 2(-CCSS+CCSS)bxby + (-CCCS-CSSS)byby ) )
; - by ( d/dx ( ( CSSS+CCCS)bxbx + 2( CCSS-CCSS)bxby + ( CCCS+CSSS)byby )
;      + d/dy ( (-SSSS+CCCC)bxbx + 2(-CSSS-CCCS)bxby + (-CCSS+CCSS)byby ) ) )
; / 
; ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; j =
; ( bx ( 2SC d(bxby)/dx + (CC-SS) d(byby)/dx - SC d(bxbx)/dy - SC d(byby)/dy )
; + by ( 2SC d(bxby)/dy + (SS-CC) d(bxbx)/dy - SC d(byby)/dx - SC d(bxbx)/dx ) )
; / 
; ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; set: a = pi/4 - atan(by,bx) so the field rotates to pi/4.
;
; notes:
;
; (1) With 'a' set as above
;
;	S = sin(a) --> ( bx - by )/( bt sqrt(2) )    bt = sqrt(bxbx+byby)
;	C = cos(a) --> ( bx + by )/( bt sqrt(2) )
;
; (2) All terms with S & C in j are ambiguity-free, do not depend on azimuth
;     resolution.
;
; (3) The denominator with 'a' above can be written
;
;	2 (SCbxbx+(CC-SS)bxby-SCbyby) --> ((bxbx+byby)*(bxbx+byby))/(bt*bt)
;
;     The bt*bt devides out with the numerator.  (bxbx+byby)^2 is
;     identical to the denominator that stabilizes Semel Skumanich
;     construct.   One can make whatever manipulations needed to get
;     the numerator in Semel Skumanich form; see (4).
;
; (4) Following is a decomposition of j above, recovering the unrotated
;     current formula.  Semel Skumanich construct decomposes in a
;     similar way.
;
; j =
; ( bx ( 2SC d(bxby)/dx + (CC-SS) d(byby)/dx - SC d(bxbx)/dy - SC d(byby)/dy )
; + by ( 2SC d(bxby)/dy + (SS-CC) d(bxbx)/dy - SC d(byby)/dx - SC d(bxbx)/dx ) )
; /
; ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; j = ( bx ( 2SC (bx d(by)/dx + by d(bx)/dx ) + (CC-SS) 2 by d(by)/dx
;     - SC 2 bx d(bx)/dy - SC 2 by d(by)/dy )
;     + by ( 2SC (by d(bx)/dy + bx d(by)/dy ) + (SS-CC) 2 bx d(bx)/dy
;     - SC 2 by d(by)/dx - SC 2 bx d(bx)/dx ) )
;     /
;     ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; j = ( d(bx)/dx ( 2 SC bxby - 2 SC bxby )
;     + d(by)/dx ( 2 SC bxbx - 2 SC byby + 2 (CC-SS) bxby )
;     + d(bx)/dy ( 2 SC byby - 2 SC bxbx + 2 (SS-CC) bybx )
;     + d(by)/dy ( 2 SC bybx - 2 SC bybx ) )
;     /
;     ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; j = d(by)/dx - d(bx)/dy
;
;------------------------------------------------------------------------------
;Following is the derivation of d(Bz)/dz in a rotated cooridinate frame
;using ambiguity free derivatives.   See above for some definitions.
;
; k = -d(bx')/dx'-d(by')/dy' = d(Bz)/dz
;
; k = -( 1/bx' d(bx'bx')/dx' + 1/by' d(by'by')/dy' )/2
;
; k = -( by' d(bx'bx')/dx' + bx' d(by'by')/dy' )/(2*bx'by')
;
; k = -( (bxS+byC) d/dx' (bxbxCC-2bxbySC+bybySS)
;      + (bxC-byS) d/dy' (bxbxSS+2bxbySC+bybyCC) )
;      /
;      ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; k = -( (bxS+byC) ( d/dx (bxbxCC-2bxbySC+bybySS) C
;                  - d/dy (bxbxCC-2bxbySC+bybySS) S )
;      + (bxC-byS) ( d/dx (bxbxSS+2bxbySC+bybyCC) S
;                  + d/dy (bxbxSS+2bxbySC+bybyCC) C ) )
;      /
;      ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; k =
; -( bx ( d/dx ( ( CCCS+CSSS)bxbx + 2(-CCSS+CCSS)bxby + ( CSSS+CCCS)byby )
;       + d/dy ( (-CCSS+CCSS)bxbx + 2( CSSS+CCCS)bxby + (-SSSS+CCCC)byby ) )
;  + by ( d/dx ( ( CCCC-SSSS)bxbx + 2(-CCCS-CSSS)bxby + ( CCSS-CCSS)byby )
;       + d/dy ( (-CCCS-CSSS)bxbx + 2( CCSS-CCSS)bxby + (-CSSS-CCCS)byby ) ) )
;  /
;  ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; k =
; -( bx ( CS d(bxbx)/dx + CS d(byby)/dx + 2CS d(bxby)/dy + (CC-SS) d(byby)/dy )
;  + by (-CS d(byby)/dy - CS d(bxbx)/dy - 2CS d(bxby)/dx + (CC-SS) d(bxbx)/dx )
;  )
;  /
;  ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; k =
; -( bx ( 2bx CS d(bx)/dx + 2by CS d(by)/dx
;  + 2CS (bx d(by)/dy + by d(bx)/dy ) + 2by (CC-SS) d(by)/dy )
;  + by (-2by CS d(by)/dy - 2bx CS d(bx)/dy
;  - 2CS (by d(bx)/dx + bx d(by)/dx ) + 2bx (CC-SS) d(bx)/dx ) )
;  /
;  ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; k = -( d(bx)/dx ( 2 SC bxbx - 2 SC byby + 2 (CC-SS) bxby )
;      + d(by)/dx ( 2 SC bxby - 2 SC bxby )
;      + d(bx)/dy ( 2 SC bxby - 2 SC bxby )
;      + d(by)/dy ( 2 SC bxbx - 2 SC byby + 2 (CC-SS) bxby ) )
;      /
;      ( 2(SCbxbx+(CC-SS)bxby-SCbyby) )
;
; k = -d(bx)/dx - d(by)/dy
;
;------------------------------------------------------------------------------
;
;	procedure:  vec_deriv
;
;	purpose:  compute derivatives of vector from ASP array row or column
;		  using IDL 'deriv.pro'.
;
;------------------------------------------------------------------------------
pro vec_deriv, vdim, vmask, dat, dlt

				    ;Make sure vectors are vectors.
vmask = reform( vmask, vdim )
dat   = reform(   dat, vdim )
				    ;Mask extend mask one invalid and one
				    ;valid point used to terminate loop.
mask = [ vmask, [ 0, 1 ] ]
				    ;Initial derivative vector.
dlt = 0.*dat
				    ;Start is one less than beginning.
v1 = -1
				    ;Loop over continuous data sequences.
while 1 do begin
				    ;Get first point of data sequence.
	v0 = v1+1 + ( where( mask(v1+1:*) ) )(0)
	if v0 ge vdim then  return
				    ;Get last data point in sequence.
	v1 = v0 + ( where( mask(v0:*) eq 0 ) )(0) - 1

				    ;Set derivatives if possible
				    ;or set mask out.
	case 1 of
	v1-v0 ge 2:  dlt(v0:v1) = deriv(dat(v0:v1))
	v1-v0 eq 1:  dlt(v0:v1) = dat(v1)-dat(v0)
	else:        vmask(v0) = 0
	end
end

end
;------------------------------------------------------------------------------
;
;	procedure:  asp_deriv
;
;	purpose:  compute derivatives of ASP array using IDL 'deriv.pro'.
;
;------------------------------------------------------------------------------
pro asp_deriv, xdim, ydim, mask0, mask, dat, dati, datj

				    ;Initial mask available derivatives.
mask = mask0
				    ;Initial derivative arrays.
dati = ( datj = 0.*dat )
				    ;Loop over rows, take x derivatives.
for j=0,ydim-1 do begin
	vmask = mask0(*,j)
	vec_deriv, xdim, vmask, dat(*,j), dlt
	dati(*,j) = dlt
	mask(*,j) = vmask
end
				    ;Loop over columns, take y derivatives.
for i=0,xdim-1 do begin
	vmask = mask0(i,*)
	vec_deriv, ydim, vmask, dat(i,*), dlt
	datj(i,*) = dlt
	mask(i,*) = vmask and mask(i,*)
end
				    ;Mask to eliminate derivatives not
				    ;availabe in both directions.
dati = dati*mask
datj = datj*mask

end
;------------------------------------------------------------------------------
;
;	procedure:  a_jz
;
;	purpose:  compute vertical to surface current from a_* files
;
;------------------------------------------------------------------------------
pro a_jz, jzl, jzt, dty_, a_str=a_str, dbzl=dbzl, dbzt=dbzt $
, nofill=nofill, bt_cutoff=bt_cutoff, continuous=continuous_, rotate=rotate_

				    ;Check number of parameters.
if n_params() eq 0 then begin
	a_jz_usage
	return
end
                                    ;Directory path; append with /.
if  n_elements(dty_) ne 0  then  dty=dty_  else  dty=''
if  dty ne ''  then  if  strmid(dty,strlen(dty)-1,1) ne '/'  then  dty=dty+'/'

				    ;Check if continuous derivatives are
				    ;used everywhere or coordinate rotation
				    ;method is used.
continuous = keyword_set(continuous_)
rotate     = keyword_set(rotate_    )

				    ;Cutoff min for tranverse field.
if  n_elements(bt_cutoff) ne 0 $
then  bt4cut = double(bt_cutoff)^4. $
else  bt4cut = 10.D^4.
				    ;Read field vectors & pointing info.
fld   = double( azam_a_image( dty+'a_fld'       , b ))
azm   = double( azam_a_image( dty+'a_azm'       , b ))*(!dpi/180.)
psi   = double( azam_a_image( dty+'a_psi'       , b ))*(!dpi/180.)
azm1  = double( azam_a_image( dty+'a_1azm'      , b ))*(!dpi/180.)
azm2  = double( azam_a_image( dty+'a_2azm'      , b ))*(!dpi/180.)
inc2  = double( azam_a_image( dty+'a_2incl'     , b ))*(!dpi/180.)
inc1  = double( azam_a_image( dty+'a_1incl'     , b ))*(!dpi/180.)
alpha = double( azam_a_image( dty+'a_alpha'     , b ))
la    = double( azam_a_image( dty+'a__latitude' , b ))*(!dpi/180.)
ln    = double( azam_a_image( dty+'a__longitude', b ))*(!dpi/180.)

				    ;Apply fill factor on keyword.
;if not keyword_set(nofill) then  fld=fld*(1.-alpha)
;  for Hinode, HMI, default has alpha as fill fract, not scatt light
if not keyword_set(nofill) then  fld=fld*alpha

				    ;Ratio (great circle)/longitude.
cosla = cos(la)
				    ;Field transverse line of sight.
tmp = fld*sin(psi)
losbx = tmp*cos(azm)
losby = tmp*sin(azm)

sin_inc1=sin(inc1)  &  sin_inc2=sin(inc2)
sin_azm1=sin(azm1)  &  sin_azm2=sin(azm2)
cos_azm1=cos(azm1)  &  cos_azm2=cos(azm2)

				    ;Projection los component field on surface.
bx  = fld*( sin_inc1*cos_azm1 + sin_inc2*cos_azm2 )/2.
by  = fld*( sin_inc1*sin_azm1 + sin_inc2*sin_azm2 )/2.

				    ;Projection transverse field on surface.
bxt = fld*( sin_inc1*cos_azm1 - sin_inc2*cos_azm2 )/2.
byt = fld*( sin_inc1*sin_azm1 - sin_inc2*sin_azm2 )/2.

				    ;Data masks.
pmask0=lonarr(b.xdim,b.ydim)  &  pmask0(b.pxy)=1
smask0=lonarr(b.xdim,b.ydim)  &  smask0(b.sxy)=1

undefine_sbsp, azm, psi, alpha, fld, tmp $
, azm1, inc1, sin_inc1, sin_azm1, cos_azm1 $
, azm2, inc2, sin_inc2, sin_azm2, cos_azm2

;------------------------------------------------------------------------------
				    ;Calculate Jzl.

				    ;Calculate derivatives.
asp_deriv, b.xdim, b.ydim, pmask0, pmask, ln, lni, lnj
asp_deriv, b.xdim, b.ydim, pmask0, pmask, la, lai, laj
asp_deriv, b.xdim, b.ydim, smask0, smask, bx, bxi, bxj
asp_deriv, b.xdim, b.ydim, smask0, smask, by, byi, byj

				    ;Jacobian transformation
				    ;(background set to 1.)
jd = (lni*laj-lnj*lai) + (1-pmask)

jbxln = bxi*laj-bxj*lai
jbyln = byi*laj-byj*lai

jbxla = bxi*lnj-bxj*lni
jbyla = byi*lnj-byj*lni
				    ;Calc derivatives.
bydx  =  jbyln/(jd*cosla)
bxdx  =  jbxln/(jd*cosla)
bxdy  = -jbxla/jd
bydy  = -jbyla/jd
				    ;Current & d(Bz)/dz.     
jzl   = bydx-bxdy
dbzl  = -bxdx-bydy

undefine_sbsp, jbxln, jbxla, bxdx, bydx, bx, bxi, bxj, pmask0
undefine_sbsp, jbyln, jbyla, bydy, bxdy, by, byi, byj, pmask

;------------------------------------------------------------------------------
				    ;Compute jzt with rotated current formula.

				    ;Squares and cross products.
bxt2  = bxt*bxt
byt2  = byt*byt
bxtyt = bxt*byt

bt4   = (bxt2+byt2)^2.
				    ;Calculate derivatives.
asp_deriv, b.xdim, b.ydim, smask0, smask, bxt2,  bxxi, bxxj
asp_deriv, b.xdim, b.ydim, smask0, smask, byt2,  byyi, byyj
asp_deriv, b.xdim, b.ydim, smask0, smask, bxtyt, bxyi, bxyj

				    ;Jacobian transformation.
jbxt2ln = bxxi*laj-bxxj*lai
jbyt2ln = byyi*laj-byyj*lai
jbxytln = bxyi*laj-bxyj*lai

jbxt2la = bxxi*lnj-bxxj*lni
jbyt2la = byyi*lnj-byyj*lni
jbxytla = bxyi*lnj-bxyj*lni
				    ;Calc derivatives.
bxt2x   =  jbxt2ln/(jd*cosla)
byt2x   =  jbyt2ln/(jd*cosla)
bxytx   =  jbxytln/(jd*cosla)

bxt2y   = -jbxt2la/jd
byt2y   = -jbyt2la/jd
bxyty   = -jbxytla/jd
				    ;Set trig except for bt*bt normalization
				    ;that devides out.  This is a rotation by
				    ;theta = pi/4 - atan( byt, bxt )
sncn    = (bxt2-byt2)/2.
cn2_sn2 = 2.*bxtyt

jzt = ( bxt * ( 2.*sncn*bxytx + cn2_sn2*byt2x - sncn*bxt2y - sncn*byt2y ) $
      + byt * ( 2.*sncn*bxyty - cn2_sn2*bxt2y - sncn*byt2x - sncn*bxt2x ) $
      ) $
      /( bt4cut > bt4 )
				    ;Semel Skumanich construct (COMMENTED OUT)
;jzt = ( bxt * ( bxt2*bxytx - .5*bxtyt*(bxt2x-byt2x) -.5*(bxt2+byt2)*bxt2y ) $
;      - byt * ( byt2*bxyty - .5*bxtyt*(byt2y-bxt2y) -.5*(bxt2+byt2)*byt2x ) $
;      ) $
;      /( bt4cut > bt4 )
				    ;d(Bz)/dz.
dbzt = -( bxt * (  sncn*bxt2x + sncn*byt2x + 2*sncn*bxyty + cn2_sn2*byt2y ) $
        + byt * ( -sncn*byt2y - sncn*bxt2y - 2*sncn*bxytx + cn2_sn2*bxt2x ) $
	) $
	/( bt4cut > bt4 )
      
undefine_sbsp, jbxt2ln, jbxt2la, bxt2x, bxt2y, bxt2,  bxxi, bxxj, sncn
undefine_sbsp, jbyt2ln, jbyt2la, byt2x, byt2y, byt2,  byyi, byyj, cn2_sn2
undefine_sbsp, jbxytln, jbxytla, bxytx, bxyty, bxtyt, bxyi, bxyj

;------------------------------------------------------------------------------
				    ;Now go over a point at a time using
				    ;the continuous method.
for j=0,b.ydim-1 do begin
for i=0,b.xdim-1 do begin
if  smask(i,j)  and  ( continuous  or  bt4(i,j) lt bt4cut )  then begin

				    ;Subscripting for vectors in (x,y).
				    ;Size 5 where available.
	x0   = 0 > (i-2)
	x1   = (i+2) < (b.xdim-1)
	xdim = x1-x0+1
	xc   = i-x0
	imask  = smask0(x0:x1,j)

	y0   = 0 > (j-2)
	y1   = (j+2) < (b.ydim-1)
	ydim = y1-y0+1
	yc   = j-y0
	jmask  = smask0(i,y0:y1)
				    ;Try to reduce to size 3 vectors.
	if  xc gt 0  and  xc lt xdim-1  then begin
	if  total( imask(xc-1:xc+1) eq 3 )  then begin
		xc    = 1
		x0    = i-1
		x1    = i+1
		xdim  = 3
		imask = [1,1,1]
	end
	end

	if  yc gt 0  and  yc lt ydim-1  then begin
	if  total( jmask(yc-1:yc+1) eq 3 )  then begin
		yc    = 1
		y0    = j-1
		y1    = j+1
		ydim  = 3
		jmask = [1,1,1]
	end
	end
				    ;Get sub vectors.
	ibxt   =   bxt(x0:x1,j)
	ibyt   =   byt(x0:x1,j)

	jbxt   =   bxt(i,y0:y1)
	jbyt   =   byt(i,y0:y1)
				    ;Get sub vectors.
	ilosbx = losbx(x0:x1,j)
	ilosby = losby(x0:x1,j)

	jlosbx = losbx(i,y0:y1)
	jlosby = losby(i,y0:y1)

				    ;From azimuth continuity in line of sight
				    ;frame, resolve azimuth to center pixel.
	whr = where( ilosbx*ilosbx(xc)+ilosby*ilosby(xc) lt 0. )
	if whr(0) ge 0 then begin
		ibxt(whr) = -ibxt(whr)
		ibyt(whr) = -ibyt(whr)
	end
	whr = where( jlosbx*jlosbx(yc)+jlosby*jlosby(yc) lt 0. )
	if whr(0) ge 0 then begin
		jbxt(whr) = -jbxt(whr)
		jbyt(whr) = -jbyt(whr)
	end
				    ;Take derivatives.
	vec_deriv, xdim, imask, ibxt, bxti  &  bxti = bxti(xc)
	vec_deriv, xdim, imask, ibyt, byti  &  byti = byti(xc)

	vec_deriv, ydim, jmask, jbxt, bxtj  &  bxtj = bxtj(yc)
	vec_deriv, ydim, jmask, jbyt, bytj  &  bytj = bytj(yc)

				    ;Jacobian transformation.
	jbxtln   = bxti*laj(i,j)-bxtj*lai(i,j)
	jbytln   = byti*laj(i,j)-bytj*lai(i,j)

	jbxtla   = bxti*lnj(i,j)-bxtj*lni(i,j)
	jbytla   = byti*lnj(i,j)-bytj*lni(i,j)

				    ;Calc derivatives.
	bytx     =  jbytln/(jd(i,j)*cosla(i,j))
	bxtx     =  jbxtln/(jd(i,j)*cosla(i,j))
	bxty     = -jbxtla/jd(i,j)
	byty     = -jbytla/jd(i,j)


				    ;Current & d(Bz)/dz.     
	jzt(i,j) = bytx-bxty
	dbzt(i,j)= -bxtx-byty
end
end
end

;------------------------------------------------------------------------------
				    ;Current G/Mm, a_* file format.
jzl  = float(  jzl(b.sxy)/700. )
jzt  = float(  jzt(b.sxy)/700. )
dbzl = float( dbzl(b.sxy)/700. )
dbzt = float( dbzt(b.sxy)/700. )
				    ;Return 2D arrays on keyword structure.
if  n_elements(a_str) eq 0  then  return

tmp = jzl
jzl = fltarr(a_str.xdim,a_str.ydim)
jzl(a_str.sxy) = tmp(a_str.vec_sxy)

tmp = jzt
jzt = fltarr(a_str.xdim,a_str.ydim)
jzt(a_str.sxy) = tmp(a_str.vec_sxy)

tmp = dbzl
dbzl = fltarr(a_str.xdim,a_str.ydim)
dbzl(a_str.sxy) = tmp(a_str.vec_sxy)

tmp = dbzt
dbzt = fltarr(a_str.xdim,a_str.ydim)
dbzt(a_str.sxy) = tmp(a_str.vec_sxy)

end