;+ ;NAME: ; GET_ALPHA ;PURPOSE: ; Determines the best value of alpha for a force-free extrapolation by ; checking the resulting error in the photospheric transverse field. ; units are PIX^-1 ;CATEGORY: ;CALLING SEQUENCE: ; best_alpha = get_alpha(mag) ;INPUTS: ; mag = magnetic field structure ;OPTIONAL INPUT PARAMETERS: ;KEYWORD PARAMETERS ; /fffalpha = return the value of alpha from fff. *Much* faster. ; Bx and By are required, so, if they do not exist in ; the structure, they are computed using a LFF field ; to crudely resolve the ambiguity. ; alimit = the max alpha (in abs) to check. ; astep = the step size in alpha for the search. ; avals = (Output) The alpha values checked. ; evals = (Output) The average residual at each of the alpha values. ; nevals = (Output) The normalized residual at each of the alpha values. ; (see bxerr, byerr below) ; sigevals = (Output) uncertainty of evals (needs bxerr,byerr) ; signevals = (Output) uncertainty of nevals (needs bxerr,byerr) ; /resolve = resolve 180 degree ambiguity at every alpha. This is the ; default if tag BX, BY, and BZ are NOT present in the mag ; structure (if BX, BY, and BZ are present, get_alpha assumes ; the ambiguity has been resolved already). This ; keyword should not be set if the 180 ambiguity has ; already been resolved. ; btcut = Compute the residual only at those points with b_transv GE btcut ; xrange = Compute the residual only within xrange [x1,x2]. ; yrange = Compute the residual only within yrange [y1,y2]. ; mask = array of same size as the field arrays. If the mask is 1 ; then the pixels is used, if zero it is neglected. ; /xyz = Use heliographic field for comparison. ; guard = passed on to fff.pro. ; /noplot = do not plot the residual vs. alpha. ; /quiet = work quietly. ; /verbose = work loudly. ; bxerr,byerr = uncertainty in bx, by ; sigalpha = (output) uncertainty in alpha, in units same as alpha (pix^-1), ; taking into account the residuals in the region of the minimum... ; bzcut = only works if /xyz is set, then uses Bz rather than Bhoriz for cut. ; ;OUTPUTS: ;COMMON BLOCKS: ;SIDE EFFECTS: ;RESTRICTIONS: ;PROCEDURE: ; Computes the force-free field from LOS field, then computes ; total(sqrt((ff.byi-byi)^2.0d0+(ff.bxi-ff.bxi)^2)) ; as the residual. Computes best alpha from a parabolic fit to the residual ; with line_cent.pro. ;MODIFICATION HISTORY: ; T. Metcalf December 1992 ; TRM 15-Apr-94 Converted find_alpha (using bff.pro) to get_alpha (using ; fff.pro). ; KDL 1998, modified to use normalized nsigma residuals ("errors"), ; KDL 2001 modified (finally!) to return the uncertainty in alpha ; given the magnitude of the residuals and the uncertainties in ; them. ; TRM 04-Apr-03 Added /shear keyword. ; TRM 21-Sep-04 Added KD's error checking. ; TRM 07-Sep-05 /resolve forces resolution of 180 ambiguity even ; if BX,BY,BZ are present. Also added guard ; keyword. ; TRM 12-May-06 Added /fffalpha keyword. ;- function get_alpha,mag,alimit=alimit,astep=astep,resolve=resolve,xyz=xyz, $ avals=avals,evals=evals,nevals=nevals,noplot=noplot,$ btcut=btcut,xrange=xrange,yrange=yrange,$ quiet=quiet,verbose=verbose,shear=shear,window=win, $ bxerr=bxerr,byerr=byerr,bzcut=bzcut,$ sigevals=sigevals,signevals=signevals,sigalpha=sigalpha, $ guard=guard,mask=maskin,fffalpha=fffalpha,_extra=_extra common get_alpha_private,plotwindow if n_elements(win) EQ 1 then plotwindow = win nx = n_elements(mag.b_long(*,0)) ny = n_elements(mag.b_long(0,*)) if n_elements(alimit) NE 1 then alimit=1.0d0 alimit=abs(double(alimit)) if n_elements(astep) NE 1 then astep=0.02d0 astep=abs(double(astep)) if NOT keyword_set(maskin) then mask = make_array(nx,ny,/integer,value=1) $ else mask = fix(reform(maskin,nx,ny)) if n_elements(btcut) NE 1 then btcut=0.0 if n_elements(bxerr) LT 1 then bxerr=btcut > 1. if n_elements(byerr) LT 1 then byerr=btcut > 1. if n_elements(xrange) NE 2 then xrange=[0,nx-1] xrange = abs(fix(xrange)) if n_elements(yrange) NE 2 then yrange=[0,ny-1] yrange = abs(fix(yrange)) x = lindgen(nx,ny) y = x / nx x = x MOD nx ;if (keyword_set(resolve) AND (NOT keyword_set(xyz))) OR $ if (keyword_set(resolve)) OR $ (tag_index(mag,'BX') LT 0 OR $ tag_index(mag,'BY') LT 0 OR $ tag_index(mag,'BZ') LT 0) then qresolve=1b $ else qresolve=0b alpha = -alimit avals = [0.0d0] evals = [0.0d0] nevals = [0.0d0] sigevals = [0.0d0] signevals = [0.0d0] p = mag.point.p b0 = mag.point.b0 if tag_index(mag,'CMD') GE 0 then lc = mag.cmd*!dtor $ else lc = mag.point.cmd if tag_index(mag,'LATITUDE') GE 0 then bc = mag.latitude*!dtor $ else bc = mag.point.lat pixsize = mag.point.pix_size if keyword_set(xyz) then begin ;good = where(sqrt(mag.bx^2+mag.by^2+mag.bz^2) GE abs(btcut) $ good = where(sqrt(mag.bx^2+mag.by^2) GE abs(btcut) $ AND x GE xrange(0) and x LE xrange(1) $ AND y GE yrange(0) and y LE yrange(1) $ AND mask NE 0) if keyword_set(bzcut) then begin good = where(abs(mag.bz) GE abs(bzcut) $ AND x GE xrange(0) and x LE xrange(1) $ AND y GE yrange(0) and y LE yrange(1) $ AND mask NE 0) endif endif else begin good = where(abs(mag.b_trans) GE abs(btcut) $ AND x GE xrange(0) and x LE xrange(1) $ AND y GE yrange(0) and y LE yrange(1) $ AND mask NE 0) endelse nbtrans = n_elements(good) if good(0) LT 0 then begin message,/info,'All data points eliminated in field cut' return,0.0 endif if keyword_set(fffalpha) then begin ; Get a fast value for alpha from fff.pro if (tag_index(mag,'BX') LT 0 OR $ tag_index(mag,'BY') LT 0 OR $ tag_index(mag,'BZ') LT 0) then begin ; Resolve the ambiguity with the potential field pot = fff(mag,_extra=_extra,/quiet,normal=0,guard=guard,/image) pxi = pot(*,*,0) ; These are image coordinates pyi = pot(*,*,1) pa_azimuth = atan(-pxi,pyi)/!dtor ; Az is CCW from North magr = fix_azimuth(mag,pa_azimuth) ; Now compute alpha and disambiguate again ralpha = get_alpha(magr,/fffalpha,/quiet,btcut=btcut,bzcut=bzcut,xyz=xyz, $ xrange=xrange,yrange=yrange,mask=mask,guard=guard, $ verbose=verbose) ; recursive call, make sure bx,by,bz are defined! pot = fff(mag,_extra=_extra,/quiet,guard=guard,/image,alpha=ralpha) pxi = pot(*,*,0) ; These are image coordinates pyi = pot(*,*,1) pa_azimuth = atan(-pxi,pyi)/!dtor ; Az is CCW from North magr = fix_azimuth(mag,pa_azimuth) endif else magr = mag bad = float(magr.bx) bad[*] = 0. bad[good] = 1. f = fff(magr,computed_alpha=calpha,_extra=_extra,/normal,/quiet,/get_alpha, $ guard=guard,bx=magr.bx*bad,by=magr.by*bad,sigalpha=sigalpha) calpha = calpha * sqrt((nx*ny)/nbtrans) ; correct for sub-range return,calpha ; skip the rest of the slow code endif if keyword_set(xyz) then begin azim = atan(-mag.bx,mag.by)/!dtor bxi = mag.bx byi = mag.by endif else begin azim = (mag.b_azim MOD 360.) bxi = -mag.b_trans*sin(azim*!dtor) byi = mag.b_trans*cos(azim*!dtor) endelse junk=check_math(/print) ; clear error status if NOT keyword_set(verbose) AND NOT keyword_set(quiet) then $ progress,0.0,label='get_alpha',/reset except = !except !except=0 REPEAT begin if keyword_set(xyz) then begin pot=fff(mag.bz,/normal,/image,b0=b0,pang=p,cmd=lc,lat=bc,pixel=pixsize, $ /quiet,alpha=alpha,bx=mag.bx,by=mag.by,guard=guard) endif else begin pot=fff(mag.b_long,b0=b0,pang=p,cmd=lc,lat=bc,pixel=pixsize,/quiet, $ alpha=alpha,/image,btrans=mag.b_trans,azim=azim,guard=guard) endelse pxi = pot(*,*,0) ; These are image coordinates pyi = pot(*,*,1) ; Fix the azimuthal angle from the force-free field. pa_azimuth = atan(-pxi,pyi)/!dtor ; Az is CCW from North azim = mag.b_azim if qresolve then begin ;pa_azimuth = atan(-pxi,pyi)/!dtor ; Az is CCW from North ;azim = mag.b_azim fixa = where(abs(arange(azim,pa_azimuth)-pa_azimuth) GT $ 90.,count) if (count gt 0) then azim(fixa) = azim(fixa) + 180. azim = (azim MOD 360.) bxi = -mag.b_trans*sin(azim*!dtor) byi = +mag.b_trans*cos(azim*!dtor) endif if good(0) GE 0 then begin if keyword_set(shear) then begin ;Check the error from the field weighted shear angle error = total(mag.b_trans[good]*abs(arange(azim(good),pa_azimuth(good))-pa_azimuth(good)))/total(mag.b_trans[good]) sigerr = stdev(mag.b_trans[good]*abs(arange(azim(good),pa_azimuth(good))-pa_azimuth(good)))/sqrt(total(mag.b_trans[good])) nerror = total(mag.b_trans[good]*abs(arange(azim(good),pa_azimuth(good))-pa_azimuth(good))/90.0)/total(mag.b_trans[good]) signerr = stdev(mag.b_trans[good]*abs(arange(azim(good),pa_azimuth(good))-pa_azimuth(good))/90.0)/sqrt(total(mag.b_trans[good])) endif else begin ;error = total(sqrt((pxi(good)-bxi(good))^2+(pyi(good)-byi(good))^2)) $ ; / nbtrans error = total(sqrt((pxi(good)-bxi(good))^2+$ (pyi(good)-byi(good))^2))/nbtrans sigerr = stdev(sqrt((pxi(good)-bxi(good))^2+$ (pyi(good)-byi(good))^2))/sqrt(nbtrans) nerror = total(sqrt(((pxi(good)-bxi(good))/bxerr(good))^2+$ ((pyi(good)-byi(good))/byerr(good))^2))/nbtrans signerr = stdev(sqrt(((pxi(good)-bxi(good))/bxerr(good))^2+$ ((pyi(good)-byi(good))/byerr(good))^2))/sqrt(nbtrans) endelse avals = [avals,alpha] evals = [evals,error] nevals= [nevals,nerror] sigevals = [sigevals,sigerr] signevals = [signevals,signerr] endif if keyword_set(verbose) then $ message,/info,strcompress('Alpha = '+string(alpha) + $ ', Error = '+string(error) + 'G') alpha = alpha + astep if NOT keyword_set(verbose) AND NOT keyword_set(quiet) then $ progress,100.*(alpha+alimit)/(2.*alimit), $ frequency=5,last=keyword_set(alpha GT alimit) endrep UNTIL (alpha GT alimit) check = check_math(mask=32) ; Ignore underflows !except = except if n_elements(avals) EQ 1 then begin message,/info,'Data cuts were too strict ... no fit' return,0. endif avals = avals(1:n_elements(avals)-1) evals = evals(1:n_elements(evals)-1) nevals = nevals(1:n_elements(nevals)-1) sigevals = sigevals(1:n_elements(sigevals)-1) signevals = signevals(1:n_elements(signevals)-1) emin = min(evals,eplace) if (abs((eplace - n_elements(evals))) le 1 OR eplace lt 1) then begin amin=0 sigalpha=0. print,'too close to edge!' goto, bomb endif ; For the shear, expect 45 deg for random orientations. if keyword_set(shear) AND emin GT 42. then $ message,/info,'Warning the minimum shear is rather large.' amin = line_cent(evals,avals,n_points=3,p_fit=pfit,i_err=lerr) ;if (pfit(0) lt 0) then begin ; amin=0 ; sigalpha=0. ; print,'negative error' ; goto, bomb ;endif if (pfit(0) lt 0 OR pfit(2) lt 0) then begin amin = line_cent(evals,avals,n_points=5,p_fit=pfit,i_err=lerr) if (pfit(0) lt 0 OR pfit(2) lt 0) then begin amin=0 sigalpha=0. print,'negative error or negative curve!' goto, bomb endif endif cc=pfit(0) & bb=pfit(1) & aa=pfit(2) ; reconstruct fit from line_cent ; uncertainty will truly be greater of that for sigevals ; or that for fit... yerr= rebin(sigevals(eplace-1:eplace+1),1) > lerr yminplusdy=cc+bb*amin + aa*amin^2 + yerr ; y-val of line-center position + residuals, ; i.e. y(min)+delta y, so now going to get ; xmin+/-delta x, delta x gives alpha uncertainty. sigalpha=sqrt((yminplusdy-cc+(bb^2/(4*aa)))/aa) if NOT keyword_set(quiet) then $ message,/info,strcompress('Minimum error at alpha = '+ string(amin)) if NOT keyword_set(noplot) then begin win = !d.window ;window,/free if n_elements(plotwindow) LE 0 then begin window,/free plotwindow = !d.window endif else begin device,window=w if plotwindow LT n_elements(w) then begin if w[plotwindow] then wset,plotwindow $ else begin if plotwindow LE 31 then $ window,plotwindow $ else $ window,/free plotwindow = !d.window endelse endif else begin window,/free plotwindow = !d.window endelse endelse if keyword_set(shear) then begin ytitle='Field Weighted Shear Angle (Deg)' title = 'Computed Field-Weighted Shear Angle vs. Alpha' endif else begin ytitle='Error (G/pix)' title = 'Error in Computed Force-Free Transverse Field vs. Alpha' endelse plot,avals,evals,/ynozero,ytitle='Error (G/pix)', $ title=title, $ xtitle=strcompress('Alpha: Best alpha is '+string(amin)) oplot,avals,pfit(0)+pfit(1)*avals+pfit(2)*avals^2,line=2 errplot,avals,pfit(0)-yerr(0)+pfit(1)*avals+pfit(2)*avals^2,$ pfit(0)+yerr(0)+pfit(1)*avals+pfit(2)*avals^2 oplot,intarr(10)+(amin-sigalpha(0)),indgen(10)*500,lines=1 oplot,intarr(10)+amin,indgen(10)*500,lines=2 oplot,intarr(10)+(amin+sigalpha(0)),indgen(10)*500,lines=1 if win GE 0 then wset,win;; else window,/free endif bomb: return,amin end