====================================================================================== The 3He GBT Data Reduction Package in IDL: a Tutorial/Cookbook --------------------------------------------------------------------------------------- HOW TO USE THE GBT IDL SINGLE DISH ANALYSIS PACKAGE written by T.M. Bania, Institute for Astrophysical Research Boston University v1.0 Jan 20 2004 v1.1 Feb 8 2004 v1.2 Feb 10 2004 --------------------------------------------------------------------------------------- v2.0 Apr 15 2004 v2.1 1 May 2004 --------------------------------------------------------------------------------------- v3.1 8 July 2004 ======================================================================================= RUNNING THE PACKAGE ======================================================================================= Run the package from the most recent release: v3.1 At Green Bank on the tbania account this is: 1. cd /users/tbania/idl/v3.1 2. source IDL.INIT.V3.1 : execute the idl.init file to set environment variables 3. IDLDE & : Use IDLDE *if* you intend to write your own procedures Use IDL for maximum clarity during execution. IDL is best for the novice. or either will execute the STARTUP.IDL.V3.1 file IDL defined in step 1. Note that the IDL route is NOT run as a background task. This starts your session. For explanations of the package you should look at the various files in the DOCS/ subdirectory of v3.1/. -------- In the discussion below, IDL package commands are written in CAPITALS. -------- ============================================================================= STARTUP.IDL.V3.1 compiles the v3.1 procedures and then runs a script to put the package into a specific configuration. The last line of this file is: start_gbt,'../data/Feb04.avg.acs.fits' This uses SDFITS file '../data/Feb04.avg.acs.fits' to MAKE_ONLINE the ONLINE data file in {gbt_data} format Files in this format have names *.gbt ==> IN ORDER TO WORK THESE SPECIFIC EXAMPLES YOU NEED TO GET A COPY OF THE SDFITS DATA FILE <== ==> ../data/Feb04.avg.acs.fits <== In Green Bank this file lives in: /users/tbania/idl/data/ At BU here is the ouput that you would see: --------------------------------------------------------------------------------------------------- Scanning contents of SDFITS file ==> ../data/Feb04.avg.acs.fits ../data/Feb04.avg.acs.fits Contains 1136 records Spectrum has 4096 data points at record 0 GBT_IDL loaded from directory= /idl/idl/v3.1 Current ONLINE file definition is: ONLINE file will be= ../data/online.gbt LUN= 104 Do you want to make the ONLINE file? (y/n) ----------------------------------------------------------------------------------------- *****> IDL has found 1146 records in the SDFITS file and all of them have 4096 channel spectral line data. Type 'y' (NO QUOTES!) to create the ONLINE data file in {gbt_data} format. Other choices exist but do this until you are more expert. N.B. This choice OVERWRITES any existing online.gbt file so if you want to keep an old one, copy it at the UNIX level before running this package or before saying 'y'.....) --------------------------------------------------------------------------------------------- Input # records per scan, i.e. # subcorrelators in configuration Input # records per scan: 8 MRDFITS: Binary table. 42 columns by 1136 rows. Constructing ONLINE data file from SDFITS file ==> ../data/Feb04.avg.acs.fits ../data/Feb04.avg.acs.fits Contains 1136 records Made ONLINE data file ../data/online.gbt with 1136 records Select scan range for searches is: 0 to 236 Current NSAVE file is: Nsave file= ../data/nsave.dat LUN= 100 NSlog file= ../data/nsave.log LUN= 101 Do you want to make a new NSAVE file? (y/n: n means use the above files) Make NSAVE? (y or n): ----------------------------------------------------------------------------------------- *****> IDL asks you to: Input # records per scan: This number depends on your GBT correlator configuration. In the mode of these test data, one TP scan produces 8 unique 4096 data point spectra: 4 different frequencies at two polarizations each. IDL has made the online.gbt file with 1136 records. The maximum scan number of the data in this file is 236. IDL gives you the opportunity to make an NSAVE file. Without special actions such as copying and renaming a previous NSAVE file, saying 'y' (NO QUOTES!) *overwrites* the old data. ----------------------------------------------------------------------------------------- Do you want to make a new NSAVE file? (y/n: n means use the above files) Make NSAVE? (y or n):y Do you want to OVERWRITE (o) these files or MAKE new ones (n) ? Overwrite or New? (o or n)o Overwriting NSAVE files Made NSAVE.DAT file ../data/nsave.dat with 256 slots Made NSAVE.LOG file ../data/nsave.log Nsave file= ../data/nsave.dat LUN= 100 NSlog file= ../data/nsave.log LUN= 101 Files currently being used are: SDFITS Data file= ../data/Feb04.avg.acs.fits LUN= 104 ONLINE Data file= ../data/online.gbt LUN= 102 OFFLINE Data file= ../data/offline.gbt LUN= 103 Nsave file= ../data/nsave.dat LUN= 100 NSlog file= ../data/nsave.log LUN= 101 PLOT file= ../figs/idl.ps Journal = ../saves/journal.dat SAVE_state= ../saves/state.dat SAVE_procs= ../saves/procs.dat STACK is now empty ONLINE data file: ../data/online.gbt selected for searches ONLINE data file contains 1136 records Select scan range for searches is: 0 to 236 IDL--> ----------------------------------------------------------------------------------------- IDL has finished initilization and has loaded data into the ONLINE file. IDL has taken an SDFITS data file named ../data/Feb04.avg.acs.fits which contains 1136 records and converted the data into the internal IDL data structure format {gbt_data}. These 1136 records each refer to a 4096 channel spectrum. They have been written to a file named ../data/online.dat which has been attached to the ONLINE data file via the command MAKE_ONLINE which was invoked by START_GBT. You can reprint what files are currently attached via the command FILES. You can change the file assignments via the ATTACH command. Issue ATTACH without any arguments to see a menu. Any number of different data files can be made via a MAKE_DATA command. These {gbt_format} files can be ATTACHED to the IDL session. In fact, one can attach two different data files simultaneously, ATTACHing them either to the ONLINE or the OFFLINE data file. The ONLINE file normally will be the real-time GBT data during observing. Two commands, ONLINE and OFFLINE toggle the package's data handling between these two files. AT PRESENT THIS PACKAGE CANNOT ACCOMODATE DATA FORMATS WITH DIFFERING NUMBERS OF SPECTRAL CHANNELS. TO CHANGE THE NUMBER OF CHANNELS YOU MUST EXIT AND RENTER THE PACKAGE WHICH AT STARTUP AUTOMATICALLY DEFINES THE NUMBER OF SPECTRAL CHANNELS FROM THE SDFITS FILE. N.B. not all of the FILES exist at startup. In particular, offline.gbt needs to be either attached from an existing .gbt file or created via MAKE_DATA and then ATTACHed. Too, neither NSAVE.DAT nor NSAVE.LOG files exist at initial startup. You must create them yourself. =========================================================================================== REAL TIME DATA ACCESS DURING OBSERVING REQUIRES UPDATING THE ONLINE.GBT DATA FILE To access data real time whilst observing: 1. In a unique X-window run: sdfits-test -online -mode=avg /home/gbtdata/AGBT_your_project_file /users/you/idl/data/fname SDFITS will run continuously while you are observing, making /users/you/idl/data/fname.avg.XXX.fits files at the end of each scan. These new data are appended to the existing .fits file above. Here XXX is the backend. This package has been tested only with the 'acs' spectrometer. Note that 'fname' is the file you invoked with sdfits-test. It is the input SDFITS file to the GBT_IDL package. 2. When running the package, when new SDFITS data are written to /users/you/idl/data/fname.avg.XXX.fits files (indicated on the Xterm running sdfits-test) you can update the online.gbt data file with the command: UPDATE It is as simple as that. Further details can be had at: ../v3.1/DOCS/REAL_TIME =========================================================================================== HARDCOPIES: See file PRINTING for details. Set printer device via SETPRINTER command. ; ; setprinter,dev select printer for hardcopy output ; -------------- ; ; ==> setprinter,dev if no device input a menu is printed ; ; setprinter,0 can specify device by integer or string ; setprinter,'lpr' both these commands do the same thing ; =========================================================================================== RECIPES "IDL-->" is the command prompt for the following example commands =========================================================================================== IDL-->list,0,10 <==== list information about the first 11 records in the ONLINE file LIST gives all records LIST,begin_rec,end_rec is a saner choice Rec# Source R.A. Dec. Vel Rx Type Pol Scan# Tsys t_intg 0 W3 02 19 4.4 +62 06 11 -40.9 rx1.1 OFF LCP 6 31.0 6.0 1 W3 02 19 4.4 +62 06 11 -40.9 rx2.1 OFF LCP 6 29.3 6.0 2 W3 02 19 4.4 +62 06 11 -40.9 rx1.2 OFF RCP 6 27.5 6.0 3 W3 02 19 4.4 +62 06 11 -40.9 rx2.2 OFF RCP 6 32.7 6.0 4 W3 02 19 4.4 +62 06 11 -40.9 rx3.1 OFF LCP 6 32.8 6.0 5 W3 02 19 4.4 +62 06 11 -40.9 rx4.1 OFF LCP 6 34.2 6.0 6 W3 02 19 4.4 +62 06 11 -40.9 rx3.2 OFF RCP 6 29.7 6.0 7 W3 02 19 4.4 +62 06 11 -40.9 rx4.2 OFF RCP 6 28.3 6.0 8 W3 02 25 44.4 +62 06 11 -40.9 rx1.1 ON: LCP 7 86.5 6.0 9 W3 02 25 44.4 +62 06 11 -40.9 rx2.1 ON: LCP 7 81.6 6.0 10 W3 02 25 44.4 +62 06 11 -40.9 rx1.2 ON: RCP 7 80.4 6.0 Note here the intrinsic data pattern of 3He data: each GBT telescope scan has produced 8 spectra (records): 4 frequency bands with two polarizations per band. We access data at the record level. (This is exactly what CLASS does.) The 'Rx' identifier flags the frequency band: ; 3He+ band rx1 ; H91 alpha band rx2 ; H115 beta band rx3 ; H92 alpha band rx4 -------------------------------------------------------------------------------------------------- IDL-->fetch,1097 ; get ON TP record # 1097 from ONLINE data file ; and calculate a TP PS pair using its associated OFF scan IDL-->show ; plot this spectrum using defaults ===> To get information on the plot in the IDL window: IDL-->hdr S209 04 11 6.7 +51 09 44 232. rx2.1 LCP PS TP Pair ===> To get a full header of information: IDL-->flagon IDL-->hdr S209 04 11 6.7 +51 09 44 232. rx2.1 LCP PS TP Pair LST +06 50 21.7 HA 2.65 ZA 30.70 AZ 307.7 EL 59.3 Fsky 8587.1613 MHz F0 8665.3000 MHz BW 50.00 MHz -> 12.2070 kHz/chan Tcal 2.99 K Tsys 30.6 K Tintg 8.0 min Vsrc -49.30 km/sec ===> A comment/history string can go here. Currently it is a string='' <=== TsysON 32.1 K TsysOFF 29.1 K TC 3.0 K Elev 59.3 ............................................................................... ===> You an toggle the header information density using FLAGON/FLAGOFF commands: IDL-->flagoff IDL-->hdr S209 04 11 6.7 +51 09 44 232. rx2.1 LCP PS TP Pair IDL-->flagon & hdr & <== Note syntax for issuing multiple IDL commands on command line S209 04 11 6.7 +51 09 44 232. rx2.1 LCP PS TP Pair LST +06 50 21.7 HA 2.65 ZA 30.70 AZ 307.7 EL 59.3 Fsky 8587.1613 MHz F0 8665.3000 MHz BW 50.00 MHz -> 12.2070 kHz/chan Tcal 2.99 K Tsys 30.6 K Tintg 8.0 min Vsrc -49.30 km/sec -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- You are looking at the H91alpha spectrum for S209. The number '232' at the top left corner is the scan number. ==> at this point you might try playing with this spectrum using your UNIPOPS instincts: think: nrset,# bb,# and gg,# for example..... Let's zoom in on the data using the SETX and SETY procs which set the axis range to display. if CURON you do this with the cursor if CUROFF input the values IDL-->curoff IDL-->setx,1000,3000 <== change x-axis range of the display IDL-->show =========================================================================================== FITTING A POLYNOMIAL BASELINE WITH THE CURSOR: =========================================================================================== IDL-->curon <== make sure cursor is on. it toggles via CURON/CUROFF so its state changes only when one of these is issued IDL-->nrset,2 <== set the NREGIONs for the baseline fit: sets 2 regions using cursor. flag from left to right in spectrum via left mouse. As in UniPOPS, NREGIONs are zones to be fitted. set nreg= 1 1027.00 2063.00 <== x-axis values of NREGIONs you defined set nreg= 2 2304.00 2972.00 % Compiled module: STDDEV. % Compiled module: MOMENT. RMS in NREGIONs = 0.0168774 <== RMS of the NREGIONs in y-axis units fit a fifth order baseline to the spectrum using the NREGIONs for the fit: you can use commands B, BB, and BBB. BBB is the super, heavy duty procedure written by John Lyon. bbb,5 ; fits 5th order polynomial and overplots the baseline fit ; issue to SHOW spectrum with baseline removed. IDL-->bbb,5 if you want to plot a zero line toggle with zlon and zloff IDL-->zlon <== ZLON/ZLOFF toggle zero line for SHOW plots IDL-->show if you want to flag a channel, flag,chan# IDL-->flag,2300 to flag radio recombination lines in this spectrum: IDL-->show IDL-->h91 =========================================================================================== FITTING A GAUSSIAN WITH THE CURSOR: =========================================================================================== Gaussian fitting commands are G, GG, and GAUSS IDL-->show Example: fit 2 gaussian components message says to click from left to right: (1) bgauss - begining of region for gaussian fit (2) center of 1st component (3) half-height down the right hand side of the line (4) center of 2nd component (5) half-height down the right hand side of the line (6) egauss - end of region for gaussian fit IDL-->show IDL-->gg,1 <==== ALAS THERE IS ONLY ONE SPECTRAL LINE IN THESE DATA .................................................................................... Click: bgauss #_gauss(center,hw_right) egauss Cursor is at Xpos= 2057.3529 Ypos= -0.0020077063 in PLOT units Cursor is at Xpos= 2192.6471 Ypos= 0.15318332 in PLOT units Cursor is at Xpos= 2219.1176 Ypos= 0.090562383 in PLOT units Cursor is at Xpos= 2313.2353 Ypos= 3.4280941e-05 in PLOT units % Compiled module: CURVEFIT. G# Height +/- sigma Center +/- sigma FWHM +/- sigma ==> kHz km/sec 0 0.167 0.003 2188.010 0.603 71.811 1.421 876.595 30.327 Center/FWHM in channels RMS in NREGIONs = 0.0168096 % Program caused arithmetic error: Floating underflow .................................................................................... ==> this error comes and goes. number is too small to be properly set in floating point. by eye the line fit looks good. do not worry. ==> i am tracking down the issue of the reported errors in the fit... SHOW is overplotted by: white line is the original data red line is the gaussian fit cavalier orange line is the difference between data - model <== can be toggled via FLAGON/FLAGOFF A SHOW command demonstrates that the original data are still in buffer 0 (!b[0] system variable structure) the Gmodel is stored in buffer 1 and the difference, Gmodel - Data, is stored in buffer 8. (See OVERVIEW for buffer definition and use.) =========================================================================================== DATA FILTERING ---- the STACK =========================================================================================== * There is a STACK for filtering on the data file and other things. The following uses the command SET to SELECT data with specific attributes from the ONLINE data file. IDL-->set Source Name= no quotes necessary Set Source Name: S209 Scan Type = ON or OFF Set Scan Type: ON Line ID = e.g. rx1.1, rx2.2, rx1, rx2, .... Set Scan ID: rx2 <== H91alpha band, both RCP+LCP Scan # = Set Scan #: IDL-->flagoff <== supress output from SELECT IDL-->select Selecting: Source= S209 Type= ON ID= rx2 Scan= * Record range for search= 0 to 1136 Scan range for search= 0 to 236 from ONLINE data file= ../data/online.gbt EOF on data file Found 106 matching records in ONLINE data file= ../data/online.gbt IDL--> ===> There are 106 ON/OFF TP pairs that match this search. The ON record numbers are stored in the STACK. To show the contents of the STACK: IDL-->tellstk STACK has !acount= 106 265 267 281 283 297 299 321 323 345 347 361 363 393 395 409 411 425 427 441 443 457 459 473 475 489 491 505 507 521 523 537 539 553 555 569 571 585 587 601 603 617 619 633 635 649 651 665 667 681 683 697 699 713 715 729 731 745 747 761 763 777 779 793 795 809 811 825 827 841 843 857 859 873 875 889 891 905 907 921 923 937 939 953 955 969 971 985 987 1001 1003 1017 1019 1033 1035 1049 1051 1065 1067 1081 1083 1097 1099 1113 1115 1129 1131 IDL--> You can get information about the STACK with CATalog: IDL-->cat Rec# Source R.A. Dec. Vel Rx Type Pol Scan# Tsys t_intg 265 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 ON: LCP 69 32.0 6.0 267 S209 04 11 6.7 +51 09 44 -49.3 rx2.2 ON: RCP 69 30.0 6.0 281 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 ON: LCP 71 32.0 6.0 283 S209 04 11 6.7 +51 09 44 -49.3 rx2.2 ON: RCP 71 30.0 6.0 297 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 ON: LCP 73 32.0 6.0 299 S209 04 11 6.7 +51 09 44 -49.3 rx2.2 ON: RCP 73 30.0 6.0 etc. STACK MANIPULATION PROCEDURES To show contents of the STACK: IDL-->tellstk STACK has !acount= 3 1017 1019 1033 To clear the STACK: IDL-->clrstk STACK is now empty IDL-->tellstk STACK is empty To delete a record from the STACK: IDL-->tellstk STACK has !acount= 3 1017 1019 1033 IDL-->sub,1019 <== SUBtract record_number from the STACK IDL-->tellstk STACK has !acount= 2 1017 1033 There are 3 approaches to filling the STACK: (1) clrstk add,1065 add,1067 etc (2) clrstk addem,[1065,1067,1081,1083,1097,1099] <-- brackets required as we are passing an array of values (3) SET and SELECT IDL-->clrstk STACK is now empty IDL-->setrec,1000,9999 Select record range for searches is: 1000 to 1135 IDL-->set Source Name= no quotes necessary Set Source Name: S209 <=== or nulls set wildcard '*' match Scan Type = ON or OFF Set Scan Type: ON Line ID = e.g. rx1.1, rx2.2, rx1, rx2, .... <=== rx2 will match both LCP+RCP Set Scan ID: rx2 Scan # = Set Scan #: Selecting: Source= S209 Type= ON ID= rx2 Scan= * IDL-->select Selecting: Source= S209 Type= ON ID= rx2 Scan= * from ONLINE data file= /idl/idl/data/online.gbt Found 18 matching records in ONLINE data file= /idl/idl/data/online.gbt IDL-->tellstk STACK has !acount= 18 1001 1003 1017 1019 1033 1035 1049 1051 1065 1067 1081 1083 1097 1099 1113 1115 1129 1131 **************************************************************************************** Selection procedures for SELECT filtering of ONLINE/OFFLINE {gbt_data} files **************************************************************************************** set -> sets parameters with walk through dialog setrec,recmin,recmax -> sets record range for search setrange,scanmin,scanmax -> sets scan range for search setsrc,'NGC3242' -> takes string argument and sets source name settyp,'ON' -> takes string argument and sets scan type setid,'rx1' -> takes string argument and sets line identification setscan,scn# -> takes integer argument and sets scan number clearset -> sets all SELECT parameters to the wildcard '*' select -> searches ONLINE or OFFLINE for matches and fills STACK with those record numbers cat -> displays information about the records currently in the STACK online/offline -> toggle to switch search data file between the ONLINE and OFFLINE {gbt_data} files =========================================================================================== DATA AVERAGING =========================================================================================== * Suppose we have the following 3 ON records in the stack: IDL-->tellstk STACK has !acount= 3 1017 1019 1033 IDL-->cat Rec# Source R.A. Dec. Vel Rx Type Pol Scan# Tsys t_intg 1017 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 ON: LCP 217 32.1 4.0 1019 S209 04 11 6.7 +51 09 44 -49.3 rx2.2 ON: RCP 217 30.3 4.0 1033 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 ON: LCP 219 32.0 4.0 * you can average these via the AVGSTK command. The nicest output comes if FLAGOFF. IDL-->avgstk ........................................................................ 3 records averaged S209 04 11 6.7 +51 09 44 217. rx2.1 LCP PS Average TsysON 31.5 K TsysOFF 28.4 K TC 3.1 K Elev 67.3 ......................................................................... IDL-->show ==> to average the two polarizations for scan 217 above: fetch,1017 <== get and calculate a ON/OFF TP pair accum <== ACCUMulate the data in averaging buffers fetch,1019 accum ave <== AVErage the ACCUMulated data show =========================================================================================== NSAVE files =========================================================================================== * As in Unipops you can save (and recall) the !b[0] buffer (header and data) in an NSAVE file. You must either create (via MAKE_NSAVE) or ATTACH the NSAVE.DAT and NSAVE.LOG files. Note that you can create many of these with different names and ATTACH at will.... At present each NSAVE file holds only 100 spectra. This is easily changed, but that might not be a good idea. NSAVE Procedures: make_nsave : run this procedure only ONCE. currently it hardwires the file ---------- name so running it again will OVERWRITE the data in this file....! putns : copies the !b[0] data structure into position !nsave in the !nsavefile ----- if !nsave_log(!nsave)=1 this slot has data in it then if !protectns=1 procedure will NOT overwrite getns : copies data in nsave(!nsave) position into !b[0] ----- nson and nsoff ; toggle value of !protectns from 1 to 0 respectively --- ----- tagid ----- tagpol ------ tagtype ------- comment : puts string into !b[0].history ------- nslog : searches NSAVE file for slots where data has been written and ----- ouputs information, example: IDL-->nslog NSAVE file /idl/idl/data/nsave.dat contains: NS# Source R.A. Dec. Vel Rx Pol Type Tsys Tintg 8 S209 04 11 6.7 +51 09 44 -49.3 rx1.1 LCP PS Average 31.4 96.1 9 S209 04 11 6.7 +51 09 44 -49.3 rx1.1 LCP PS Average 30.6 776.8 10 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 LCP PS Average 30.4 96.1 11 S209 04 11 6.7 +51 09 44 -49.3 rx2.1 L+R PS Average 29.6 776.8 12 S209 04 11 6.7 +51 09 44 -49.3 rx3.1 LCP PS Average 33.2 96.1 13 S209 04 11 6.7 +51 09 44 -49.3 rx3.1 LCP PS Average 32.6 776.8 14 S209 04 11 6.7 +51 09 44 -49.3 rx4.1 LCP PS Average 33.4 96.1 15 S209 04 11 6.7 +51 09 44 -49.3 rx4.1 LCP PS Average 32.7 776.8 =========================================================================================== ===========================================================================================