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