;+ ; NAME: ; atr2hxa_dbase ; PURPOSE: ; generate database that contains the solar radii as seen by the ; HXA detector. database is used for the HXA limb reconstruction ; in hxa2hxaxy. ; CATEGORY: ; CALLING SEQUENCE: ; atr2hxa_dbase,weeks,years ; INPUTS: ; weeks = list of week numbers to be read ; years = year of the week(s) to b read. Either scalar or vector ; of same size as weeks ; KEYWORDS (INPUT): ; dir = directory for ATR files (default is DIR_GEN_ATR) ; dbname = name of the database file to be updated (or created) ; OUTPUTS: ; KEYWORDS (OUTPUT): ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; Calculates average values for each week. ; To catch a fast transition, such as a HXA gain change, ; a search for up to one significant transition within the ; week is done. ; MODIFICATION HISTORY: ; JPW, 17-sep-96 ; JPW, 12-sep-98 changed dfmax from 0.05 to 0.25 ;- pro atr2hxa_dbase,weeks,years,dir=dir,dbname=dbnam,verbose=verb t2nitmin = 1.0 ; only consider values at least xxx min. before sunset dfmax = 0.25 ; maximum deviation of fx,fy from 1.0 to be considered n_min = 1001 ; minimum number of datapoints to be considered for database mdif = 0.002 ; minimum difference in fx,fy to trigger a "signif. change" if n_elements(dbnam) eq 0 then $ dbnam = concat_dir('DIR_GEN_STATUS','hxa_dbase.genx') ;if n_elements(dbnam) eq 0 then dbnam = '$ys/ucon/idl/wuelser/hxa_dbase.genx' dbn = findfile(dbnam) if dbn(0) ne '' then begin restgen,db,file=dbnam spawn,'cp '+dbnam+' '+dbnam+'.BAK' endif if n_elements(years) ne n_elements(weeks) then $ years=intarr(n_elements(weeks))+years(0) hxa_parms, hpix=hpix, h_xx=h_xx, h_yy=h_yy for i=0,n_elements(weeks)-1 do begin rd_atr,weeks(i),years(i),atr,indir=dir hxa = (atr.hxa > 0) < 2047 ; create bit pattern of hidden limbs. zon = intarr(n_elements(atr)) forzones,xzon,yzon whe = where(xzon(hxa(0,*))) if whe(0) ge 0 then zon(whe) = zon(whe) or 1 whe = where(xzon(hxa(1,*))) if whe(0) ge 0 then zon(whe) = zon(whe) or 2 whe = where(yzon(hxa(2,*))) if whe(0) ge 0 then zon(whe) = zon(whe) or 4 whe = where(yzon(hxa(3,*))) if whe(0) ge 0 then zon(whe) = zon(whe) or 8 xoff = (hxa(0,*)+hxa(1,*))/2.0d - h_xx yoff = (hxa(2,*)+hxa(3,*))/2.0d - h_yy xrad = (hxa(1,*)-hxa(0,*))/2.0d yrad = (hxa(3,*)-hxa(2,*))/2.0d sunr = get_rb0p(atr,/radius)/hpix ; solar radius in hxa pixels fx = sqrt(yoff*yoff+xrad*xrad) / sunr fy = sqrt(xoff*xoff+yrad*yrad) / sunr tim2orbit,atr,tim2night=t2nit,/simple ;offmax = 90.0 ;offmax = 150.0 offmax = 120.0 whe = where(zon eq 0 and abs(xoff) lt offmax and abs(yoff) lt offmax $ and abs(1.0-fx) lt dfmax and abs(1.0-fy) lt dfmax $ and t2nit gt t2nitmin,nn) ;print,offmax,n_elements(whe) ; for testing if nn ge n_min then begin tim = anytim2ints(atr(whe)) fx = fx(whe) fy = fy(whe) ; calc. average afx = total(fx)/nn afy = total(fy)/nn atim = anytim2ints(tim(0), $ off=total(int2secarr(tim,tim(0)))/nn) ; search for (up to one) significant change ; coupling of x and y assumes that changes in x and y are either ; correlated or uncorrelated but not anticorrelated. This should ; be true for either degradation or proper HXA adjustment. if nn ge 2*n_min then begin fxy = smooth(fx,n_min) + smooth(fy,n_min) fxy = fxy(n_min/2:nn-n_min/2-1) fmi = min(fxy,max=fma) ; "clouds" of low and high values: wl = where(fxy lt fmi+0.3*(fma-fmi),nwl) + n_min/2 ; re-add offset wh = where(fxy gt fma-0.3*(fma-fmi),nwh) + n_min/2 ; of fxy:n_min/2 ; redefine clouds in terms of "early" and "late" cloud ; also reverse early cloud to have element 0 at the end if total(wl)/nwl le total(wh)/nwh then begin wmi = rotate(wl,2) nwmi = nwl wma = wh nwma = nwh endif else begin wmi = rotate(wh,2) nwmi = nwh wma = wl nwma = nwl endelse ;!p.multi = [0,1,3] ;utplot,tim,fx,psym=3,charsiz=1.8,/ynoz ;utplot,tim,fy,psym=3,charsiz=1.8,/ynoz ;utplot,tim,fxy,psym=3,charsiz=1.8,/ynoz ;outplot,tim(wmi),fxy(wmi-n_min/2),psym=1 ;outplot,tim(wma),fxy(wma-n_min/2),psym=4 ; find "edges" of clouds, or best guess of separation point if ; they overlap -> wmi0, wma0 nwia = nwmi < nwma dt = abs(int2secarr(tim(wmi(0:nwia-1)),tim(wma(0:nwia-1)))) wdt = where(dt eq min(dt)) wmi0 = wmi(wdt(0)) wma0 = wma(wdt(0)) ; redefine clouds again, everything before/after separation point: wmi = lindgen(wmi0+1) nwmi = wmi0+1 wma = lindgen(nn-wma0)+wma0 nwma = nn-wma0 ; calculate separate values before/after change: afx2 = [total(fx(wmi))/nwmi,total(fx(wma))/nwma] afy2 = [total(fy(wmi))/nwmi,total(fy(wma))/nwma] ; atim2 = anytim2ints(tim(0), $ ; off=[total(int2secarr(tim(wmi),tim(0)))/nwmi, $ ; total(int2secarr(tim(wma),tim(0)))/nwmi]) atim2 = tim([wmi0,wma0]) ;outplot,[atim2(0),atim2(0),atim2(1),atim2(1)],[0.9,1.1,1.1,0.9] ; apply a set of conditions to see if we really think there ; is a significant change if (abs(afx2(0)-afx2(1)) gt mdif or abs(afy2(0)-afy2(1)) gt mdif) $ and nwmi gt n_min and nwma gt n_min then begin afx = afx2 afy = afy2 atim = atim2 endif endif ;print,weeks(i),afx,afy,' ',fmt_tim(atim) ; now update the database: na = n_elements(afx) if n_elements(db) eq 0 then begin db = replicate({hxa_rad_db,time:0L,day:0,hxa_fx:10000,hxa_fy:10000},na) endif else begin ss = sel_timrange(db,atr(0),atr(n_elements(atr)-1),/between,/boolean) ww = where(ss eq 0) if ww(0) ge 0 then db = [replicate(db(0),na),db(ww)] $ else db = replicate(db(0),na) endelse ii = indgen(na) db(ii).day = atim.day db(ii).time = atim.time db(ii).hxa_fx = (afx+0.00005)*10000 db(ii).hxa_fy = (afy+0.00005)*10000 endif ; sort the database according to time and rewrite it if n_elements(db) gt 1 then db = db(sort(int2secarr(db))) if n_elements(db) ge 1 then savegen,db,file=dbnam if keyword_set(verb) then $ print,weeks(i),years(i),n_elements(atr),n_elements(db) endfor end