;+ ; Main program for SPEX - Spectral Analysis Executive ; ; mod, ras, 12-apr-1997, change call from PS to SPS for ; graphics. ; Paul Bilodeau, 3-April-2002, changed calls to wsetshow to use ; TIMEHISTORY and SPECTRAl keywords. ;- @function_com.pro @spex_commons.pro @winup_common.pro ; common window_array, windows_index, windows_title common op_com_com, com_stack ;command stack for op_com initialize: setup_spex, com_stack, style_sel, need, prompt, comments, $ script_opt1, script_opt2, standard_types, bdata, $ dir_queue, edg_units, last_found, last_read, wtitles, gaction, $ gdaction, goaction, plotbuffer, bang_enter, setup_sys, $ nowidget=nowidget, merge_list=merge_list, last_plot=last_plot, $ handles=handles setup_params, prm_spex, prm_spex_info not_first_time: parameters = prm_spex.nam prm_spex_sav = prm_spex.val for i = 0, n_elements(prm_spex) - 1 do begin nv = prm_spex(i).numval if nv gt 1 then temp = fltarr(nv) else temp = 0. r = execute(parameters(i) + '=' + prm_spex(i).typ + '(temp)') r = execute('reads,prm_spex(i).val,' + parameters(i)) endfor families = 'main' ; family of parameters to display in opcom opcom: alpha_page ;Temporary patch - some variables have been changed by the action of the ;procedure, such as uts,ute,apar,energy_bands,erange, set their ;corresponding parameters ;However, if they have been changed, then do nothing. ; for i = 0, n_elements(prm_spex) - 1 do begin temp = where( prm_spex(i).val ne prm_spex_sav(*,i), r) if r eq 0 then r = execute('prm_spex(i).val=' + parameters(i)) endfor ;for i = 0, n_elements(prm_spex) - 1 do then $ ; r = execute('prm_spex(i).val=' + parameters(i)) ;Enter here if parameter structure doesn't need updating opcom_fast: ;print,'opcom_fast: input_line = ',fcheck(input_line,'') ;more,['COM_STACK',com_stack] op_com, prm_spex, prm_spex_info, input_line=input_line, families=families,$ prompt=prompt, noreadinput=fcheck(noreadinput,0), $ com_sel=command, com_param=com_param, com_line=com_line, $ nolist=nolist, which=ipar,error=error,com_del='!!' nolist = 1 if error then goto, opcom ;LOOK FOR A CHANGE IN MODEL FUNCTION AND SET UP DEFAULTS i_f_model = (where(prm_spex.nam eq 'F_MODEL'))(0) if prm_spex(i_f_model).val(0) ne prm_spex_sav(0,i_f_model) then begin ;the model has changed, make other changes temp = model_components(prm_spex(i_f_model).val(0), number_param, param_ranges, free, apar) prm_spex((where(prm_spex.nam eq 'FREE'))(0)).numval = number_param prm_spex((where(prm_spex.nam eq 'FREE'))(0)).val = free prm_spex((where(prm_spex.nam eq 'APAR'))(0)).numval = number_param prm_spex((where(prm_spex.nam eq 'APAR'))(0)).val = apar prm_spex((where(prm_spex.nam eq 'RANGE_LO'))(0)).numval = number_param prm_spex((where(prm_spex.nam eq 'RANGE_HI'))(0)).numval = number_param prm_spex((where(prm_spex.nam eq 'RANGE_LO'))(0)).val = param_ranges(*,0) prm_spex((where(prm_spex.nam eq 'RANGE_HI'))(0)).val = param_ranges(*,1) endif ; Move new values of parameters into variables that have the same names ; as the parameter name. parameters = prm_spex.nam result = (where(prm_spex.val ne prm_spex_sav,nv) / n_elements(prm_spex_sav(*,0))) if nv ge 1 then begin; some value has changed, must change corresponding parameter result= result(uniq(result,sort(result))) for i = 0, n_elements(result) - 1 do begin j = result(i) ;change these parameters nv = prm_spex(j).numval if nv gt 1 then temp = fltarr(nv) else temp = 0. r = execute(parameters(j) + '=' + prm_spex(j).typ + '(temp)') r = execute('reads,prm_spex(j).val,' + parameters(j)) endfor endif prm_spex_sav = prm_spex.val if n_elements(edges) gt 2 and need.wuse then $ wuse = where( edges(1,*) gt erange(0) and edges(0,*) lt erange(1)) range = reform( [range_lo(*),range_hi(*)], n_elements(range_lo), 2) processor: if ipar ne -1 then goto, opcom_fast ;ipar is -1 if it is a command wgraph = where( command(0) eq gaction, nwgraph) if nwgraph eq 1 then plotbuffer(0) = [gdaction(wgraph),plotbuffer(0)] ;temporary fix for plotbuffer to take care of 1 specific case, ras, 2-jun-95 if command(0) eq 'photon_spectrum' and strupcase(com_param(0)) eq 'CAL_TEST' then $ plotbuffer(0)='photon_spectrum,cal_test' case 1 of command(0) eq '/': families='main' command(0) eq '/time_history': families='time' command(0) eq '/interval': families='inte' command(0) eq '/drm': families='drm' command(0) eq '/fitting': families='fit' command(0) eq '/summary': families='summary' command(0) eq 'multidrm': goto, drm command(0) eq 'accum': goto, accum command(0) eq 'merge_batse': goto, merge_batse command(0) eq 'time_history': goto, time_history command(0) eq 'count_spectrum' or command(0) eq 'photon_spectrum': goto, set_spectrum command(0) eq 'list': nolist=0 command(0) eq 'exit': goto, done command(0) eq 'fitting': goto, fitting command(0) eq 'execute' or command(0) eq 'idl': goto, xecute command(0) eq 'stop': goto, stop_pause command(0) eq 'zoom': goto, zoom command(0) eq 'display_intervals': goto, display_intrvl command(0) eq 'select_interval': goto, select_intrvl command(0) eq 'background_interval': goto, background command(0) eq 'initialize': goto, initialize command(0) eq 'read_drm': goto, drm command(0) eq 'build_drm': goto, drm command(0) eq 'recalibrate': goto, drm command(0) eq 'cal_save': goto, drm command(0) eq 'cal_restore': goto, drm command(0) eq 'graph': goto, graph command(0) eq 'steps_help': goto, steps_help command(0) eq 'save_event': goto, save_event command(0) eq 'restore_event': goto, restore_event command(0) eq 'erase': if !d.name eq 'TEK' then erase command(0) eq 'tek': set_plot,'tek' command(0) eq 'create_ps': goto, createps command(0) eq 'pspreview': goto, pspreview command(0) eq 'hard': goto, printhard command(0) eq 'find_files': goto, find_datafiles command(0) eq 'ch_bands': goto, ch_bands command(0) eq 'spec_bands': goto, spec_bands command(0) eq 'nospec_bands': iegroup=0 command(0) eq 'bspectrum': goto, bspectrum ;plot background spectrum command(0) eq 'restore_dflts': goto, restore_dflts ;restore, dflts_file command(0) eq 'save_dflts': goto, restore_dflts ;save_dflts ; command(0) eq 'enrange': goto, erange ;fitting energy range command(0) eq 'clear_ut': goto, clear_ut command(0) eq 'kleanplot': cleanplot command(0) eq 'force_apar': goto, force_apar command(0) eq 'comments': goto, comments command(0) eq 't_hist_mode': goto, t_hist_mode command(0) eq '?': goto, spex_help command(0) eq 'script': goto, script command(0) eq 'preview': goto, preview command(0) eq 'online?': goto, online command(0) eq 'batse_menu': goto, batse_menu command(0) eq 'reset_event': begin ;if last_read.id is an impossible value, then time_history must read event file! last_read.id = -1 last_found.id= -1 end else: endcase if strpos(command(0),'/') eq 0 then nolist=0 goto, opcom ;#------------------------------------------------------------ comments: for i=0,n_elements(comments)-1 do printx, comments(i) printx,' ' goto, opcom_fast ;#------------------------------------------------------------ online: printx,'The database has these events corresponding to DIR_QUEUE and DATA_TIPE:' if strupcase(data_tipe(0)) eq 'BATSE' then begin if n_elements(datalist) eq 0 then begin datalist = [''] ;for i=1,n_elements(dir_queue)-1 do $ ; datalist=[datalist,findfile(concat_dir(dir_queue(i),'sherb*.*'))] datalist = [datalist, loc_file( path=dir_queue, /all, 'sherb*.*')] break_file, datalist(1:*), disk_list, dir_list datalist = lonarr(n_elements(dir_list)) for i=0,n_elements(datalist)-1 do $ datalist(i) = call_function( $ "burst_flare",/burst,long(strmid(dir_list(i),strpos(strupcase(dir_list(i)),'DB_')+3,5))) endif printx, 'These are the FLARE numbers for the BATSE events which are online:' printx, datalist endif else printx,'Database searches only valid for BATSE now.' goto, opcom ;#------------------------------------------------------------ spex_help: if com_param(0) ne '' then begin help_spex, prm_spex, list=com_param, outcommand=outcommand if outcommand ne '' then help, name=outcommand+'*' endif else begin line = '' printx,'For help, Enter the name(s) of the parameter, variable, or command '+ $ 'with comma delimiters' printx,'Or for a list of all parameters, variables, and commands .' read, 'Enter name(s) or : ',line if line ne '' then help_spex,prm_spex, list = line else begin printx,'Enter ?,item for more information at SPEX>.' help_spex, /short endelse endelse goto, opcom_fast ;#------------------------------------------------------------ restore_dflts: ;determine file name if strupcase(dflts_file) eq 'AUTO' then $ dflts_file = strlowcase(arr2str(data_tipe,delim='_')+'_'+strtrim(flare,2)+'_'+strtrim(det_id,2)) break_file, dflts_file,disk_log, dir, filnam, ext if ext eq '' then dflts_file = dflts_file + '.dflt' if command(0) eq 'save_dflts' then goto, save_dflts ; temp = loc_file( path=([curdir(),sumdirectory])[keyword_set(sumdirectory)], $ strupcase(dflts_file), count=r) if r eq 0 then begin temp = loc_file(path=([curdir(),sumdirectory])[keyword_set(sumdirectory)], $ strlowcase(dflts_file), count=r) if r eq 0 then $ printx, 'Cannot find defaults file, '+ $ concat_dir(chklog(sumdirectory), strlowcase(dflts_file)) endif if r ge 1 then begin restore, temp(0) setut, utstart=uts, utend=ute endif goto, opcom ;#------------------------------------------------------------ save_dflts: ;result = execute('save, file="'+concat_dir(chklog(sumdirectory), strlowcase(dflts_file))+ $ ; '",'+arr2str([parameters]) ) extension ='.dat' temp = form_filename(strlowcase(dflts_file), extension,directory = chklog(sumdirectory)) result = execute('save, file="'+ temp + $ '",'+arr2str([parameters]) ) printx, 'Writing defaults file, ' + temp goto, opcom ;#------------------------------------------------------------ xecute: line = '' if com_param(0) eq '' then begin printx,'Enter an IDL statement to execute.' read,prompt='=> ',line endif else line = arr2str(com_param,delim=',') if strpos(strupcase(strtrim(line,2)), 'EXIT') ne -1 then begin answer = 'NO' read, 'Do you really want to exit IDL? (Def = No)', answer if strtrim(answer,2) eq '' then answer='NO' if strpos(strupcase(strtrim(answer,2)), 'Y') ne 0 then goto, opcom endif printx, !prompt + line result = execute(line) goto, opcom ;#------------------------------------------------------------ stop_pause: printx, 'Type .c to continue' stop goto, opcom ;#------------------------------------------------------------ t_hist_mode: printx, 'Enter mode for time history plots,' printx, ' COUNTS "Counts/s" or FLUX "Counts/s/keV/cm2"' line = '' read, 'Mode = ', line if strmid(strupcase(line),0,1) eq 'C' then begin t_hist_mode='COUNTS' endif else begin t_hist_mode='FLUX' endelse goto, opcom ;#------------------------------------------------------------ zoom: if need.tplot then begin printx,'No plot to zoom on.' case last_plot(0) of 'graph': goto,graph 'time_history': goto,time_history else: goto,graph endcase endif if !d.name eq xdevice() then wsetshow, /TIMEHISTORY getut, uts=uts, ute=ute zoom_limits = [uts,ute] - getutbase() spex_intervals, ut, zoom_limits, xinput=zoom_limits, $ style='b', v_styles=['b'], query=check_defaults,color=fcolor(7), $ xoffset= (indgen(3)*200+200),yoffset = (indgen(3)*200+200), $ message='SELECT XRANGE FOR ZOOM', error=error, input_variable='zoom_limits' zoom_limits = limits(zoom_limits) uts = atime(zoom_limits(0)+getutbase()) ute = atime(zoom_limits(1)+getutbase()) setut, uts=uts, ute=ute input_line = 'graph' goto, opcom ;#------------------------------------------------------------ batse_menu: getut,uts=old_uts, ute=old_ute, utb=old_utb,/string if !d.name eq xdevice() then window,2,title='BATSE ARCHIVE' bangx = !x if com_param(0) ne '' then call_procedure, "batse_menu", com_param(0) else $ call_procedure, "batse_menu" !x = bangx setut,uts=old_uts, ute=old_ute, utb=old_utb goto, opcom ;#------------------------------------------------------------ clear_ut: ;reset the utstart and end to the limits of the dataset setut, uts=min(ut + getutbase(0)), ute=max(ut+getutbase()) getut, uts=uts, ute=ute,/s uts_test = uts & ute_test = ute goto, opcom ;#------------------------------------------------------------ merge_batse: ;If the command "merge_batse, current" is given, then the ;values for the detectors and/or files in the merge_list structure are ;used to generate the merge script, otherwise they are generated from ;the optimum detector ids ; ; "merge_batse, norun" means the script is generated but the ; it isn't sent to the command processor ; spex_merge_control, com_param, merge_list, mscript, input_line goto, opcom ;#------------------------------------------------------------ script: ;using the command buffer capability of op_com, execute a series ;of SPEX commands if strupcase(strmid(com_param(0),0,1)) eq 'C' then $ ;command parameter is createps script_commands = script_opt2 else script_commands = script_opt1 printx,['This is the series of commands to be executed consecutively:',script_commands] input_line = arr2str(script_commands,delim = '!!') goto, opcom ;#------------------------------------------------------------ createps: ;repeat the plotbuffer to the postscript file wgraph = where( plotbuffer(0) eq goaction, nwgraph) if nwgraph eq 1 then $ redo_plot = [plotbuffer(1),plotbuffer(0)] $;last command was overplot else redo_plot = plotbuffer(0) pbuffer = strlowcase(parse_comline(redo_plot(0))) timecode = (strmid(strtrim(anytim(!stime,/ex)+100,2),1,2))([6,5,4,0,1,2]) psfile = pbuffer(0)+'_'+ arr2str(timecode(0:2),delim='') + $ '_'+arr2str(timecode(3:5),delim='') case psstyle of 'GIF': begin if com_param(0) eq '' then psfile = psfile +'.gif' else psfile=com_param(0) tvlct, red, green, blue, /get write_gif, psfile, tvrd(), red, green, blue end else: begin if com_param(0) eq '' then psfile = psfile +'.ps' else psfile=com_param(0) redo_plot = ['idl,sps,/'+psstyle,'idl,device,file=psfile', $ redo_plot,'idl,device,/close & x'] printx,'This is the series of commands to be executed consecutively:' for i=0,n_elements(redo_plot)-1 do printx, redo_plot(i) input_line = arr2str(redo_plot,delim = '!!') end endcase goto, opcom ;#------------------------------------------------------------ printhard: if strpos(fcheck(psfile ,''),'.ps') lt 0 then printx, 'CREATE PLOTFILE FIRST!' else $ psplot,file=psfile goto, opcom ;#------------------------------------------------------------ pspreview: if strpos(fcheck(psfile ,''),'.ps') lt 0 then printx, 'CREATE PLOTFILE FIRST!' else $ xps,psfile goto, opcom ;#------------------------------------------------------------ ; Describe steps in spectral processing steps_help: read_seqfile, steps, concat_dir('SSWDB_SPEX','spex_steps.txt') for i=0,n_elements(steps) - 1 do printx, steps(i) goto, opcom ;#------------------------------------------------------------ ;Locate the filenames consistent with the specified directory and filename ;save them in the string array, found_files. ;This is mainly to interactively select the strings used to find data and drm files find_datafiles: printx,'Using Preview command to locate data files, use Build to create DRMs' preview: in_files = [_1file,_2file] spex_preview, data_tipe, flare, det_id, in_files, directory, text_out=text, $ dir_queue=dir_queue, standard_types=standard_types, source_radec=source_radec, $ files, event_time= avg( utime([start_time_data,end_time_data])), $ det_cosines, dsel, last_found, command, batse_burst=batse_burst checkvar, text, '' if n_elements(files) eq 0 then files='' if fcheck(files(0),'') ne '' then begin text = [text, string(/print,'Data file name returned: ',files) ] endif if n_elements(text) ge 1 then printx, text if flare ne 0 and strupcase(data_tipe(0)) eq 'BATSE' then begin call_procedure,"read_flare", flare, fld def_tb1= fld.bkg1_start_secs+[0.,fld.bkg1_duration] def_tb2= fld.bkg2_start_secs+[0.,fld.bkg2_duration] endif else begin def_tb1= fltarr(2) def_tb2= fltarr(2) endelse def_tb=reform([def_tb1,def_tb2],2,2) if command(0) eq 'time_history' or command(0) eq 'graph' or $ command(0) eq 'background_interval' or command eq 'select_interval' $ then goto, after_preview goto, opcom ;#------------------------------------------------------------ time_history: graph: ; ; Check to see if data specifications have changed before plotting, if so then read ; the data file. ; Check to see if data specifications have changed before trying to find the matching ; file(s). if !d.name eq xdevice() then wsetshow, /TIMEHISTORY if (last_read.data(0) ne strupcase(data_tipe(0)) ) or $ (last_read.data(1) ne strupcase(data_tipe(1))) or $ (last_read.file ne _1file) or $ (last_read.flare ne flare) or $ (last_read.id ne det_id) or $ (need.read eq 1) then begin if (strupcase(data_tipe(0)) eq 'HXRS' and last_read.file ne _1file) or $ (last_read.data(0) ne strupcase(data_tipe(0)) ) or $ (last_read.data(1) ne strupcase(data_tipe(1))) or $ (last_read.flare ne flare) or $ (last_read.id ne det_id) then begin start_time_data='' end_time_data='' endif goto, preview ;find the data files after_preview: command_out=com_param read_4_spex, instrument=data_tipe(0), dformat=data_tipe(1), files=files, $ start_time = start_time_data, end_time=end_time_data, $ det_id=det_id, id_style='ID', flare=flare, p_read_data=p_read_data, $ flux=flux, eflux=eflux, ut=ut, units=uflux, area=area, command_out=command_out, $ ltime=ltime, edges=edges, delta_light=delta_light, wchan=wchan, id_use=id_use, $ mode=datamode, efficiency=eff, cosines=det_cosines, $ auto=fcheck(auto_cal,1), noplot=fcheck(noplot_cal,1), title=title ;Record the last data read last_read = last_found need.back = 1 ;background needs accumulation use_band = -1 ;first background uses all energy channels need.wuse = 1 ;need some channel/energy range default for fitting iegroup = 0 ;grouping factor, no grouping set for new data need.accum = 1 dodrm = getenv('DODRM') if keyword_set(dodrm) then dodrm = fix(dodrm) need.drm = 1 and dodrm need.fit =1 need.tplot=1 need.splot=1 if fcheck(inhibit_drm, 0) then begin ;Remove drm command from command_out wdrm=where(strpos(strlowcase(command_out), 'drm') eq -1,ndrm) if ndrm ge 1 then command_out = command_out(wdrm) else command_out='' endif if command(0) ne '' then $ com_stack = [command_out, com_stack] else command = command_out checkvar, det_cosines, fltarr(1)+1.0 cosine = det_cosines(fcheck(id_use,0) < (n_elements(det_cosines)-1)) setutbase, atime(ut(0), /date) ut = ut - getutbase() start_time_data = anytim( start_time_data,/yohkoh) end_time_data = anytim( end_time_data,/yohkoh) edge_products, edges, mean=emdge, width = wedge count_2_flux = delta_light * area rate = flux ;allow graph to work with no bkg subtracted erate = eflux ;Set default energy bands ;wchan are the energy channels that are allowed ;if a channel is not in wchan it is usually below or above a ;discriminator threshold external to the PHA ;Set Erange after reading data erange = limits( edges(*,wchan) ) printx,string(/print,'Note: Erange reset after reading new data to ==>>',erange) ;Divide the range into 4 bands default_ebands = (indgen(5)*n_elements(wchan)/4 + wchan(0)) < max(wchan) edge_products, default_ebands, edges_2=default_ebands default_ebands = [edges(0,default_ebands(0,*)),edges(1,default_ebands(1,*)-1)] if n_elements(wchan) lt 4 then begin default_ebands = default_ebands*0.0 + 1.1*max(edges) default_ebands(0) = edges(*) endif default_ebands(1,0:2)= default_ebands(1,0:2)*.95 + .05*default_ebands(0,0:2) energy_bands = default_ebands(*) ;make sure uts and ute are set to initial plot boundaries getut, uts=uts_test, ute=ute_test ut_lim = limits(ut)+getutbase() printx,string(/print, 'ut_lim= ',atime(ut_lim),'uts_test,ute_test= ',atime([uts_test,ute_test])) if $ ut_lim(0) le uts_test and ut_lim(1) ge ute_test then begin uts=anytim(uts_test, /yohkoh) ute=anytim(ute_test, /yohkoh) endif else begin uts=anytim(ut_lim(0), /yohkoh) ute=anytim(ut_lim(1), /yohkoh) endelse setut, uts = uts, ute=ute need.read = 0 endif spex_thistory, ut, rate, erate, command, data_tipe, edges, $ count_2_flux, uflux, area, title, $ t_hist_mode, energy_bands, scale_bands, th_yrange, th_ytype, $ real_eband = real_eband, edg_units=edg_units, eplots=eplots, $ sampav=fcheck(tavg,iavg), ltime=ltime, linestyle=tlinestyle need.tplot = 0 wsetshow, /noshow ;help,command last_plot(0)=command(0) ;Send back through opcom if there is more to do, but set input_line to the command! if command(0) ne 'time_history' and command(0) ne 'graph' then input_line =command(0) goto, opcom ;#------------------------------------------------------------ background: if need.read or need.tplot then goto, time_history if !d.name eq xdevice() then wsetshow, /TIMEHISTORY spex_background, ut, flux, eflux, rate, erate, back, eback, ltime, $ tb, wback, com_param, check_defaults, command, $ back_order=back_order, title=title, def_tb=def_tb, $; use_band= use_band, need = need, energy_bands=energy_bands, edges = edges ;Special temporary insert for Yohkoh HXT which only has low channel ;during non-flare if strupcase(arr2str(data_tipe)) eq 'YOHKOH,HXT' then $ for i=1,3 do rate(i,*) = datamode*rate(i,*) if command(0) eq 'graph' then begin input_line=command(0) goto, opcom_fast endif goto, opcom ;#------------------------------------------------------------ display_intrvl: if need.select then begin printx, 'No intervals selected. Enter select_interval command.' goto, opcom endif if need.read or need.tplot then goto, time_history if !d.name eq xdevice() then wsetshow, /TIMEHISTORY graphics_page checkvar, xselect, fltarr(2,1) for i=0,n_elements(xselect(0,*))-1 do begin oplot, xselect(0,i)+fltarr(2), crange('y'), linest=1, $ color = fcolor(5) oplot, xselect(1,i)+fltarr(2), crange('y'), linest=2, $ color = fcolor(4) endfor goto, opcom_fast ;#------------------------------------------------------------ select_intrvl: if need.read or need.tplot then goto, time_history if !d.name eq xdevice() then wsetshow, /TIMEHISTORY need.accum = 1 spex_intervals, ut, xselect ,xinput=xselect, sampav=iavg, $ iselect=iselect, style=style_sel, v_styles=['c','d','r'], query=check_defaults,color=fcolor(9), $ xoffset= (indgen(3)*200+200),yoffset = (indgen(3)*200+200),$ message='SELECT INTERVALS FOR FITTING', error=error, input_variable='xselect' ; printx, 'The vector XSELECT is used for boundaries,' ; printx, 'The format for XSELECT is 2 x n_elements( iselect(0,*) ),' ; help, xselect ; if n_elements(xselect) eq 0 then begin ; printx,'Error in XSELECT ' ; goto, opcom ; endif ; xinput = xselect wint = indgen( n_elements( iselect(0,*) )) ilast= n_elements( iselect(0,*) )-1 ifirst = 0 accum, select = iselect, ltime=ltime, in_obs=flux, in_eobs=eflux, in_back=back, in_eback=eback, $ ;inputs live=live, obsi=obsi,convi=convi, eobsi=eobsi, backi=backi, ebacki=ebacki apar_check, apar, n_elements(wint), apar_arr need.accum = 0 need.select = 0 goto, opcom ;#------------------------------------------------------------ ch_bands: spec_bands: erange: if need.splot then begin printx, 'No spectral plot drawn yet.' goto, bspectrum ;Could have used last_plot(1), ras 9 April 1996 endif if !d.name eq xdevice() then wsetshow, /SPECTRAL graphics_page if command(0) eq 'spec_bands' then goto, spec_bands1 if command(0) eq 'enrange' then goto, erange1 ;""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" ;""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" if strupcase(com_param(0)) eq 'A' then begin ;A stands for AUTO if strupcase(arr2str(data_tipe)) eq 'YOHKOH,HXT' then $ energy_bands = edges(*)*[1.1,.9,1.1,.9,1.1,.9,1.1,.9] goto, opcom endif if n_elements(energy_bands) ne 8 then begin xeselect = fltarr(8) xeselect(0) = energy_bands(*) energy_bands = xeselect endif spex_intervals, edges, xeselect, xinput=reform(energy_bands,2,4), $ xoffset= (indgen(3)*200+200),yoffset = (indgen(3)*200+200), $ iselect=ieselect, style=style_sel, v_styles=['c','d'], query=check_defaults,color=fcolor(9), $ message='SELECT ENERGY BANDS FOR TIME HISTORIES.', error=error, input_variable='energy_bands' energy_bands = ([xeselect(*),fltarr(8)])(0:7) goto, opcom ;#------------------------------------------------------------ spec_bands1: spex_intervals, edges, xeselect, xinput=iegroup, $ xoffset= (indgen(3)*200+200),yoffset = (indgen(3)*200+200), $ iselect=ieselect, style=style_sel, v_styles=['c','d','r'], query=check_defaults,color=fcolor(9), $ message='SELECT ENERGY BANDS FOR SPECTRAL PLOTS.', sampav=ieavg, $ error=error, input_variable='iegroup' iegroup = ieselect goto, opcom ;#------------------------------------------------------------ erange1: If keyword_set( com_param ) then begin spex_energy_ranges, float( com_param ) endif else begin spex_intervals, edges, xeselect, $ xoffset= (indgen(3)*200+200),yoffset = (indgen(3)*200+200), $ iselect=ieselect, style=style_sel, v_styles=['d'], query=check_defaults,color=fcolor(9), $ message='SELECT THE VALID ENERGY RANGE FOR SPECTRAL FITTING.', error=error wuse = ieselect(0,0) + indgen(ieselect(1,0)-ieselect(0,0)) if n_elements(ieselect(0,*)) gt 1 then begin for i=1,n_elements(ieselect(0,*))-1 do $ wuse = [wuse, ieselect(0,i) + indgen(ieselect(1,i)-ieselect(0,i))] wuse = wuse( uniq(wuse, sort(wuse) )) endif wl = limits(wuse) > 0 < (n_elements(edges(0,*))-1) erange = [edges(0,wl(0)),edges(1,wl(1))] endelse need.wuse = 0 goto, opcom ;#------------------------------------------------------------ drm: if need.read then begin printx, 'No data read in yet. Enter time_history or graph command.' goto, opcom endif if strupcase(data_tipe(0)) eq 'BATSE' then call_procedure,"set_pendleton", cont_pendleton command_drm = command(0) delta_light_in = delta_light drm_4_spex, command_drm, data_tipe=data_tipe, det_id=det_id, com_param=com_param, need_drm=need.drm, $ high_drm=high_drm, test_drm=test_drm, efficiency=eff, error=error, $ cosine=cosine, flare=flare, dir=ddirectory, sfile=sfile, dfile=dfile, $ edges=edges, e_in=e_in, drm=drm, area=area, $ gain_offset=gain_offset, gain_mult=gain_mult, $ delta_light=delta_light, p_drm_reader=p_drm_reader if total( abs( delta_light_in-delta_light ) ) ne 0.0 then begin nspec = n_elements(ut(0,*)) for i=0,nspec-1 do begin if n_elements(rate) ne 0 then begin rate(*,i) = f_div(rate(*,i)*delta_light_in,delta_light) erate(*,i) = f_div(erate(*,i)*delta_light_in,delta_light) endif if n_elements(back) ne 0 then begin back(*,i) = f_div(back(*,i)*delta_light_in,delta_light) eback(*,i) = f_div(eback(*,i)*delta_light_in,delta_light) endif if n_elements(flux) ne 0 then begin flux(*,i) = f_div(flux(*,i)*delta_light_in,delta_light) eflux(*,i) = f_div(eflux(*,i)*delta_light_in,delta_light) endif if n_elements(obsi) ne 0 then begin obsi(*,i) = f_div(obsi(*,i)*delta_light_in,delta_light) eobsi(*,i) = f_div(eobsi(*,i)*delta_light_in,delta_light) endif if n_elements(backi) ne 0 then begin backi(*,i) = f_div(backi(*,i)*delta_light_in,delta_light) ebacki(*,i) = f_div(ebacki(*,i)*delta_light_in,delta_light) endif endfor endif if not error then begin edge_products, e_in, gmean=em_in, wid=w_in edge_products, edges, mean=emdge, width = wedge ;get_edges, e_in=e_in, edges=edges, /send endif need.drm = 0 goto, opcom ;#------------------------------------------------------------ save_event: spex_summary,/save, sumfile,data_tipe,flare,dflts_file,$ apar_last=apar_last, chi=chi, xselect=xselect, iselect=iselect, $ tb=tb, prm_spex=prm_spex,style_sel=style_sel,$ input_line=input_line,check_defaults=check_defaults goto, opcom restore_event: spex_summary,/restore,sumfile,data_tipe,flare,dflts_file,$ apar_last=apar_last,chi=chi, xselect=xselect, iselect=iselect, $ tb=tb, prm_spex=prm_spex,style_sel=style_sel,$ input_line=input_line,check_defaults=check_defaults goto, opcom ;#------------------------------------------------------------ accum: ;Done automatically after select! ;Sum the intervals selected together accum, select = iselect, ltime=ltime, in_obs=flux, in_eobs=eflux, in_back=back, in_eback=eback, $ ;inputs live=live, obsi=obsi,convi=convi, eobsi=eobsi, backi=backi, ebacki=ebacki need.fit = 1 apar_check, apar, n_elements(wint), apar_arr, apar_last=apar_last need.accum = 0 goto, opcom ;#-------------------------------------------------------------- Bspectrum: if need.back then begin printx, 'No background intervals selected. Enter background command.' goto, opcom_fast endif if !d.name eq xdevice() then wsetshow, /SPECTRAL ;Plot the background spectrum at the time of the first background accum. ytitle = 'Counts' + uflux spex_spec_plot, edges=edges, group=iegroup, cflux=back(*,avg(wback)),$ xtype=sp_xtype, ytype=sp_ytype, delta_light=delta_light, $ units= ytitle, yrange=spyrange, xrange=spxrange, title='Background '+title need.splot = 0 wsetshow, /noshow last_plot(1) = 'bspectrum' ;Send back through opcom if there is more to do, but set input_line to the command! if command(0) ne 'bspectrum' then input_line =command(0) goto, opcom_fast ;#-------------------------------------------------------------- Force_apar: ;Reset the values in the first parameter array, apar_arr to the input parameter field apar_check, apar, n_elements(wint), apar_arr,apar_last=apar_last apar_arr(*,fcheck(ifirst,0)) = apar goto, opcom ;#-------------------------------------------------------------- Set_spectrum: ;entry point photon spectral plotting ;Disable parameter selection for a command = photon_spectrum, noset if strupcase(strmid(com_param(0),0,1)) eq 'N' $ then set = 0 else set = 1 ;enables graphic parameter selection if strupcase(com_param(0)) eq 'CAL_TEST' $ then cal_test = 1 else cal_test = 0 ;enables graphic parameter selection Count_spectrum: ;entry point for count rate spectral plotting if need.read or need.back or need.select or need.drm then begin printx, 'Need to read data, select background and intervals, and do detector' printx, 'response calculations first.' goto, opcom endif if command(0) eq 'count_spectrum' then begin count=1 last_plot(1)='count_spectrum' endif else begin count=0 last_plot(1)='photon_spectrum' endelse ; ;y_spectitle controls whether spectral plots are in photon or count space!!! ; if count then y_spectitle = 'Counts'+ uflux else y_spectitle = 'Photons'+ uflux if need.accum then begin accum, select = iselect, ltime=ltime, in_obs=flux, in_eobs=eflux, $ in_back=back, in_eback=eback, $ ;inputs live=live, obsi=obsi,convi=convi, eobsi=eobsi, backi=backi, ebacki=ebacki need.accum = 0 need.fit = 1 endif if strpos(strupcase(com_param(0)),'NOBA') ne -1 then backscale=0 else backscale=1 apar_check, apar, n_elements(wint), apar_arr, apar_last=apar_last apar = apar_arr(*,ifirst) apar_in = apar ix = ifirst if !d.name eq xdevice() then wsetshow, /SPECTRAL if need.fit then begin chisqr = 0.0 sigmaa = 0.0 endif else begin chisqr = chi(ix) sigmaa= apar_sigma(*,ix) endelse if strupcase(data_tipe(0)+data_tipe(1)) eq 'BATSESHERS' then begin uncert_vec = obsi(*,ix)*0.0 + uncert wkedge = where( edges(0,*) ge 30 and edges(0,*) le 45.,nwkedge) if nwkedge ge 1 then uncert_vec(wkedge) = uncert_vec(wkedge) * 4. ;relax fit around k_edge endif else uncert_vec = uncert spex_spec_plot, edges=edges, drm=drm, delta_light=delta_light, e_in=e_in, $ apar=apar, group=iegroup, cflux=obsi(*,ix), $ ecflux=eobsi(*,ix), eplot=eplots, $ bflux=backi(*,ix)*backscale, $ ebflux=ebacki(*,ix)*backscale, $ trange = xselect(*,ix), interval=ix, uncert=uncert_vec, $ sigmaa=sigmaa, chisqr=chisqr, more_info=more_info, $ cal_test=cal_test, set=set, units=y_spectitle, $ xtype=sp_xtype, ytype=sp_ytype, wuse=wuse, $ yrange=spyrange, xrange=spxrange, title=title wsetshow, /noshow if total(abs(apar_in-apar)) ne 0 then begin ; apar has changed so replot set = 0 ;replot but don't set query for new parameters apar_arr(*,ix) = apar goto, count_spectrum endif need.splot = 0 set = 0 goto, opcom ;#-------------------------------------------------------------- FITTING: if need.read or need.back or need.drm or need.select then begin printx, 'Need to read data, select background and intervals, and do detector' printx, 'response calculations first.' goto, opcom endif if need.accum then begin accum, select = iselect, ltime=ltime, in_obs=flux, in_eobs=eflux, $ in_back=back, in_eback=eback, $ ;inputs live=live, obsi=obsi,convi=convi, eobsi=eobsi, backi=backi, ebacki=ebacki need.accum = 0 need.fit = 1 endif apar_check, apar, n_elements(wint), apar_arr, apar_last=apar_last ; ITER - Number of times in outside loop where conversion factors are ; recomputed ; NITER - Number of times in loop inside Curvefit where parameter ; increments are calculated ; QUIET - If set then informational print statements are suppressed. checkvar, iter, 4 checkvar, niter, 4 checkvar, nolambda, 0 dofit = 1 dospecplot = 1 ;Plotting and Fitting are defaults here, if one is requested, the other is skipped case strupcase(strmid(com_param(0),0,1)) of 'P': dofit = 0 ;if plot is requested, disable fitting 'F': dospecplot = 0 ;i fit is requested, disable plotting else: ;no action take endcase if n_elements(chi) ne n_elements(apar_last(0,*) ) then chi= (apar_last(0,*)*0.0)(*) if n_elements(wint) ne n_elements(apar_last(0,*) ) then wint = indgen( n_elements( iselect(0,*) ) ) if n_elements(apar_sigma) ne n_elements(apar_last) then apar_sigma = apar_last if !d.name eq xdevice() then wsetshow, /SPECTRAL printx, string(/print,'Free Parameters: ', free) for i=ifirst,ilast <(n_elements(wint)-1) do begin ix = wint(i) if sumfit then apar = apar_arr(*,ix) else $ apar = apar_arr(*,ix-1>ifirst) if not dofit then apar = apar_last(*,ix) if strupcase(data_tipe(0)+data_tipe(1)) eq 'BATSESHERS' then begin uncert_vec = obsi(*,ix)*0.0 + uncert wkedge = where( edges(0,*) ge 30 and edges(0,*) le 45.,nwkedge) if nwkedge ge 1 then uncert_vec(wkedge) = uncert_vec(wkedge) * 4. ;relax fit around k_edge endif else uncert_vec = uncert ; if iter + niter ge 1 then begin sigmaa = apar(where(free)) * 0.0 run_curvefit, e_in=e_in, w_in=w_in, drm=drm, f_model=f_model, $ e_out=edges, count_2_flux=count_2_flux, live=live(*,i), wuse=wuse, $ observe=obsi(*,ix), eobserve=eobsi(*,ix),uncert=uncert_vec, $ use_mod_sig=fcheck(use_mod_sig,0), $ background=backi(*,ix), eback=ebacki(*,ix),conversion=conversion, $ iter=iter, niter=niter, quiet=quiet, nolambda=nolambda, nofit=1-dofit, $ chisqr=chisqr, apar=apar, free=free, range=range, sigmaa=sigmaa convi(*,ix) = conversion apar_arr(*,ix) = apar apar_last(*,ix) = apar apar_sigma(*,ix) = apar_sigma(*,ix)*0.0 apar_sigma(where(free),ix) = sigmaa chi(ix) = chisqr ;reduced chi square need.fit = 0 ;all fit parameters have been defined to something! endif if dospecplot then begin graphics_page ;Flash intervals on time history plot ;Remove intervals on time history plot if i gt ifirst then spex_flash, wint(i-1), xselect, /off spex_flash, ix, xselect, /on spex_spec_plot, edges=edges, drm=drm, delta_light=delta_light, e_in=e_in, apar=apar,$ group=iegroup, cflux=obsi(*,ix), ecflux=eobsi(*,ix), bflux=backi(*,ix), eplot=eplots, $ ebflux=ebacki(*,ix), trange = xselect(*,ix), interval=ix, uncert=uncert_vec, $ set=set, units=y_spectitle, xtype=sp_xtype, ytype=sp_ytype, $ wuse=wuse, yrange = spyrange, sigmaa=apar_sigma(*,ix), chisqr=chisqr, more_info=more_info, $ xrange =spxrange, title=title, /residuals wsetshow, /noshow need.splot = 0 alpha_page endif endfor apar = apar_arr(*,ifirst) ;set parameters to first interval of set last_plot(1) = 'fitting' goto,opcom done: wsetshow, /TIMEHISTORY, /noshow ;reset the system variables only if control is really going back to the outside world if not fcheck(noreadinput,0) then begin !p= bang_enter.p !x = bang_enter.x !y = bang_enter.y !z = bang_enter.z need.tplot = 1 need.splot = 1 endif if keyword_set(gui) and n_elements(main_menu) gt 0 then begin ;printx,'LEAVING SPEX' widget_control, main_menu , sensitive=1 endif end