;+ ;NAME: ; hxtimg_accum ;PURPOSE: ; plots hxt data, prompts for start and end times for data ; accumulation and background, and returns an array with the ; accumulated counts, for each interval and channel. ;CATEGORY: ; Yohkoh HXT ;CALLING SEQUENCE: ; Hxtimg_accum, index, data, index1, hxi_index, counts, sigma, bcounts, $ ; bsigma, img_trange = img_trange, bck_trange=bck_trange, $ ; dt_resolution=dt_resolution, alt_tsel=alt_tsel, $ ; accum_cnts=accum_cnts, ch_accum=ch_accum, cnts_str=cnts_str,$ ; bkgd_str=bkgd_str, hdafile=hdafile, channel=channel, $ ; noback=noback, chisqr=chisqr, st_times = st_times, $ ; int_times = int_times ;INPUT: ; index=data info from hda file ; data=counts in each subcollimator for each channel ;KEYWORDS: (INPUT) ; img_trange = a time range for the images ; bck_range = a time_range for the background accumulation ; dt_resolution = number of 1/2 second (flare mode) or 4 second ; (quiet mode) frames per image, default is to ; prompt the user. If set to -1, then the whole chosen ; interval is used. NOTE THAT THIS KEYWORD IS IGNORED ; IF MULTIPLE INTERVALS HAVE BEEN SELECTED EITHER ; IN PLOT_LCUR, OR IF THEY HAVE BEEN PASSED IN. ; alt_tsel = use alt_tim_sel to obtain interval, rather than plot_lcur ; Useful if you don't have x-windows ; accum_cnts = if set, the value of average counts/SC in a given ; that you want for each set of images, the routine ; will accumulate counts until the number of counts ; in the given channel are greater than or equal to ; this value. If set with a slash, ie. the value=1 ; then the default is 200 (a good value). ; accum_cnts_ch = the channel you want to use for the accum_cnts ; structures = if set, the data will be passed out in structures, ; duplicating the use of MK_HXT_SUM for use in HXTPRO. ; The structure hxi_index need not be passed in if ; the structures keyword is set. THIS IS IMPORTANT: ; DO NOT USE ANY OF THE OPUTPUT VARIABLES IF YOU SET ; THIS KEYWORD, I.E. THE CALLING SEQUENCE ONLY INCLUDES ; INDEX AND DATA AS PARAMETERS, THE REST MUST BE KEYWORDS, ; LIKE THIS: Hxtimg_accum, index, data, /struct, cnts_str=cnts_str, $ ; bkgd_str=bkgd_str, other keywords.... ; hdafile= the name of the hda file used, the structures want this ; i have no clue why, if not set we make one up. ; channel= a channel number for the plots, default is to prompt the user ; noback= if set, the background is set to zero. Requested by HSH 21-jun-95 ; /chisqr= if set, sigma(counts) is a weighted sum of total(sqrt(counts)) ; and sqrt(total(counts)), where the total is taken over the ; accumulation intervals. This may be more useful for cases ; with low signal-to-noise. This used to be the default, ; but not anymore, 6-mar-96, the default is to use sqrt(total(counts)) ; st_times= an array of accumulation start times, in any yohkoh format ; int_times= interval times corresponding to st_times, may be a scalar, ; if this is not set, and st_times is set, the routine will ; use the differences between st_times. This must be an integer ; number of the appropriate resolution in seconds, if the ; value of int_times is less than a frame value. ;KEYWORDS: (OUTPUT) ; cnts_str = the counts in a structure, the structures keyword must ; be set for this to be calculated, for use by HXTPRO ; bkgd_str = the background in a structure, the structures keyword must ; be set for this to be calculated. ;OUTPUT: ; index1= index structure for the first frame of each image interval ; hxi_index.chan= channel used ; =0 L0 (13-23 keV) ; =1 M1 (23-33 keV) ; =2 M2 (33-53 keV) ; =3 HI (53-93 keV) ; hxi_index.actim= image accumulation time in 10ths of secs ; hxi_index.interval= actim ; hxi_index.bkg_time= start time for background accumulation ; (millisec OF day) ; hxi_index.bkg_day= day of start of background accumulation ; (days since 1-Jan-79) ; hxi_index.bkg_actim= background accumulation time in 10ths of secs ; hxi_index.bkg_interval= actim ; counts= the counts accumulated over the chosen intervals ; sigma= sqrt(counts) ; bcounts= the counts accumulated over the background interval ; bsigma= sqrt(bcounts) ;DESCRIPTION: ; Uses the function plot_lcur to select the data to be imaged, ; and the background. ;MODIFICATION HISTORY: ; Written Oct-93 by JMM, ; Changed to give full resolution, 19-jul-94, JMM ; 9/21/94, account for 4.0 second time lag in hda files, jmm ; 1/9/95, changed calculation of sigma to be the same as ; in Tom Metcalf's hxtaccumulate.pro, to be consistent ; in pixon calculations. ; Added the ability to do non-flare mode data, 2/12/95, jmm ; Speed things up using new-style TOTAL commands, 2/13/95, jmm ; 13-mar-95 JPW : img_trange can now be an array of ranges (2,n) ; 16-mar-95 JPW : can now use multiset option in plot_lcur ; Stripped out accumulation stuff, replaced it with calls to ; CHOOSE_INTERVAL, 5-apr-95, jmm ; added cnts_str, bkgd_str keyword outputs, 6-apr-95, jmm ; Added the call to setgra, which gives the subscripts for ; the subcollimator switch, from HXTSBT_MULT, 6-apr-95, jmm. ; This accomodates the structures... ; I don't think anyone will notice. ; Subtract bcakground before choosing the imaging intervals, ; You can see where it's less than zero, and ACCUM_CNTS now ; refers to Excess counts, 9-jun-95 ; Added /noback option, 22-jun-95 ; Fixed bug in calculation of SC uncertainties for single ; interval accumulations, 1-Aug-95, Jmm ; Changed default sigmas back to poisson statistics, added ; /chisqr keyword, 6-mar-96, jmm ; Added st_times, and int_times keywords, 3-jan-1997, jmm ; Call to hxt_chsums is changed, 6-apr-1999, jmm ;- PRO Hxtimg_accum, index0, data, index1, hxi_index, counts, sigma, bcounts, $ bsigma, Img_trange=img_trange, bck_trange=bck_trange, $ dt_resolution=dt_resolution, alt_tsel=alt_tsel, $ accum_cnts=accum_cnts, ch_accum=ch_accum, structures=structures, $ Cnts_str=cnts_str, Bkgd_str=bkgd_str, Hdafile=hdafile, $ channel=channel, noback=noback, chisqr=chisqr, $ st_times=st_times, int_times=int_times ;The structure hxi_index may not be passed in, if the structures keyword is set IF(N_ELEMENTS(hxi_index) GT 0) THEN qhxi_index = 1 ELSE qhxi_index = 0 IF(qhxi_index) THEN serr = hxi_index(0).syserr ELSE serr = 0.0 IF(KEYWORD_SET(chisqr)) THEN qchisqr = 1 ELSE qchisqr = 0 IF(KEYWORD_SET(structures) AND n_params() GT 2) THEN BEGIN message, /info, 'THE STRUCTURES KEYWORD HAS BEEN SET' message, /Info, 'DO NOT USE ANY NON-KEYWORDS IN THE FUNCTION CALL OTHER THAN INDEX AND DATA' RETURN ENDIF ;sum SC's to plot and you'll need a time array for each data point ;Also deal with any Cal mode data. zz = hxt_chsums(index0, data, /half_sec, tim_arr = tim_arr1, $ tim_sec = tim_arr, dt = dx4, index_out = index, $ data_out = data_x) time0 = tim_arr1(0) nc0 = 4 ;no. of channels nc01 = nc0-1 ndset = N_ELEMENTS(index) ;no. of data sets nds1 = ndset-1 ;For full resolution, 4 minor frames for each data set nds4 = ndset*4 & nds41 = nds4-1 ;reform data to get into the proper form data_x = reform(data_x, nc0, 64, nds4, /overwrite) ;If i need to choose intervals by plotting no_time_selected = NOT (KEYWORD_SET(img_trange) OR KEYWORD_SET(st_times)) no_bck_selected = NOT (KEYWORD_SET(bck_trange) OR KEYWORD_SET(noback)) IF(no_time_selected OR no_bck_selected) THEN BEGIN IF(N_ELEMENTS(channel) GT 0) THEN ich = fix(channel(0)) ELSE BEGIN print, ' INPUT CHANNEL NUMBER FOR THE INTERVAL>>' print, ' 0 FOR LO CHANNEL' print, ' 1 FOR M1 CHANNEL' print, ' 2 FOR M2 CHANNEL' print, ' 3 FOR HI CHANNEL' read, ich ich = fix(ich) ENDELSE ENDIF ELSE ich = 0 y0 = reform(zz(ich, *))/(dx4*64.0) ;y0 is the curve used to find the intervals IF(KEYWORD_SET(noback)) THEN BEGIN print, 'BACKGROUND SET TO ZERO, NOT RECOMMENDED' bcounts = fltarr(64, 4) bsigma = bcounts ib1 = 0 & ib2 = 3 ; i may need this stuff, for the structures ib14 = ib1/4 & ib24 = ib2/4 bactim = total(dx4(ib1:ib2)) ENDIF ELSE BEGIN ;choose the background interval print, 'BACKGROUND SELECTION::' choose_interval, tim_arr1, ib1, ib2, data = y0, plot_title = 'HXT COUNTS/SECOND/SC', $ time_range = bck_trange, dt_resolution = -1, alt_tsel = alt_tsel ib1 = ib1(0) & ib2 = ib2(0) ; Check time interval for non-flare mode data ib14 = ib1/4 & ib24 = ib2/4 not_flare = where(gt_dp_mode(index(ib14:ib24)) NE 9) IF(not_flare(0) NE -1) THEN BEGIN message, /info, 'WARNING, NON-FLARE MODE BACKGROUND CHOSEN, IS ZERO FOR HI, M1, AND M2 CHANNELS' ENDIF ; time interval, bactim = total(dx4(ib1:ib2)) ; Counts array tcounts = data_x(*, *, ib1:ib2) ;Speed things up, using total commands on 1 dimension nibs = ib2-ib1+1 IF(nibs GT 1) THEN BEGIN ; new stuff, see explanation below ; bsigma = sqrt(bcounts) is not true anymore ; error1 is total(sqrt(counts)), bsigma is sqrt(total(counts)) bcounts = transpose(total(tcounts, 3)) bsigma = transpose(sqrt(total(tcounts, 3)+ $ total((serr*tcounts)^2, 3))) IF(qchisqr) THEN BEGIN ; jmm, 6-mar-96 error1 = transpose(sqrt(total(sqrt(tcounts), 3)^2 $ +total((serr*tcounts)^2, 3))) mn = avsig_1(tcounts, sdev, dim = 3) ;mean and st deviation mn = mn > 0 xxx = where(sdev GT 0.0) mn(xxx) = ((1.0-sqrt(mn(xxx))/sdev(xxx)) > 0.0) < 1.0 wt = transpose(mn) bsigma = error1*wt + bsigma*(1.0-wt) ENDIF ELSE bsigma = bsigma ENDIF ELSE BEGIN ;no need to total bcounts = transpose(tcounts) bsigma = transpose(sqrt(tcounts+(serr*tcounts)^2)) ;oops, jmm,1-aug-95 ENDELSE ENDELSE ; Ok, now you choose the image interval, subtract the background first ; you need a sum for the background IF(KEYWORD_SET(noback)) THEN bcps = fltarr(4) ELSE BEGIN IF(ib2 NE ib1) THEN bcps = total(zz(*, ib1:ib2), 2) ELSE bcps = zz(*, ib1) ENDELSE bcps = bcps/(64.0*bactim) ;average background rate per SC print, 'IMAGE INTERVAL SELECTION::' IF(KEYWORD_SET(accum_cnts)) THEN BEGIN IF(accum_cnts EQ 1) THEN acc_cnts = 200.0 ELSE acc_cnts = accum_cnts(0) IF(N_ELEMENTS(ch_accum) GT 0) THEN ich0 = ch_accum(0) ELSE ich0 = 0 y0 = reform(zz(ich0, *)/64.0) y0 = y0-bcps(ich0)*dx4 choose_interval, tim_arr1, ss_st, ss_en, data = y0, accum_cnts = acc_cnts, $ time_range = img_trange, plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SC', $ st_times = st_times, int_times = int_times ENDIF ELSE BEGIN y0 = y0-bcps(ich) choose_interval, tim_arr1, ss_st, ss_en, data = y0, $ plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SEC/SC', $ time_range = img_trange, dt_res = dt_resolution, alt_tsel = alt_tsel, $ st_times = st_times, int_times = int_times ENDELSE ; Check time interval for non-flare mode data ix14 = min(ss_st)/4 & ix24 = max(ss_en)/4 not_flare = where(gt_dp_mode(index(ix14:ix24)) NE 9) IF(not_flare(0) NE -1) THEN BEGIN message, /info, 'WARNING, NON-FLARE MODE DATA IN CHOSEN INTERVAL IS ZERO FOR HI, M1, AND M2 CHANNELS' ENDIF n_intv = N_ELEMENTS(ss_st) ;You'll not have an index structure for each minor frame, so for ;each minor frame, find the appropriate index, and reset .time ;and .day ss_st_index = ss_st/4 & index1 = index(ss_st_index) index1.gen.time = tim_arr1(ss_st).time index1.gen.day = tim_arr1(ss_en).day ; now accumulate data, dti will be the time interval length, in seconds counts = fltarr(64, nc0, n_intv) IF(qchisqr) THEN BEGIN ;newer stuff, jmm, 6-mar-96 wt = counts ;new stuff, 1/9/95, jmm error1 = counts ENDIF sigma = counts dti = fltarr(n_intv) FOR j = 0, n_intv-1 DO BEGIN dti(j) = total(dx4(ss_st(j):ss_en(j))) tcounts = data_x(*, *, ss_st(j):ss_en(j)) dnn = ss_en(j)-ss_st(j)+1 IF(dnn GT 1) THEN BEGIN counts(*, *, j) = transpose(total(tcounts, 3)) sigma(*, *, j) = transpose(sqrt(total(tcounts, 3) $ +total((serr*tcounts)^2, 3))) IF(qchisqr) THEN BEGIN ;newer stuff, jmm, 6-mar-96 error1(*, *, j) = transpose(sqrt(total(sqrt(tcounts), 3)^2 $ +total((serr*tcounts)^2, 3))) mn = avsig_1(tcounts, sdev, dim = 3) ;mean and st dev mn = mn > 0 xxx = where(sdev GT 0.0) mn(xxx) = ((1.0-sqrt(mn(xxx))/sdev(xxx)) > 0.0) < 1.0 wt(*, *, j) = transpose(mn) ENDIF ELSE BEGIN error1 = 0.0 wt = 0.0 ENDELSE ENDIF ELSE BEGIN ;no need to total counts(*, *, j) = transpose(tcounts) sigma(*, *, j) = transpose(sqrt(tcounts+(serr*tcounts)^2)) error1 = 0.0 wt = 0.0 ENDELSE ENDFOR ; now calculate the uncertainty=sqrt(counts) ; The error would be total(sqrt(counts)) if each observation ; were independent. The error would be sqrt(total(counts)) ; if each observation were the same. Reality is in between. ; Use the standard deviation compared to the error on a single ; observation to weight each type of error. ; error1 is total(sqrt(counts)), sigma is sqrt(total(counts)) sigma = error1*wt + sigma*(1.0-wt) ; Replicate and fill hxi_index if necessary, IF(qhxi_index) THEN BEGIN hxi_index.bkg_time = tim_arr1(ib1).time hxi_index.bkg_day = gt_day(tim_arr1(ib1)) hxi_index.bkg_actim = fix(10.0*bactim) hxi_index.bkg_interval = fix(10.0*(tim_arr(ib2)-tim_arr(ib1)+dx4(ib2))) hxi_index = replicate(hxi_index(0), nc0, n_intv) FOR i = 0, nc01 DO hxi_index(i, *).chan = i FOR j = 0, n_intv-1 DO BEGIN hxi_index(*, j).actim = fix(10.0*dti(j)) hxi_index(*, j).interval = fix(10.0*(tim_arr(ss_en(j))-tim_arr(ss_st(j))+dx4(ss_en(j)))) ENDFOR ENDIF data_x = 0.0 ; hxt_setgra gives the subscripts for the subcollimator switch hxt_setgra, sn2gra, gra2sn bcounts = bcounts(gra2sn, *) counts = counts(gra2sn, *, *) bsigma = bsigma(gra2sn, *) sigma = sigma(gra2sn, *, *) ;For structure, this should be moderately easy IF(KEYWORD_SET(structures)) THEN BEGIN ss_en_index = ss_en/4 set_hxt_struct, cnts_str ;set counts structure cnts_str.gen.id = ' HXTIMG_ACCUM Output' ;and start filling cnts_str.inf_raw.comment = 'No Comment' IF(KEYWORD_SET(hdafile)) THEN cnts_str.inf_raw.rf_file = hdafile $ ELSE cnts_str.inf_raw.rf_file = 'hda'+ex2fid(anytim2ex(index(0))) cnts_str = replicate(cnts_str(0), n_intv) ;the times look formatted. This looks like a mess.... t7 = anytim2ex(tim_arr1(ss_st)) t7_str = strcompress(string(t7), /remove_all) lt_10 = where(t7 LT 10) IF(lt_10(0) NE -1) THEN t7_str(lt_10) = ' '+t7_str(lt_10) lt_100 = where(t7(3,*) LT 100) IF(lt_100(0) NE -1) THEN t7_str(3, lt_100) = ' '+t7_str(3, lt_100) t71 = anytim2ex(anytim2ints(tim_arr1(ss_st), offset=dti)) t71_str = strcompress(string(t71), /remove_all) lt_10 = where(t71 LT 10) IF(lt_10(0) NE -1) THEN t71_str(lt_10) = ' '+t71_str(lt_10) lt_100 = where(t71(3,*) LT 100) IF(lt_100(0) NE -1) THEN t71_str(3, lt_100) = ' '+t71_str(3, lt_100) time_str = t7_str(6, *)+'/'+t7_str(5, *)+'/'+t7_str(4, *) $ +' '+t7_str(0, *)+':'+t7_str(1, *)+':'+t7_str(2, *)+'.'+t7_str(3, *) $ +' - '+t71_str(0, *)+':'+t71_str(1, *)+':'+t71_str(2, *)+'.'+t71_str(3, *) FOR j = 0, n_intv-1 DO BEGIN cnts_str(j).inf_raw.range(0, *) = ss_st_index(j) ;the subscript in index for the start cnts_str(j).inf_raw.range(1, *) = ss_st(j)-4*ss_st_index(j) ;start minor frame cnts_str(j).inf_raw.range(2, *) = ss_en_index(j) ;the subscript in index for the end cnts_str(j).inf_raw.range(3, *) = ss_en(j)-4*ss_en_index(j) ;end minor frame cnts_str(j).inf_raw.time(*) = time_str(j) cnts_str(j).data.acc_tim(*) = dti(j) cnts_str(j).data.cts = counts(*, *, j) cnts_str(j).data.sigma = sqrt(counts(*, *, j)) ENDFOR set_hxt_struct, bkgd_str ;set counts structure bkgd_str.gen.id = ' HXTIMG_ACCUM Output' ;and start filling bkgd_str.inf_raw.comment = 'No Comment' IF(KEYWORD_SET(hdafile)) THEN bkgd_str.inf_raw.rf_file = hdafile $ ELSE bkgd_str.inf_raw.rf_file = 'hda'+ex2fid(anytim2ex(index(0))) ;the times look formatted. t7 = anytim2ex(tim_arr1(ib1)) t7_str = strcompress(string(t7), /remove_all) lt_10 = where(t7 LT 10) IF(lt_10(0) NE -1) THEN t7_str(lt_10) = ' '+t7_str(lt_10) IF(t7(3) LT 100) THEN t7_str(3) = ' '+t7_str(3) t71 = anytim2ex(anytim2ints(tim_arr1(ib1), offset = bactim)) t71_str = strcompress(string(t71), /remove_all) lt_10 = where(t71 LT 10) IF(lt_10(0) NE -1) THEN t71_str(lt_10) = ' '+t71_str(lt_10) IF(t71(3) LT 100) THEN t71_str(3) = ' '+t71_str(3) bkgd_str.inf_raw.time(*) = t7_str(6)+'/'+t7_str(5)+'/'+t7_str(4) $ +' '+t7_str(0)+':'+t7_str(1)+':'+t7_str(2)+'.'+t7_str(3) $ +' - '+t71_str(0)+':'+t71_str(1)+':'+t71_str(2)+'.'+t71_str(3) bkgd_str.inf_raw.range(0, *) = ib14 ;the subscript in index for the start bkgd_str.inf_raw.range(1, *) = ib1-4*ib14 ;start minor frame bkgd_str.inf_raw.range(2, *) = ib24 ;the subscript in index for the end bkgd_str.inf_raw.range(3, *) = ib2-4*ib24 ;end minor frame bkgd_str.data.acc_tim(*) = bactim bkgd_str.data.cts = bcounts bkgd_str.data.sigma = sqrt(bcounts(*, *)) ENDIF RETURN END