NAME:
acorexample2 - Looking at calibration data by month.
At the end of each month, all of the calibration scans for the
month are processesed and stored in an idl archive file. This includes
x102, x101,x113,gui calibrate button runs, etc.. . The processing includes
the standard mm_mueller0 processing (which applies the current mueller
matrix to correct the data). It does not run mueller2_5 (these routines
would recompute the mueller matrix elements).
The output of the processing is an array mm[] of {mueller}
structures (one for each pattern done). See mmrestore() for a
description of this structure.
The data is stored in the directory:
/share/megs/phil/x101/x102/runs
The files are:
cyymmdd1_yymmdd2.dat - these are ascii files with a list of the files
and scans within each file for the month.
yymmdd1 is the first date of the file,
yymmdd2 is the last date of the file
cyymmdd1_yymmdd2.sav - these are the idl save files that contain the
processed data.
rcvrN.dat - An ascii file containing a list of all of the scans for
receiver N for all of the data files.
As of 15apr02 the list of save files was:
c010101_010131.sav c010501_010531.sav c010901_010930.sav c020101_020131.sav
c010201_010228.sav c010601_010630.sav c011001_011031.sav c020201_020228.sav
c010301_010331.sav c010701_010731.sav c011101_011130.sav c020301_020331.sav
c010401_010430.sav c010801_010831.sav c011201_011231.sav c020401_020415.sav
Use the routine mmgetarchive to retrieve all or a subset of this data.
eg:
nrecs=mmgetarchive(010501,020501,mm,rcvnum=5)
You can then create a subset of the data with mmget.
You can plot the data with:
mmplotgtsb .. gain,tsys,sefd, beamwidth
mmplotcsme .. coma, sidelobes,mainbeam efficiency
mmplotpnterr .. pointing error
(See /pkg/rsi/local/libao/phil/Cor2/acorexample2.pro)
NAME:
arch_getdata - get cor data using the archive tbl arrays
SYNTAX: n=arch_getdata(slAr,slfileAr,indAr,b,type=type,incompat=incompat,
han=han,missing=missing)
ARGS:
slAr[l] : {sl} or {corsl} returned by arch_gettbl
slFileAr[m] : {slInd} returned by arch_gettbl
indAr[] : long indices into slar to return
KEYWORDS:
type : int 0 first rec (hdr and data) of each scan1
: 1 hdrs from first rec of each scan
: 2 all recs each scan (hdr and data)
: 3 average rec (hdr and data) each scan
: 4 all hdrs from scan
RETURNS:
l : int number of elements in b
b : depending on type it can be {corget} or {hdr}
slind[n]: long the index into slar for each element returned in b.
incompat[p] long indices in indAr that were not returned because the
datatype differs from that of the first record
missing[m] long indices in indAr that were not in the file location
that the database had recorded.
DESCRIPTION:
After using arch_gettbl() and possibly where(), call this routine
to read the header and data from disc. What you get back is determined by
the keyword type and ind. The number of elements in b[] can be greater
than the number of indices in ind[] (eg you asked for all of the records
of the scans, or you are returning just headers). The slind keyword
array has the same number of elements as b. It contains the index into
slAr for each elements of b.
EXAMPLES
;get all data for jan02->apr02 cband
nscans=arch_gettbl(020101,020430,slAr,slFileAr,rcv=9)
; get the start of all the on/off pairs with cal on /off
pattype= 1 ; on/off position with cal on,off
nfound=corfindpat(slar,indar,pattype=pattype)
n=arch_getdata(slar,slfilear,indar,bon,type=3,incompat=incompat)
n=arch_getdata(slar,slfilear,indar+1,boff,type=3,incompat=incompat)
; get just the headers from the on scan of each calon/off pair
pattype= 6 ; on/off position with cal on,off
nfound=corfindpat(slar,indar,pattype=pattype)
n=arch_getdata(slar,slfilear,indar,hdrcon,type=1,incompat=incompat,
slind=slind)
NOTE:
When returning just headers, each header of each board is returned
as a separate entry in b[]. Use slind to figure out which scan each
belongs to.
(See /pkg/rsi/local/libao/phil/Cor2/arch_getdata.pro)
NAME:
arch_getmap - input a map from the archive
SYNTAX: nstrips=arch_getmap(yymmdd1,yymmdd2,projId,srcNm,m,cals,$
calSpc,procNm=procNm,polABrd=polABrd,polBBrd=polBBrd,$
rcvNum=rcvNum,smpPerStrip=smpPerStrip,norev=norev,$
avgsmp=avgsmp,han=han,verb=verb,$
useslar=useslar,slar=slar,slfilear=slfilear,badscans=badscans
ARGS:
yymmdd1: long starting date in the archive to search:eg 030501 .
yymmdd2: long ending date in the archive to search:eg 030801 .
projId : string project id you used for the map (eg 'a1731').
srcNm : string source name used for the map.
KEYWORDS:
procNm: string procedure name used for taking data. Supported names
are: 'cormap1','cormapdec','cordrift'. The default is
'cormap1'
polABrd: int correlator board number to take polA (1 thru 4). default
is 1.
polBBrd: int correlator board number to take polB (1 thru 4). default
is 1. (This assumes two pols per board).
rcvNum: int restrict search to this receiver number. The default
is all receivers found.
smpPerStrip:long samples per strip in map. Default is to read it from
the first strip found. Use this parameter if the
first strip found is not complete.
norev: When you drive both directions, the routine will normally
reverse the odd strips so that the data,header, and cals of
adjacent strips will line up in RA. Setting norev to true
will override this default. The data will be returned in
the order it is taken. The default is false.
avgsmp: long Number of samples to average together. Handy if you've
oversampled each strip. If this number does not divide
evenly into the samples/strip, then samples at the end
will be dropped. avsmp=0 or 1 is the same as no averaging
The default is no averaging.
han: if set then hanning smooth the data on input. The default
is no hanning smoothing.
useslar: if set the routine will use the slar,slfilear passed in
by the user (rather than the default archive).
slar[]: {corsl} slar passed in by user (if useslar is set) or passed
back by this routine
slfilear[]: {slind} slfilear passed in by user (if useslar is set) or
passed back by this routine.
badscans[]: long if present then exclude these scan. Helps if strip with no cal..
RETURNS:
nstrips: long The number of strips returned. This may be more than the
requested number of strips in the file if the user
redid some of the strips.
m[2,smpPerStrip,nstrips]:{} array containg the map info. There is 1 entry
per strip taken. See cormapinp() for a description of
this structure.
cals[n*nstrips]:{} cal info for each strip. n is 1 for 1 cal per strip
2 for 2 cals per strip (but note the warning below
for 2 cals per strip).
calSpc[n*nstrips]:{corget} if supplied then return the actual cal spectra
for each strip (calOn followed by calOff). Use this for
recomputing the cal values with a different mask.
DESCRIPTION:
On the fly mapping can be done with cormap1,cormapdec, or the cordrift
routines. Users could use cormapinp() or cormapinplist() to input the
map from one or more files. The routine arch_getmap() provides the
same functionality as cormapinplist() with the added benefit that it
locates all of the strips for the map from the data archive.
To identify the map data in the archive you need to enter:
1.The range of dates in the archive to search (yymmdd1,yymmdd2).
2.The project id used for the map.
3.The source name used for the map.
4.If the map was not taken with cormap1 you need to supply
the keyword procNm='cormapdec' or procNm='cordrift'
5.If the same source was mapped with more than one receiver (eg 430,lb)
then you need to enter the receiver number to use. If not it will
combine all of the data it finds. receiver numbers are:
1=327,2=430,3=610,5=lbw,6=lbn,7=sbw,8=sbh,9=cb,11=xb,12=sbn.
6.By default the routine will take polA, and polB from correlator board
1. You can change this with the keyword polABrd=2,polBBrd=3
arch_getmap searches the archive for the locations of the requested
dataset. It creates a list of files and starting scannumbers that it
then passes to cormapinp() (via cormapinplist). cormapinp() will stop
reading a file if it finds a strip with fewer than the requested number of
strips, or it has input the last strip of the map. arch_getmap() is
careful in searching the complete file looking for any strips that the
user has redone. These strips will be included in the returned data.
This makes it possible for the returned array m[] to have more strips
than the number requested in the map. You can find the stripnumber of
each strip from: stripNum=m[0,0,*].h.proc.iar[4]
EXAMPLE:
the source 'ca' was mapped by project a1731 at 430 using cormap1 (which
drives in ra) and cormapdec (which drives in dec). The data was taken
May through august of 2003. Input both maps.
src='ca'
polABrd=1
polBBrd=2
rcvnum=2 ; the 430 dome receiver
yymmdd1=030501
yymmdd2=030801
projId='a1731'
procNm='cormap1'
nstrips=arch_getmap(yymmdd1,yymmdd2,projId,srcNm,mra,procNm=procNm,$
polABrd=polABrd,polBBrd=polBBrd,rcvnum=rcvnum,$
/verb)
procNm='cormapdec'
nstrips=arch_getmap(yymmdd1,yymmdd2,projId,srcNm,mdec,procNm=procNm,$
polABrd=polABrd,polBBrd=polBBrd,rcvnum=rcvnum,$
/verb)
SEE ALSO:
cormapinp(),cormapinplist().
Note:
The routine will now search backwards for the first calon/off of the
map (if there was one).
You need to have an idea how much memory the maps will take so you
don't ask for more memory that your computer can provide. As
an example: the cormapdec map above used:
(1024chan)*(4bytes/chan)*(2pol)*(257smp/strip)*(65 strips) is 137 megaBytes
Remember that the routine will input all of the strips done (including
any repeated strips).
(See /pkg/rsi/local/libao/phil/Cor2/arch_getmap.pro)
NAME:
arch_getmapinfo - get info on maps from the archive
SYNTAX: nmaps=arch_getmapinfo(yymmdd1,yymmdd2,mapI,projId=projid,srcNm=srcNm,$
procNm=procNm,rcvNum=rcvNum,useslar=useslar,$
slar=slar,slfilear=slfilear
ARGS:
yymmdd1: long starting date in the archive to search:eg 030501 .
yymmdd2: long ending date in the archive to search:eg 030801 .
KEYWORDS:
projId : string if provided then limit the maps to this project id.
srcNm : string if provided then limit the maps to this source name.
procNm: string if provided then limit the maps to this data taking
procedure name. Supported names are:
'cormap1','cormapdec','cordrift'.
rcvNum: int restrict search to this receiver number.
useslar: if set the routine will use the slar,slfilear passed in
by the user (rather than the default archive).
slar[]:{corsl} slar passed in by user (if useslar is set) or passed
back by this routine
slfilear[]:{slind} slfilear passed in by user (if useslar is set) or
passed back by this routine.
RETURNS:
nmaps: long The number of maps found
mI[nmaps]:{} array of mapinfo structures, one foreach map
DESCRIPTION:
Mapping can be done with the datataking routines: cormap1,cormapdec,
and cordrift. The routine cormapinp() can be used to input a single
map if the user knows the filename and scannumber. The routine
arch_getmap() can be used to retrieve all map data given a number of
constraints (date range, projid, srcname, etc..). This data is returned
as an array of strips. For some projects, the amount of mapping data is
too large to recall at once. In this case, arch_getmapinfo() can be
used to get a summary of the maps that are in the archive, without
retrieving the actual data. You can than use the summary info
to decide which maps to input.
arch_getmapinfo() defines a map as a set of scans in a mapping procedure
that are contiguous in time and have increasing strip numbers. It returns
an array of map Info strucutures (one for each map found).
If a map has 120 samples per strip and 36 strips and the user
did one complete map and a second map of strips 1 through 32. This
routine would return a mapInfo structures with two entries.
The mapInfo structure contains:
** Structure <8266244>, 12 tags, length=1052, data length=1050, refs=1:
SRCNM STRING 'MRK335' Name of source
PROCNM STRING 'cordrift' Name of datataking procedure
SMPSTRIP INT 120 requested samples in a strip
REQSTRIPS INT 36 Requested strips in a map
RACENHR FLOAT 0.105417 Ra center in hours
DECCENDEG FLOAT 19.9528 Dec center in degrees
FIRSTSTRIP INT 1 First strip returned(count from 1)
NUMSTRIPS INT 36 number of strips returned
FNAME STRING '/proj/a1552/corfile.15nov01.a1552.1'
SCANNUM LONG 131900007 scan number 1st sample,1st strip
COMPLETEMAP INT 1 map is complete
HDR STRUCT -> HDR Array[1] hdr brd 1, 1st sample,1st strip
EXAMPLE:
the source 'MRK335' was mapped by project a1552 at lband using cordrift.
The data was taken nov01,dec01. To get the mapinfo data:
srcnM='MRK335'
projId='a1552'
procNm='cordrift'
yymmdd1=011101
yymmdd2=031231 ; search over a larger date range
nmaps=arch_getmapinfo(yymmdd1,yymmdd2,mI,projId=projId,srcNm=srcNm,$
procNm=procNm)
help,mI
MI STRUCT = -> Array[29] ; there were 29 maps
It turns out that the maps had two separate centers. You could find
them by looking at racenhr,deccendeg:
SEE ALSO:
arch_getmap(),cormapinp(),cormapinplist().
(See /pkg/rsi/local/libao/phil/Cor2/arch_getmapinfo.pro)
NAME:
arch_getonoff - get on/off -1 data from archive
SYNTAX: nfound=arch_getonoff(slAr,slfileAr,indar,bout,infoAr,$
incompat=incompat,han=han,scljy=scljy,verbose=verbose)
ARGS:
slAr[n] : {sl} or {corsl} returned by arch_gettbl
slFileAr[m] : {slInd} returned by arch_gettbl
indar[] : long indices into slar for scans to return
KEYWORDS: (for corposonoff)
han: if set then hanning smooth
scljy: if set then scale to jy using gain curve for date of scan.
if not set then return in Kelvins.
verbose:if set then call correcinfo for each scan we process
RETURNS:
nfound : long number of on/off-1 returned in bout
infoAr[2,n]: long info on each on/off -1 found. Values are:
info[0,n] status of the returned on/off's
1 - on/off -1 returned as requested.
2 - on/off -1 but only the off was in ind[]. This could
occur if you asked for all the scans within the radius
of some ra/dec and the off fell within the radius
but not the on. You might then look for a negative
going galaxy in the on/off-1.
<0- no cal was found so the units are Tsys.
returnes -1 or -2
info[1,n] this is the index into indar[] where this on/off
was found (could be the on or the off).
bout[n] : on/off-1 scaled to jy, K, or tsys (see info[]). There is
one elemetn for each on/off pair found
incompat[p] long indices in indAr that were not returned because the
datatype differs from that of the first on/off pair.
Since bout[] is an array, each element must have the same
format (eg. number of boards, lags/sbc, etc..).
DESCRIPTION:
After using arch_gettbl() and possibly where(), call this routine
to process all of the on,off position switch records included in the
subset slar[ind]. It will search for all of the on or offs in ind[] and
process them via corposonoff(). If only an ON or an OFF is found in ind[]
it will still try to process the pair. By default it will use the cal
to return the data in kelvins. The /scljy will return the data in janskies.
If no cal scan is found, then the data is returned in units of Tsys.
The infoar[2,*] returns information on the returned data. The first
index tells the format of the data, the second index is the ptr into
slar where the onPos was found.
EXAMPLES
;get all data for jan02->apr02 cband
nscans=arch_gettbl(020101,020430,slAr,slFileAr,rcv=9)
; select all the on/off pairs
pattype= 1 ; on/off position with cal on,off
nfound=corfindpat(slar,indar,pattype=pattype
n=arch_getonoff(slar,slfilear,indar,b,type=3,incompat=incompat)
NOTE:
The routine expects to find the on,off, and cals in the slar. You should
pass the entire slar with indar as the thing that does the subsetting.
Don't pass a subset of slar into this routine: eg..
arch_gettbl(010101,011230,slar,slfilear,rcv=6)
ra =131215.1
dec=101112.
dist=cmpcordist(ra,dec,slar=slar)
; find all points within 30 arcminutes of ra,dec from projid a1199
; ..wrong..
ind=where(dist[2,*] lt 30.,count)
slarn=slar[ind]
ind=where(slarn.projid eq 'a1199',count)
npts=arch_getonoff(slarn,slfilear,ind,/scljy)
; .. correct:
ind=where((dist[2,*] lt 30.) and (slar.projid eq 'a1199'),count)
npts=arch_getonoff(slar,slfilear,ind,/scljy)
(See /pkg/rsi/local/libao/phil/Cor2/arch_getonoff.pro)
NAME:
arch_gettbl - input table of scans from cor archive
SYNTAX: nscans=arch_gettbl(yymmdd1,yymmdd2,slAr,slfileAr,rcvnum=rcvnum,
freq=freq,proj=proj,cor=cor)
ARGS:
yymmdd1 : long year,month,day of first day to get (ast)
for interim correlator data yy is two digits.
for was data, yy is 4 digits 2004mmdd
yymmdd2 : long year,month,day of last day to get (ast)
for interim correlator data yy is two digits.
for was data, yy is 4 digits 2004mmdd
KEYWORDS:
rcvnum : long .. receiver number to extract:
1=327,2=430,3=610,5=lbw,6=lbn,7=sbw,8=sbw,9=cb,$
10=xb,12=sbn,100=430ch
freq[2] : float Require cfr of band to be between freq[0] and freq[1]
in Mhz. At least 1 sbc of scan must match this.
proj : string project number to match (eg. 'a1473')
cor : If set then return the cor scanlist (see below)
RETURNS:
slAr[count] : {sl} or {corsl} scanlist array for each scan
slFileAr[n] : {slInd} one entry per file found
nscans : long number of scans found
DESCRIPTION:
This routine will scan the correlator archive and return a table
(an array of scanlist structures) for all scans in the archive that meet
the requested criteria. Two arrays are returned:
1. slAr[] holds 1 entry per scan
2. slfileAr[m] holds 1 entry per file found.
-------------------------------------------------------------------
For Interim correlator data:
Each slAr entry is an {sl} or {corsl} structure that contains:
scan : 0L, $; scannumber this entry
bytepos : 0L, $; byte pos start of this scan
fileindex : 0L, $; lets you point to a filename array
stat : 0B ,$; not used yet..
rcvnum : 0B ,$; receiver number 1-16
numfrq : 0B ,$; number of freq,cor boards used this scan
rectype : 0B ,$;type of pattern used to take the data. The
values are listed below (see getsl() under
generic idl routines for the most recent list).
1-calon,2-caloff,3-posOn,4-posOff,5-coron
6-cormap1,7-cormapdec,8-cordrift,
9-corcrossch (used by calibration routines)
10-x111auto (rfi monitoring),11-one(bml),
12-onoffbml
... see getsl() in GENEral idl
numrecs : 0L ,$; number of groups(records in scan)
freq : fltarr(4),$;topocentric freqMhz center each subband
julday : 0.D,$; julian date start of scan
srcname : ' ',$;source name (max 12 long)
procname : ' '};procedure name used.
If the /cor keyword is included then the slAr will be an {corsl}
structure and will contain the following additional fields:
projId : '' ,$; from the filename
patId : 0L ,$; groups scans beloging to a known pattern
secsPerrec : 0. ,$; seconds integration per record
channels : intarr(4),$; number output channels each sbc
bwNum : bytarr(4),$; bandwidth number 1=50Mhz,2=25Mhz...
lagConfig : bytarr(4),$; lag config each sbc
lag0 : fltarr(2,4),$; lag 0 power ratio (scan average)
blanking : 0B ,$; 0 or 1
azavg : 0. ,$; actual encoder azimuth average of scan
zaavg : 0. ,$; actual encoder za average of scan
raHrReq : 0.D,$; ra hours requested , start of scan
decDReq : 0.D,$; dec degrees requested, start of scan J2000
Delta ScanEnd-scanStart.. real angle for motion
raDelta : 0. ,$; delta ra Arcminutes real angle
decDelta : 0. ,$; delta dec Arcminutes
azErrAsec : 0. ,$; avg azErr Asecs great circle
zaErrAsec : 0. $; avg zaErr Asecs great circle
-------------------------------------------------------------------
For was2 data:
The return slar data structure is:
a={slwas ,$
scan : 0L, $; scannumber this entry
rowStart : 0L, $; row in fits file start of scan 0 based.
fileindex : 0L, $; lets you point to a filename array
stat : 0B ,$; not used yet..
rcvnum : 0B ,$; receiver number 1-16, 17=alfa
numfrq : 0B ,$; number of freq,cor boards used this scan
rectype : 0B ,$;1-calon,2-caloff,3-posOn,4-posOff
numrecs : 0L ,$; number of groups(records in scan)
freq : fltarr(8),$;topocentric freqMhz center each subband
julday : 0.D,$; julian day start of scan
srcname : ' ',$;source name (max 12 long)
procname : ' ',$;procedure name used.
stepName : ' ',$;name of step in procedure this scan
projId : '' ,$; from the filename
patId : 0L ,$; groups scans beloging to a known pattern
secsPerrec : 0. ,$; seconds integration per record
channels : intarr(8),$; number output channels each sbc
bw : fltarr(8),$; bandwidth used Mhz
backendmode: strarr(8),$; lag config each sbc
lag0 : fltarr(2,8),$; lag 0 power ratio (scan average)
blanking : 0B ,$; 0 or 1
azavg : 0. ,$; actual encoder azimuth average of scan
zaavg : 0. ,$; actual encoder za average of scan
encTmSt : 0. , $; secs Midnite ast when encoders read
; start of scan
raHrReq : 0.D,$; requested ra , start of scan J2000
decDReq : 0.D,$; requested dec, start of scan J2000.
; Delta end-start real angle for requested position
raDelta : 0. ,$; delta ra last-first recs. Amins real angle
decDelta : 0. ,$; delta dec (last-frist)Arcminutes real angle
pntErrAsec : 0. ,$; avg great circle pnt error
; alfa related
alfaAngle: 0. , $; alfa rotation angle used in deg
alfaCen :
You can use the slAr and slFileAr to then access the actual data files
without having to search through the file.
EXAMPLES
The following examples will read the scanlist array for all the data
beteen jan02 and apr02. The examples will then select different pieces
of information that are in this scanlist array. The final step will
extract some of the actual datascans.
;get tbl for all data for jan02->apr02.
nscans=arch_gettbl(020101,020430,slAr,slFileAr)
; work with this scanlist dataset..
; make a list of the unique source names
srcnames=slar[uniq(slar.srcname,sort(slar.srcname))].srcname
print,srcnames
;
; find all of the on/off position switch data for cband (rcv=9)
; (note corfindpat does not yet work for was2 data..)
n=corfindpat(slar,indar,pattype=1,rcv=9)
; find source NGC2264
indar=where(slar.srcname eq 'NGC2264',count)
now extract the scan averaged data. It will extract all data scans
that match the dataformat of the first entry of indar[0] (see arch_getdata().
n=arch_getdata(slar,slfilear,indar,b,type=3,/han,incompat=incompat)
SEE ALSO:getsl, arch_getdata,corfindpat.
NOTES:
22aug04. this is now implemented for was2 data. The cor= keyword is
ignored. The data stucture returned is a little different than the
interim corr structure.
(See /pkg/rsi/local/libao/phil/Cor2/arch_gettbl.pro)
NAME:
arch_openfile - open a file using archive arrays
SYNTAX: stat=arch_openfile(slAr,slfileAr,ind,lun)
ARGS:
slAr[l] : {sl} or {corsl} {slwas} returned by arch_gettbl
slFileAr[m] : {slInd} returned by arch_gettbl
ind : long index into slar for file to open
RETURNS:
stat: int 1 opened ok, 0 could not open
lun: int or {desc} for the open file
DESCRIPTION:
arch_openfile() will open a file that contains a given scan.
The user inputs slar[], slfilear[] that are returned by arch_gettbl()
and an index into slar[] for the scan that you want. The routine will
then find the file that this scan is in and open it. It returns the
lun (or descriptor if this is a was file) that can then be used to
do i/o on the file. The routine does not read any data from the file.
When done be sure and free the lun/desc using free_lun,lun or
wasclose,lun
EXAMPLES
For was data..
;get table for data august 20 thru august 31 2004.
nscans=arch_gettbl(20040820,20040831,slAr,slFileAr)
; find scan number=423314477L
scan=423314477L
ind=where(slar.scan eq scan,count)
if count eq 0 then begin
print,'scan:',scan,' not found'
goto,done
endif else begin
istat=slar_openfile(slar,slfilear,ind[0],desc)
if istat ne 1 then goto,done
endelse
;
; now input the scan
;
istat=corinpscan(desc,b,scan=scan)
(See /pkg/rsi/local/libao/phil/Cor2/arch_openfile.pro)
NAME:
chkcalonoff - check that hdrs are a valid cal on,off.
SYNTAX: istat=chkcalonoff(hdrOn,hdrOff)
ARGS:
lun: assigned to file we are reading
hdrOn: hdr from cal on.
hdrOff: hdr from cal off.
RETURNS:
istat: 1 it is a valid cal onoff pair
: -1 The hdrOn passed in is not a cal on
: -2 The hdrOff passed in is not a cal off
DESCRIPTION:
Check that the two headers passed in belong to a calonoff pair.
The routine checks the procedure name and that car[*,0] contains 'on' and
'off' respectively.
EXAMPLE:
input 2 records and check that they are from a cal on,off
print,posscan(lun,32500007L,1)
print,corgetm(lun,2,bb)
istat=chkcalonoff(bb[0].b1.h,bb[1].b1.h)
(See /pkg/rsi/local/libao/phil/Cor2/chkcalonoff.pro)
NAME:
coraccum - accumulate a record in a summary rec
SYNTAX: coraccum,binp,badd,new=new,scl=scl,array=array,brd=brd
ARGS: binp[]: {corget} input data
badd: {corget} accumulate data here
KEYWORDS:
new : if keyword set then this is the first call, alloc badd.
scl[] : float if provided, then scale the binp data by scl before
adding. This can be used to weight data by g/t.
if scl is an array it should equal the number of output
boards requested (see brd keyword). If it is a scalar
then that value will be used for all the boards.
array : if set then treat badd as an array and add element wise
the array binp[n] to badd[n]. Use with the new keyword.
By default if badd is an array then binp will be added
element wise.
brd[] : int if specified then only process the specified boards.
numbering is 1,2,3,4.
DESCRIPTION:
Accumulate 1 or more records worth of data into badd. If keyword
/new is set then allocate badd before adding. The header values in
badd will come from the first record added into badd.
Each element of binp[i] will be added to badd with the following weights:
1. If scl= is not defined, it defaults to 1.
2. if binp[i].b1.accum is 0., it is set to 1.
3. if binp[i].b1.accum is < 0 it is multplied by -1.
(This happens after coravg has been called on an accumlated badd.
It divides by badd.b1.accum and then sets badd.b1.accum to its negative.)
4 sclFact= binp[i].b1.d*scl*binp[i].b1.accum
5. badd.b1.d+=sclFact*binp[i].b1.d
6. badd.b1.accum+=sclFact
When corplot is called, it will scale the data by 1./badd.b1.accum. before
plotting.
When calling coraccum with the new keyword, you can include the
array keyword. This will allocate badd to be the same dimension as
binp. All future calls using badd will add binp element wise to badd.
This can be used when accumulating multiple maps.
Accumulated data must be of the same type (numlags, numbsbc, bw,etc..).
If you have observations to accumulate with only partial overlap of the
data types, you can use the brd keyword to specify which boards to accum. The
accumlated data must still be of the same type.
Example:
print,corget(lun,b)
coraccum,b,badd,/new
print,corget(lun,b)
coraccum,b,badd
corplot,badd
; Add n scans together element wise:
for i=0,n-1 do begin
print,corinpscan(lun,b)
coraccum,b,bsum,new=(i eq 0),/array
endfor
;
; input an entire scan and then plot the average of the records
; (this can also be done directly by corinpscan).
print,corinpscan(lun,b,scan=scan)
coraccum,b,bsum,/new
corplot,bsum
;
; let scan 1 be:4 brds,2pol, 1024 lags
; let scan 2 be:2 brds,2 pol,512 lags followed by 2 brds,2 pol, 1024 lags.
; To accumulate brds 1,2 of scan 1 with brds 3,4 of scan 2:
print,corinpscan(lun,b,scan=scan1)
coraccum,b,badd,/new,brd=[1,2]
print,corinpscan(lun,b,scan=scan2)
coraccum,b,badd,brd=[3,4]
help,badd.b1,/st
** Structure <3957c8>, 4 tags, length=9200, refs=2:
H STRUCT -> HDR Array[1]
P INT Array[2]
ACCUM DOUBLE 2.00000
D FLOAT Array[1024, 2]
(See /pkg/rsi/local/libao/phil/Cor2/coraccum.pro)
NAME:
coracf - compute the acf from the spectra
SYNTAX: bacf=coracf(b,/spectra)
ARGS:
b[n] : {corget} spectra to convert
KEYWORDS:
spectra: if set, then the input is acf's and they want spectra.
RETURNS:
bacf[n]: {corget} spectra converted to acf
DESCRIPTION:
Convert the spectra to acf's.
(See /pkg/rsi/local/libao/phil/Cor2/coracf.pro)
NAME:
corallocstr - allocate an array of structures.
SYNTAX: barr=corstostr(b,len)
ARGS: b: {corget} structure to duplicate
len: long length of the array to allocate
RETURNS:
barr[n]: {corget} array of structures to return.
DESCRIPTION:
Use corallocstr to allocate an array of {corget} structures. This routine
is necessary since each corget structure returned is an anonymous structure;
(See /pkg/rsi/local/libao/phil/Cor2/corallocstr.pro)
NAME:
coravg - average correlator data.
SYNTAX: bavg=coravg(binp,pol=pol,max=max,min=min,norm=norm,onlypol=onlypol)
ARGS: binp[]: {corget} input data
KEYWORDS:
pol : if set then average polarizations together
onlypol : only average polarizations. if binp[] is an array
do not average together records.
max : return max of each channel rather than average
min : return min of each channel rather than average
norm : if set then normalize each returned sbc to a mean of 1
(used when making bandpass correction spectra).
RETURNS:
bavg : {corget} averaged data.
DESCRIPTION:
coravg() will average correlator data. It has 3 basic functions:
1. average multiple records in a array together.
2. compute the average of an accumulated record (output of coraccum)
After doing this, it will set the bavg.b1.accum value to be the
negative of what was there (so corplot does not get confused).
3. average polarizations. This can be done on a single record, or
data from steps 1 or 2 above. If the /onlypol keyword is set then
it will average polarizations but it will not average together records.
The data is returned in bavg.
Polarization averaging will average the two polarizations on all boards
that have two polarization data. If polarizations are on separate correlator
boards then the routine will average all boards that have the same setup:
nlags, freq, and bandwidth. It will not combine boards that have two
polarizations per board.
If polarization averaging is used then the following header fields will
be modified.
b.(i).h.std.grptotrecs to the number of boards returned
b.(i).h.std.grpcurrec renumber 1..grptotrecs
b.(i).h.cor.numsbcout 2->1 for dual pol sbc
b.(i).h.cor.numbrdsused to the number of returned brds after averaging.
b.(i).h.cor.numsbcout to the number of returned sbc in each brd after
averaging
b.(i).h.dop.freqoffsets in case we reduce the number of boars from
pol averaging.
If binp contains multiple records then the max or min keyword can be
used to return the max or minimum of each channel. This will not work
on data from coraccum since the sum has already been performed.
Example:
; average polarizations:
print,corget(lun,b)
bavg=coravg(b,/pol)
;
; average a scan and then average polarizations
;
print,corinpscan(lun,b)
bavg=coravg(b,/pol)
;
; average the data accumulated by coraccum
; then average polarizations
;
print,corinpscan(lun,b)
coraccum,b,baccum,/new
print,corinpscan(lun,b)
coraccum,b,baccum
bavg=coravg(baccum,/pol)
;
; return max value from each channel
;
print,corinpscan(lun,b)
bmax=coravg(b,/max)
; average polarizations but do not average records. An
; example would be where pfposonoff() processes all of the
; on/off pairs in a file (possibly different sources), and you
; want to average the pols of each source.
n=pfposonoff(file,bar,tsys,cals) ; bar[n]
baravg=coravg(bar,/onlypol) ; baravg[n] with pols averaged
(See /pkg/rsi/local/libao/phil/Cor2/coravg.pro)
NAME:
coravgint - average multiple integrations.
SYNTAX: bavg=coravgint(bin)
ARGS: bin[] : {corget} input data
RETURNS:bavg : {corget} averaged data
DESCRIPTION:
If you have input an entire scans worth of data (eg corinpscan), you can
use coravgint to compute the average of all of the records of the scan. It
returns a {corget} structure that is the average of the N records input.
The header returned will be from the first record of the data (bin[0]).
Example:
assume scan 32500029 has 300 1 second integrations. Then:
print,corinpscan(lun,bin,scan=32500029L)
bavg=coravgint(bin)
bavg will then contain the average of the 300 records.
NOTES:
Many routines will do the averaging for you automatically if you request
it (eg corinpscan,coronoffpos,..)
SEE ALSO:
cormedian - to compute the median rather than the average of a scan.
(See /pkg/rsi/local/libao/phil/Cor2/coravgint.pro)
NAME:
corbl - baseline a correlator data set
SYNTAX: istat=corbl(bdat,blfit,maskused,deg=deg,mask=mask,$
edgefract=edgefract,auto=auto,sub=sub,svd=svd)
ARGS:
bdat[n]:{corget} input data to baseline. n can be >=1
KEYWORDS :
deg: int user can pass in polynomial degree for fit.
mask: {cormask} user can pass in cormask structure to use for mask
(see cormask()).
edgefract[n]: float The fraction of each edge of the bandpass to not use
in the fit. if n=1 then use the same fraction for both
sides. If n=2 then use [0] for the left edge and [1] for
the right edge. If both mask= and edgefr= are specified
then mask= will be used.
auto: if set then automatically do the baselining from the
user's input data, don't query the user for the parameters.
.. don't forget to set edgefract or pass in mask
sub : if set then return the data-blfit rather than the fit.
svd : if set then use svdfit rather than poly_fit for the
fitting. This works better for higher order polynomials
but it is slower.
RETURNS :
blfit[n]:{corget} the baselined fits or bdat-bfit if /sub set.
maskused:{cormask} the mask used for the fits.
istat : int 1 if all the fits done, 0 if one or more sbc not fit.
DESCRIPTION:
corbl will baseline a correlator data set. bdat can be a single
corget structure or an array of corget structures. By default the routine
will interactively baseline the spectra for each correlator board. The user
supplies the board, mask, and fit order from a menu.
On exit, the routine will return the fit values in a corget structure.
If the /sub keyword is set then the difference bdat-bfit will be returned.
The mask that was used will be returned in maskused. If bdat is an array
then the averages will be used for the plotting (data and fits) but the
fits will be done separately for each record.
The auto keyword will do the fit without prompting the user. In this
case you must enter the keywords:
deg= as the polynomial degree for fit. It will use the same for all boards.
You must also specify the mask to use by passing in mask or setting
edgefract to the fraction of the bandpass to ignore on each edge (the same
value will be used for all boards).
For interactive use the program menu is:
KEY ARGS FUNCTION
m .. define mask
c b .. copy mask from board b to current brd
f n .. fit polynomial degree n
h h1 h2 .. change horizontal scale for plot to h1,h2
v v1 v2 .. change vertical scale for plot to v1,v2
b brd .. switch to board..1->nboards
q .. quit
current board: 1
brdsLeftToDo: Brd1 Brd2 Brd3 Brd4
The plots are color coded as:
White : first polarization of this board
Red : 2nd polarization of this board (if there is one).
Green : fits to the data.
yellow: the mask that is used for the fit
An example of using the program is:
print,corbl(bdat,bfit,/usesvd)
1. Adjust the horizontal and vertical range of the plot with:
h h1 h2
v v1 v2
2. define the mask for the current board:
m
then use the left button to start,stop a range and the left button to
quit.
3. try fitting various orders to the data and mask:
f 3
when the fit has been done, an extra line is output giving stats of the fit:
FitInfo deg: 3 usesvd:1 mask%: 75.7 Rms: 0.0602 0.0581
We used a 3rd order fit, the mask is 75.7% of the entire range, and the
residual rms's are .06 and .058 in whatever units your spectrum is.
Trying a 7th order fit would give:
f 7
FitInfo deg: 7 usesvd:1 mask%: 75.7 Rms: 0.0461 0.0455
Here the residuals have decreased to .046 .
4. Move to another board:
b 3
5. define the mask for this board then fit.
6. repeat for all boards. To exit enter q
If you do not used the /sub keyword, then the fits are returned. The
baseline can then be removed using:
bldat=cormath(bdat,blfit,/sub)
Idl has 2 polynomial fitting routines: poly_fit that uses matrix inversion
and svdfit that uses singular value decomposition. poly_fit is faster while
svdfit is more robust for larger order fits. The fits are done in double
precision using the xrange 0. to 1.
SEE ALSO: cormask, blmask (in general routines),cormath
(See /pkg/rsi/local/libao/phil/Cor2/corbl.pro)
NAME:
corblauto - automatic baselining of a correlator data set
SYNTAX: istat=corblauto(bdat,blfit,maskused,coef,$
edgefract=edgefract,masktouse=masktouse,$
nsig=nsig,ndel=ndel,maxloop=maxloop,$
deg=deg,sub=sub,verbose=verbose,fsin=fsin,
double=double,raw=raw,plver=plver,badfit=badfit
ARGS:
bdat:{corget} input data to baseline.
KEYWORDS :
edgefract[n]: float The fraction of each edge of the bandpass to not use
in the computations. if n=1 then use the same fraction for
both sides. If n=2 then use [0] for the left edge and [1] for
the right edge. The default value is .08.
The keyword masktouse will override this option.
masktouse:{cormask} if provided,then use this rather than edgefract for the
starting mask. This allows to you to force certain
areas to not be included in the fit. see cormask() for
creating the mask. If the raw keyword is used, then this
is just a float array instead of a cormask structure.
nsig: int The rms loop will throw out all outliers greater than
nsig*sigma. By default nsig is 3.
ndel: int The number of extra channels to remove whenever we
throw out a channel because of large rms. def=1
maxloop: int When removing outliers in the rms loop, it will continue
looping until no points are removed or maxloop iterations
is hit. By default maxloop is 15.
deg: int maximum degree of polynomial fit.
The routine starts at order 1 and iterates the
rmsloop,fit iterating the degree after each pass
until the last fitorder done is deg. By default deg is 1.
sub : if set then return the data-blfit rather than the fit.
verbose : 1 plot status after each fit, stop after last fit of plot
2 plot status each fit, stop after each plot
3 plot status each rms loop. stop each plot
-1 same as 1 but don't wait..
fsin : int the order of cos(nx),sin(nx) to also fit.
double : if set then fit using double precision
raw: if raw keyword set, then the input/output is a float
array instead of a correlator data struct. Raw should
be set to the number of elements the the bdata array.
The maskused argument will be returned as if 1 brd, 1 pol
was used. This assumes that the data in bdat are
evenly spaced.
plver[2]:float vertical range [min,max ]to use for plotting bandpasses
(before subtraction). the default is autoscaling.
badfit : int badfit=0 . if a fit fails, the routine returns without
doing the follwowing bands. this is the default
badfit=1 . just skip this band, continue with the other
bands
RETURNS :
blfit :{corget} the baselined fits or bdat-bfit if /sub set.
maskused:{cormask} the mask used for the fits.
coefInfo:{} structure holding the coef's for the fit of each sbc
istat : int 1 if all the fits done, 0 if one or more sbc not fit.
DESCRIPTION:
corblauto will automatically create a mask and then use it to baseline
a correlator data set. The input is a {corget} structure of data. It will
return the masks created, the coefs of the fits, and the fits. It will
optionally remove the fits from the input data and return the difference
data-fit. The verb= option will plot the data as the fitting is done.
The routine can be used for mask creation (just ignore the fits)
or mask creation and baselining.
The /raw keyword lets you input an array of floats rather than a
{corget} structure. It will process the data and return the mask,fits
as if there was one correlator board with 1 polarization.
The parameters that control the algorithm are:
edgefract: This is the fraction of points on each edge of the data that
are not used. The default is .08 . For a sbc 1024 points this
would ignore 82 points on each edge.
masktouse: This is a mask with 0's and 1s. The 0 parts will not be
used in the fits. See cormask() on how to create the mask.
deg: The maximum degree for the fit. The default is 1. Values larger
than 12 or 13 may have pivot errors in the fitting.
ndel: The number of extra channels to throw out on the left and right
of each new channel we delete. This insures that skirts of
lines don't get included in the mask.
nsig : The rms computation removes all points greater than nsig*sigma.
By default nsig is set to 3.
maxLoop: The rms removal loops until there are no points within Nsigma*sigma
or maxloop iterations have occured. By default maxLoop is 15.
---------------------------------------------------------------------------
ALOGORITHM:
0. compute the starting points to use from the edge fraction or the masktouse.
The excluded points will never be included in the mask selection or fits.
let x0 by the x points, y0 be the y points after this selection.
1. set fitOrder=0, set yfit=mean(y0)
2. set yresiduals= y0 - yfit
the rms loop starts here
3. compute the rms of yresiduals
4. remove all points in y with yresiduals gt nsig*sigma .
For each point removed, also remove ndel adjacent points (determined
by the ndel keyword).
5. if the number of points removed in 4 is gt 0 and we've looped lt maxloop
times then goto 3.
fitting .. we now have a set of points within nsig of the mean.
6. fitdeg=fitdeg+1
7. fit a polynomial of fitDeg and a harmonic function of order fsin
to the points left in y.
8. evaluate the fitted function for the entire dataset y0
yfit=poly(x0,coef).
9. If fitdeg lt keyword deg goto 2.
---------------------------------------------------------------------------
The coef structure returns the fit results:
{corcoef}:
pol[2,4]: int 1-> this sbc/pol is used, 0--> it isn't
deg : int of polynomial fit
fsin : int of harmonic fit
nlags[4]: int number lags each sbc
coefAr[(deg+1)+fsin*2,2,4]: float coef for each fit. the polynomial
coef are followed by the cos,sin coef.
rms[2,4]: float of the fitted region within the mask
maskFract[2,4]: float npointsmask/nlags
To use the coefficients to recompute the fit, the x values should go
from -pi to +pi If Nchn is the number of channels in the spectra then
x=((findgen(Nchn)/(Nchn-1.) - .5)*2.*pi
EXAMPLES:
1. Use a polynomial fit of 3, remove the fit from the data, plot the
data and stop after each sbc done:
istat=corblauto(bdat,bfit,maskused,coefinfo,deg=3,/sub,verb=1)
; plot the datafit with the mask overplotted
corplot,bfit,cmask=maskused
2. Use a polynomial fit of 2 and a harmonic fit of 2. Create a mask
to use before calling this. Plot the results without stopping for keyboard
input
cormask,bdat,mask
istat=corblauto(bdat,bfit,maskused,coefinfo,deg=2,fsin=2,/sub,$
masktouse=masktouse,verb=-1)
On return :
coefinfo.coefar[0:2, are the c0,c1,c2 of the polynomial fit and
coefinfo.coefar[3:4] are the cos(x) sin(x) amplitudes while
coefinfo.coefar[5:6] are the cos(2x) sin(2x) amplitudes
SEE ALSO: cormask, corbl
(See /pkg/rsi/local/libao/phil/Cor2/corblauto.pro)
NAME:
corblautoeval - evaluate corblauto fit
SYNTAX: yfit=corblautoeval(x,coef,sbc=sbc,pol=pol)
ARGS:
x[n]: float xvalue to fit . Full scale should be 0. to 1.
coefI: {} coef struct returned by corblauto
KEYWORDS:
sbc: int sbc to eval. default is 1. count 1..n
pol int pol to eval. default is 1. count 1,2
RETURNS:
yfit[n]: float the fit evaluated at x
DESCRIPTION:
corblauto fits a polynomial and harmonic function returning the
fit coeficients in coef. This routine will evaluate the fit at the
specified points x[n]. x should go 0 to 1. for full scale.
SEE ALSO: corblauto
(See /pkg/rsi/local/libao/phil/Cor2/corblautoeval.pro)
NAME:
corcalcmp - given calon,off buffers, compute scale to kelvins
SYNTAX: istat=corcalcmp(calOn,calOff,calScl,datToScl,han=han,$
useinpm=useinpm,bla=bla,blda=blda,edgefract=edgefract,$
mask=mask,caldif=caldif,calValAr=calValAr,Tsys=Tsys
ARGS:
calOn[n]: {corget} cal on data. If n is greater than 1, then
the n calOn recs will be averaged together.
calOff[n]: {corget} cal off data. If n is greater than 1, then
the n caloff recs will be averaged together.
KEYWORDS:
han: if set then hanning smoooth the cal data
useinpm: if set then user passed mask in via mask keyword.
bla: If set then corblauto will use the caldeflection
to compute the mask of channels to use
blda: If set then corblauto will be passed
(calon-caloff)/caloff . It will fit and compute
the mask using this value.
edgefract[2]:float fraction of the bandpass to ignore on each
edge. Default is .08 . If 1 number is provided,
it is will be used for both edges.
extraParms: any extra parameters will be passed to
corblauto (if /bla is set). In particular you
can used deg=,fsin=,verb= .
RETURNS:
istat: int 1 ok, 0 trouble getting cal value
calScl[2,numbrds]: float factor that converts correlator counts to
kelvins for each pol of each board [pol,brd].
[0,n] is the first pol of each board (whether
it is polA or polB). If a board has only 1 pol
then calScl[1,n] will be 0.
See WARNIG below..
datToScl[m]: {corget} If provided, this data will be scaled to
kelvins for you.
mask: {cormask} The mask used when computing the cal. This is
computed in the following order:
1. if bla set, then mask comes from corblauto
2. if mask is provided then it is used as
passed in:
3. if edge= provided, then mask is computed
ignoring the edgefraction specified.
4. edgefraction is set to .08 and then use 3.
keyword.
caldif: {corget} return calon-caloff in correlator counts.
calValAr: fltar(2,nbrds) return calvalues for each board. index
[0,*] is the first spectra of board, [1,*] is the
2nd spectra of the board.
Tsys: fltar(2,nbrds,2) return Tsys as measured from the spectra
[pol,nbrds,calOn/calOff]
DESCRIPTION:
Given the cal on and off records, compute the scaling factors to go from
correlator counts to kelvins. It returns the results in calScl. If the
argument, datToScl is provided, then scale this data to kelvins using
the computed values (it should be {corget} structures).
The routine computes (where .i are the individual channels):
calDif.i= calon.i - calOff.i
It will then average the cal over a section (or mask) of the bandpass.
The mask is derived in the following order:
1. /useinpm is set. The user passes in the mask
2. /bla is set. corblauto will determine the mask from the
caldeflection.
3. edgefract is provided. It will make a mask ignoring this fraction
of channels from each edge.
It then computes (for each pol on each brd):
calScl[sbc,brd]= calK/mean(calDefl.i )
where the mean is summed over the non-zero mask elements.
If datToScl is provided then it will be converted to kelvins and returned.
The mask used in the cal computation is returned (if you do a bandpass
correction you will probably want to normalize to this set of
channels).
EXAMPLE:
This routine can be used if the cal onoff was taken with a
non stanard routine. Suppose you have a cal Off followed by a
cal on record:
print,corgetm(lun,b,2,/han) ; assume b[0] is cal off ,b[1] is cal on
istat=corcalcmp(b[1],b[0],calscl,b,/bla,verb=-1,deg=4,fsin=4,mask=mask)
Notes:
Calon,caloff are assumed to be from the same calon/off measurement.
It uses the first entry of calOn to find the rcvr, freq to get the
cal values.
datToScl should be the same setup as calOn,calOff (eg same number
of boards, and pols per board.
WARNING:
calscl is returned as a fltarr(2,nbrds) where the first index
is the sbc per board. If the board has only 1 sbc, then
calscl[1,x] is 0. Be careful passing calscl to the cormath() routine
to do the multiplication. Cormath wants 1 scaler for each sbc. It does
not want the extra zeros.
the WRONG way..
bk=cormath(b,smul=calscl) ... This will fail if a brd has 1 sbc.
the CORRECT way..
calsclOk=calscl[where(calscl ne 0.)] .. gets rid of the zeros
bk=cormath(b,smul=calsclok) .. this will work ok..
(See /pkg/rsi/local/libao/phil/Cor2/corcalcmp.pro)
NAME:
corcalib - intensity calibrate a single spectra
SYNTAX: istat=corcalib(lun,bdat,bcal,bfit,scan=scan,calscan=calscan,$
calinp=calinp,datinp=datinp,avg=avg,maxrecs=maxrecs,$
han=han,sl=sl,edgefract=edgefract,mask=mask,$
bpc=bpc,fitbpc=fitbpc,smobpc=smobpc,blrem=blrem,svd=svd,nochk=nochk
ARGS:
lun : int file descriptor to read from
KEYWORDS
scan : long scan number for data. default is current position
calscan : long scan number for cal on scan. def: scan following
the data scan.
calinp[2]:{corget} pass in the caldata rather than read it from lun
datinp[n]:{corget} pass in the data rather than read it from lun
avg: if set then return the averaged source record
maxrecs: long maximum number of recs to read in. default is 300.
han: if set then hanning smooth the data.
sl[]: {sl} array used to do direct access to scan.
edgefract[1/2]: float fraction of bandpass on each side to not use during
calibration. default .1
mask:{cormask} mask structure created for each brd via cormask routine
use this rather than edgefract.
note: currently mask.b1[1024,2] will use the mask for the
first pol of each brd for both entries of the board.
bpc: int 1 band pass correct with cal off
2 band pass correct with calon-caloff
3 band pass correct (smooth or fit) with data spectra
The default is no bandpass correction
fitbpc: int fit a polynomial of order fitbpc to the masked
version of the band pass correction and use the
polynomial rather than the data to do the bandpass
correction. This is only done if bpc is specified.
smobpc: int smooth the bandpass correction spectra by smobpc channels
before using it to correct the data. The number should be
an odd number of channels. This is only done if bpc is
specified.
blrem: int Remove a polynomial baseline of order blrem. Fit to the
masked portion of the spectra. This is done after
any bandpass correction or averaging.
svd : If baselining is done (blrem) then the default fitting
routine is poly_fit (matrix inversion). Setting svd
will use svdfit (single value decomposition) which is
more robust but slower.
nochk : if set then don't bother to check the cal records
to see if they are valid (in case they were written
with a non-standard program.
RETURNS:
bdat: {corget} intensity calibrated data spectra
bcal: {corget} intensity calibrated cal spectra
bfit: {corget} if supplied then return the smoothed or fit data that was
used for the bandpass correction.
istat: 1 ok
: 0 hiteof
:-1 no cal onoff recs
:-2 could not get cal value
:-3 cal,data scans different configs
:-5 sbc length does not match mask length
DESCRIPTION:
For those people who do not do position switching, corcalib allows you
to scale spectra from a src only scan to kelvins and optionally bandpass
correct it. The routine uses a src scan and a cal on,off pair. The data can
be read from disc or input directly to this routine (via calinp, datinp).
On output bdat and bcal will be in Kelvins. By default the individual records
of the scan are returned. If the /avg keyword is used, the average of all
of the src records will be returned. If the bfit argument is supplied then
the fit or smoothed version used for the bandpass correction will also be
returned. It will be scaled to the median bdat value so you can overplot
them.
If the data is input from disc then lun should be the file descriptor
from the open command. By default it will start reading the src scan from
the current file position. The scan=scan keyword lets you position to the
src scan before reading. By default the calscans will be the two scans
following the src scan. The calscan=calscan keyword lets you position
to the cal on scan before reading them. If the scans on disc have more than
300 records you need to use the maxrecs= keywords so the entire scan will
be input.
By default 10% of the bandpass on each edge is not used for the calibration
( when computing the conversion factor:
CalInKelvins/ <Calon[maskFrqChn]-calOff[maskFrqChn]>).
You can increase or decrease this with the edgefract keyword. The mask
keyword overrides the edgefract keyword and allows you to use a mask for
each sbc (use cormask to create the mask before calling corcalib).
The calibration will then only use the channels within the mask when computing
the gain calibration factors. This mask can be used to exclude rfi or
spectral lines.
Bandpass correction can be done with the cal off scan, the
calon-caloff difference spectrum, the data scan, or not at all. These
can be divided into the data scan as they are (although dividing the
data scan into itself is not very interesting!) or you can smooth or
fit a polynomial to the bandpass correction spectrum and then use the
fit/smooothed spectrum for the bandpass correction.
The keyword fitbpc=n will fit an n'th order polynomial to the bandpass
selected by the bpc keyword. Only the area within the mask (or edgefraction)
will be used for the fit.
The keyword smobpc=n will smooth the bandpass selected by the keyword
bpc and use it to do the bandpass correction (n should be an odd number >= 3).
You can pass in the data and/or calscans directly by using the
datinp, calinp keywords.
THE PROCESSING:
let Src be the src spectral data
let CalOn be the calOn spectra
let CalOff be the calOff spectra
let < > average over selected channels.The names will then have
Avg appended to them (eg calOnAvg=<calOn>)
let IndFrq be the set of frequency channels selected to use for the
calibration
The calibration consists of:
1. choose indFrq (the channels to use) in the following order:
a. The mask for each board from the mask keyword
b. Use edgefract to throw out this fraction of channels at each edge.
c. Use an edgefraction of .1 (10%)
2. compute CalOnAvg=<calOn[IndFrq]>,calOffAvg=<calOffAvg[IndFrq]>
3. Scale to Kelvins using: CntToK=CalValK/(calOnAvg-calOffAvg)
CalOnK =calOn*CntToK
CalOffK=calOff*CntToK
SrcK =Src*CntToK
4. If band pass correction is selected (bpc=1 or 2) then:
bpc=1: bpcN= calOff/<calOff[IndFrq]>
bpc=2: dif= calOn - calOff
bpcN= (dif)/<dif[IndFrq]>
If fitBpc > 0 then
bpcN=polyfit order fitbpc to bpcN[IndFrq]
else if smobpc > 2 then
bpcN=boxcar smooth (bcpN, smooth=smobpc)
then
SrcK=SrcK/bpcN
When deciding on the mask or edge fraction to use, you should have
a region where the calon-calOff is relatively flat (no filter rolloff and
no rfi).
EXAMPLE:
Suppose we have the following scans:
corlist,lun
SOURCE SCAN GRPS PROCEDURE STEP LST RCV
B1641+173 210200235 5 on on 17:42:08 5
B1641+173 210200236 1 calonoff on 17:47:10 5
B1641+173 210200237 1 calonoff off 17:47:21 5
40.42+0.70 210200238 5 on on 18:02:53 5
40.42+0.70 210200239 1 calonoff on 18:07:56 5
40.42+0.70 210200240 1 calonoff off 18:08:07 5
To process the first two sets:
--rew,lun
--print,corcalib(lun,bdat,bcal); will process the first set
--print,corcalib(lun,bdat,bcal,/han); will process the 2nd set with hanning
To process the 2nd set directly with an edgefraction=.12:
--print,corcalib(lun,bdat,bcal,scan=210200238L,edgefract=.12)
To input the data first, interactively create a mask, and then process
the data with a mask
--print,corinpscan(lun,bdatinp,scan=210200238L,/han)
--print,corgetm(lun,2 ,bcalinp,/han) ; reads the next two records
--cormask,bcalinp,mask ; interactively creates the mask
--print,corcalib(lun,bdat,bcal,calinp=bcalinp,datinp=bdatinp,mask=mask)
Use the same cal for multiple data scans:
--print,corgetm(lun,2 ,bcalinp,scan=210200236L/han);
--print,corcalib(lun,bdat1,bcal1,calinp=bcalinp,scan=210200235L)
--print,corcalib(lun,bdat2,bcal2,calinp=bcalinp,scan=210200238L)
Use the cal off for the bandpass correction. Use a 3rd order polynomial fit
to the cal off for the bandpass correction.
--print,corcalib(lun,bdat,bcal,scan=210200238L,bpc=1,fitbpc=3)
The bandpass correction is a bit tricky and depends on what type of
observations you are doing. The integration time for the off is usually
a lot less than the on positions so you need to use either the bandpass
fit or smoothing. It would probably be a good idea to add an option for
the user to input a bandpass to use for the correction (from an off src
position).
SEE ALSO:
cormask,corbl,corlist
(See /pkg/rsi/local/libao/phil/Cor2/corcalib.pro)
NAME:
corcalonoff - process a cal on/off pair
SYNTAX: istat=corcalonoff(lun,retDat,calOnrec,scan=scan,calval=calval,
sl=sl,swappol=swappol)
ARGS:
lun: unit number to read from (already opened).
retData[nbrds]: {cals} return data as array of structures.
calOnRec: {corget}optional parameter. If present it is the first
calon rec (in case you already read it).
KEYWORDS:
scan: if set, position to start of this scan before read
calval[2,nbrds]: cal value polA,polB for each board
sl[]: {sl} scan list array returned from getsl. If sl is
included then direct access positioning is available.
swappol: if set then swap polA,polB cal values. This can be
used to correct for the 1320 hipass cable reversal or
the use of a xfer switch in the iflo.
RETURNS:
istat - 1 ok, 0 did not complete
DESCRIPTION:
Read the calOnOff records and the cal values for the receiver
and frequency. Return the cal info in the retData[nbrds] array of
cal structs (one entry in retData for each correlator board).
The stucture format is:
retDat[i].h: {hdr} from first cal on
retDat[i].calval[2]:float calval in Kelvins for each sbc of this board
retDat[i].calscl[2]:float scale factor to use when scaling from
correlator counts to Kelvins: calT/(calon-caloff)
h.cor.calOn will contain the total power calOn (correlator units)
h.cor.calOff will contain the total power calOff (correlator units)
calon-caloff is the total power of the cal.
If a particular board has only 1 sbc, then the data will be in the
first entry [0] (whether or not it is polA or polB).
The routine will average multiple dumps if necessary.
You can scale your spectra to Kelvins via
b.bN.d[*,sbcn]= retDat[N-1].calscl[sbcn]*b.bN.d[*,sbcn]
(See /pkg/rsi/local/libao/phil/Cor2/corcalonoff.pro)
NAME:
corcalonoffm - process a cal on/off pair with a mask
SYNTAX: istat=corcalonoffm(lun,m,retDat,spc,scan=scan,calval=caltouse,sl=sl,$
edgefract=edgefract,han=han,swappol=swappol)
ARGS:
lun: unit number to read from (already opened).
m: {cormask} holds masks for the data 0 or 1.
It is created by the cormask routine.
Note. if m.b1[1024,2] has two entries per board, it uses
the first entry for both pols on the board.
retData[nbrds]: {cals} return data as array of structures.
spc[2]: {corget} if provided then return the calon (spc[0]) and
caloff spectra (spc[1]) after converting them to K.
KEYWORDS:
scan: if set, position to start of this scan before read
caltouse[2,nbrds]: cal value polA,polB for each board
sl[]: {sl} returned from getsl routine. If provided then
position to scan will be direct access.
edgefract: float if provided then create the mask ignoring m. Ignore
edgefract*lagsSbc lags from each side of the bandpass.
Return the mask in m (whatever is there gets overwritten).
han: if set then hanning smooth the spectra
swappol: if set then swap polA,polB cal values. This can be
used to correct for the 1320 hipass cable reversal or
the use of a xfer switch in the iflo.
RETURNS:
istat - 1 ok, 0 did not complete
DESCRIPTION:
Read the calOnOff records and the cal values for the receiver
and frequency. Compute the cal using the non zero lags in the mask m.
If edgefract is provided, then create the mask ignoring whatever is in m.
Zero edgefract*lagsbc lags from each side of the bandpass. Return this mask
in m (overwriting whatever was there).
The cal values and the scale factors (correlator counts to Kelvins) are
returned in the data structure retDat[n] (1 for each correlator board).
The data stucture format is:
retDat[i].h: {hdr} from first cal on
retDat[i].calval[2]:float calval in kelvin for each sbc of this board
retDat[i].calscl[2]:float to scale to kelvins calT/(calon-caloff)
h.cor.calOn will contain the total power calOn (correlator units)
h.cor.calOff will contain the total power calOff (correlator units)
calon-caloff is total power.
If a particular board has only 1 sbc, then the data will be in the
first entry [0] (whether or not it is polA or polB).
If the spc parameter is provided, the spectral data for the cal
on (spc[0]) and the cal off (spc[1]) will be returned after converting
them to Kelvins.
You can scale any other spectra to Kelvins using:
b.bN.d[*,sbcn]*= retDat[N-1].calscl[sbcn]
SEE ALSO:
corcalonoff
(See /pkg/rsi/local/libao/phil/Cor2/corcalonoffm.pro)
NAME:
corchkstr - check buffers for same structure
SYNTAX: istat=chbufstr(b1,b2)
ARGS:
b1: {corget} correlator data buf
b2: {corget} correlator data buf
RETURNS:
istat: 1 buffers are compatible
0 buffers are not compatible
DESCRIPTION:
Some routines store correlator data structures i{corget} in arrays
use corstostr. This is legal if the two buffers meet the following
requirements.
1. They have the same number of brds.
2. Each corresponding board has the same number of subcorrelators
3. Each sbc of each corresponding board has the same number of lags.
EXAMPLE:
istat=corget(lun,b1)
istat=corget(lun,b2)
same=corchkstr(b1,b2)
(See /pkg/rsi/local/libao/phil/Cor2/corchkstr.pro)
NAME:
corcmbsav - combine save files with multiple sources.
SYNTAX - npat=corcmbsav(savelistinp,b_all,b_sec,b_rms,b_info,usrcname,$
savefileout=savefileout,/filelist)
ARGS:
savelistinp[]: string list of save files to use as input
These files should contain:
bf[n]: an array of {corget} structures.
KEYWORDS:
filelist: if set then savelistinp is the data filelist used to input the
data. In this case this routine will generate the standard
save filenames to read:
corfile.ddmonyy.projid.n -> ddmonyy.projid.n.sav
It will assume that the .sav files are in the current directory.
savefileout: string name of save file to store the combined arrays in.
RETURNS:
npat : long number of patterns processed
b_all[]:{corget} the combined on/off -1 from the save files
b_sec[]:{corget} the combined secperchn from the save files
b_rms[]:{corget} the combined rms/mean from the save files
b_info[]:{} struct telling which kind of data is valid for
each b_rms/sec[m]
usrcname[l] string array holding the uniq source names present.
DESCRIPTION:
Process a set of idl save files each containing an array bf[n] of
corget structures. Combine these into a single array and return
them to the caller. optionally save them to an idl save file.
The routines pfposonoff(), pfposonoffrfi() process all of the
on/off patterns in a file (as long as they are the same data size).
They then save this info to an idl save file. This would typically be
1 days observations.
This routine's job is to combine all of the info from multiple
save files (days) concatentating them into 1 large array.
The input data in the save files are:
bf[n] {corget} holds the n on/off's processed for 1 file
arSecPerChn[n] {corget} secs/chan integration for each bf[n]. This is
output by pfposonoffrfi(). After rfi excision there may
not be the same integration time in each freq channel
arRmsByChn[n] {corget} rms/mean for each on/off after rfi excision
(actually after the rms fit, before the total power
test).
These will be combined into:
b_all[M] holds the combined bf[n]
b_sec[M] holds the combined arsecperchn[n]
b_rms[M] holds the combined arrmsbychn[n]
b_info[M] a struct holding source name for each entry and
whether or not there is b_sec and b_rms info for each
entry:
b_info[i].srcname ; source name
b_info[i].usesec ; 1 if it has secPerchn info from rfi processing
b_info[i].userms ; 1 if it had rms/mean info from rfi processing
The new arrays are returned to the caller. They can also be saved
to an idl save file with the savefileout keyword.
EXAMPLE:
Suppose project a1721 generated the following corfiles (either in
/share/olcor or /proj/a1721). You could process all the data with the
following code:
; The pfposonoff processing:
flist=['corfile.08may03.a1721.2',$
'corfile.08may03.a1721.3',$
'corfile.09may03.a1721.1']
dir=['/share/olcor/','/proj/a1721/'] ; directories to look for the files
voff=.03 ; plotting increment
han=1 ; hanning smooth
hor
ver
for i=0,n_elements(flist)-1 do begin &$
npairs=pfposonoff(flist[i],bf,tf,calsf,han=han,dir=dir,/scljy,/sav) &$
; plot the data so we can see what's up. offset each strip by voff
ver,-voff,(npairs+1)*voff &$
corplot,bf,off=voff &$
endfor
; The processed on/off-1 are now in the save files (one for each corfile)
; Use corcmbsav() to combine them:
; 1. input all the data from the daily save files.
; 2. save all of these arrays in the savbysrc filename.
; 3. the /filelist keyword flags that the string array has the
; original datafilenames and this routine should generate the
; input save file names.
;
savefileout='onoff_cmb.sav' ; place to store the new data
istat=corsavcmb(flist,b_all,b_sec,b_rms,b_info,$
savefileout=savefileout,/filelist)
;
WARNING:
This works as long as the data structure is the same type for
all of the on/offs. numboards, pol/brd, lags/pol
(See /pkg/rsi/local/libao/phil/Cor2/corcmbsav.pro)
NAME:
corcmpdist - compute distance (ra/dec) for scans
SYNTAX: dist=corcmpdist(ra,dec,b=b,slar=slar,hms=hms)
ARGS :
ra: float/double ra Hours J2000
dec: float/double dec Degrees J2000
KEYWORDS:
b[n]: {corget} array of correlator recs/scans to use for distance
sla[n]: {corsl} array of archive recs to use for distance
hms : if set then ra is in hms.s dec is in dms.s rather than hour,degrees
RETURNS:
dist[3,n]: float distance from ra,dec to each measurement.
[0,*] = ra - measured ra Arcminutes (greatcircle angle)
[1,*] = dec - measured dec arcminutes (angle)
[3,*] = total distance arcminutes (greatcircle angle)
DESCRIPTION:
Compute the distance from ra,dec to each measured entry. The measured
entries are passed in via b=b or slar=slar (you can use only one of these).
The requested ra,dec positions are used rather than the actual.
If b= is used then it is the requested position for each record.
If slar= is used then it is the average requested position of each scan.
The ra,dec arguments should also be in j2000 coordinates.
(See /pkg/rsi/local/libao/phil/Cor2/corcmpdist.pro)
NAME:
corcumfilter - cumfilter a correlator dataset.
SYNTAX: cmask=corcumfilter(b,limit,lengthfract,edgefract=edgefract)
ARGS:
b[n]: {corget} data to cumfilter
limit: float max deflection to allow
lengthfract:float fration of channels to use for filtering.
KEYWORDS:
edgefract: float fraction of channels to ignore on each side
(in case band pass still present)
RETURNS:
cmask[n]:{cormask} hold results of cumfiltering.
DESCRIPTION:
cumfilter is a routine written by carl heiles to remove outliers
from a dataset. The routine:
1. sorts the data by amplitude
2. selects length points (length=lengthfract*numpoints)
about the center of the sorted list
3. keeps all the points within limit*maxdeviation of the length
points from the center value.
The routine returns a cormask data structure that has a 0 for
points that were rejected and a 1 for points that were kept.
If b[n] is an array of structures then n cormask structures will
be returned (1 for each element of b).
EXAMPLE:
limit=2.
length=.5
print,corget(lun,b)
cmask=corcumfilter(b,rnage,length,edgefract=.08)
cmask is a single structure
You can then use cmask in basline fitting or operations:
istat=corbl(b,blfit,maskused,/auto,deg=1,mask=cmask,/sub)
print,corgetinpscan(lun,b)
cmask=corcumfilter(b,limit,length,edgefract=.08)
cmask is an array of cmask structures.
WARNING:
if b[] is an array , then corcumfilter will create and array of
masks. Some of the other cor routines (corplot,corstat, etc..) will
only accept single element cormasks so you need to be careful passing
any arrays of cmasks to other routines.
(See /pkg/rsi/local/libao/phil/Cor2/corcumfilter.pro)
NAME:
cordfbp - return the digital filter bandpasses.
SYNTAX: bp=cordfbp(b,force=force)
ARGS:
b: {corget} correlator data to make banpasses for.
KEYWORDS:
force: if set then force the recomputation of the bandpasses
rather than taking them from the common block if they are
already there.
RETURNS :
bp : {corget} holds the normalized digital filter bandpasses.
DESCRIPTION:
The normalized digital filter bandpasses are returned in a standard
{corget} data structure. If the input (b) is an array, then only
bandpasses for the first element of b are returned. The computed
bandpasses are stored in a local common block so that repeated calls with the
same type of data will go faster.
50 Mhz bandpasses are returned as all 1's (since there is no digital
filtering).
Stokes data returns the digital filter bandpasses for the
original PolA, polB data rather than I and Q. The last two polarized band
passes are returned as all 1's.
EXAMPLE:
; Input a data scan and divide each record by the digital filter bandpass.
; Note that the digital filter bandpasses are already normalized to unity.
;
istat=corinpscan(lun,b,scan=scan)
dfbp=cordfbp(b)
bpc=cormath(b,dfbp,/div)
The current filtering scheme seems to be:
filter on upconvert - all
filter 1 for interpolation after upconvert : double nyquist
filter 2 for interpolation after upnconvert: 12.5 Mhz and below
SEE ALSO: dfhbfilter, cormath
(See /pkg/rsi/local/libao/phil/Cor2/cordfbp.pro)
NAME:
corfindpat - get the indices for the start of a pattern
SYNTAX: nfound=corfindpat(sl,indar,pattype=pattype,rcv=rcv)
ARGS:
sl[]: {getsl} scan list array from getsl
KEYWORDS:
pattype: int Type of pattern we are looking for.
1 - on/off position switch with cal on/off
2 - on/off position switch whether or not cal there
3 - on followed by cal on ,off
4 - heiles calibrate scan two crosses
5 - heiles calibrate scan 1 or more crosses
6 - cal on,off
7 - x111auto with calonoff
8 - x111auto with or without cal
9 - heiles calibrate scan 4 crosses
If not provided than pattype 1 is the default.
dosl: lun if keyword dosl is set to a valid open file, then this
routine will do the sl=getsl(lun) call and return the
scan list in sl.
rcv : int if supplied then only find patterns for this receiver.
1..12 (see helpdt feeds) for a list of the feeds
RETURNS:
indar[npat]: long indices into sl[] for start of the pattern
npat : long number of patterns found
DESCRIPTION:
corfindpat() is used for the automatic processing of entire datafiles.
It processes a scanlist array (returned by getsl()) and returns an array of
pointers for the start of all the completed patterns of a particular type
that are located in the datafile. Patterns can be: on/off position switching,
heiles calibration runs, on with cal on/off,..etc.
The requirements for a completed scan depend on the pattern type:
type order needed requirements
1 posOn,posOff,calon,caloff. Number of on recs must equal number of
off recs
2 posOn,posOff. Number of on recs must equal number of
off recs.
3 posOn,calon,caloff.
4 calon,caloff,heiles calibrate scan with 2 crosses.
nrecs in each cross is the same.
5 calon,caloff,heiles calibrate scan with at least 1 cross.
nrecs in each cross is the same.
6 calon,caloff
7 x111auto (60recs) calon,caloff
8 x111auto (60recs)
9 calon,caloff,heiles calibrate scan with 4 crosses 120 samples/strip
nrecs in each cross is the same.
You can call sl=getsl(lun) once prior to this routine, or you can
set the keyword dosl=lun and it will create the sl array and return it.
It is also possible to limit the pattern to a particular receiver using
the rcv=rcvnum keyword.
EXAMPLE:
openr,lun,'/share/olcor/corfile.23aug02.x101.1',/get_lun
sl=getsl(lun)
; get poson,off,followed by cal on,off
npat=corfindpa(sl,indfound)
for i=0,npat-1 do begin
scan=sl[indfound[i]].scan
.. process this scan
; get ons followed by cal on,off. have the routine do the sl search.
; only get data for cband (rcv=9)
openr,lun,'/share/olcor/corfile.23aug02.x101.1',/get_lun
npat=corfindpa(sl,indfound,dosl=lun,pattype=3,rcv=9)
SEE ALSO: arch_gettbl,mmfindpattern,getsl (in general routines)
NOTE:
There is also a record type in the {sl} structure that lists the
name of the active pattern when the data was taken. It's coding and the
coding for this routine are not the same (sorry..). It may or may
not be accurate (for some test data, the pattern type is not set so the
last one supplied is used). Corfindpat differs in that it will
try and verify that the pattern is what it says it is and that it is
complete (eg calon followed by caloff).
(See /pkg/rsi/local/libao/phil/Cor2/corfindpat.pro)
NAME:
corfrq - compute the freq/vel array for a spectra
SYNTAX: retArr=corfrq(hdr,retvel=retvel,retrest=retrest)
ARGS:
hdr: header for this board
KEYWORDS:
retvel : if not equal to zero, then return velocity rather
than frequency array
retrest: if not equal to zero, then return the rest frequency
for the object being observered rather than the
topocentric frequency.
RETURNS:
retArr: Array of floating point frequencies or velocities.
DESCRIPTION:
Compute the topocentric frequency array (in Mhz) for the correlator board
corresponding to the header (hdr) passed in. If the keyword retvel
is set (not equal to zero) then return the velocity (optical definition).
If the keyword retrest is set, then return the rest frequency of the object
rather than the topocentric frequency.
The array returned (retArr) will have the same number of elements
as a single output sub correlator.
The order of the data assumes that the spectral channels are in
increasing frequency order (corget always returns the data this way).
If the spectra are spDat[2048] and then retAr[0] will be the lowest
frequecy or the highest velocity.
EXAMPLE:
.. assume 2 boards used, pola,b per board
corget,lun,b
frqArr=corfrq(b.b1.h)
frqArrRest=corfrq(b.b1.h,/rest)
velArr=corfrq(b.b1.h,/retvel)
plot,frqArr,b.b1.d[*,0]
(See /pkg/rsi/local/libao/phil/Cor2/corfrq.pro)
NAME:
corget - input next correlator record from disc
SYNTAX: istat=corget(lun,b,scan=scan,noscale=noscale,sl=sl,han=han)
ARGS:
lun: logical unit number to read from.
RETURNS:
b: structure holding the data input
istat: 1 ok
: 0 hiteof
:-1 i/o error..bad hdr, etc..
KEYWORDS:
noscale : if set, then do not scale each sub correlator to the
9level corrected 0 lag.
scan : if set, position to start of scan first..
han : if set, then hanning smooth the data
sl[] : {sl} array used to do direct access to scan.
This array is returned by the getsl procedure.
DESCRIPTION:
Read the next group of correlator data from the file pointed to
by lun. If keywords scan is present, position to scan before reading.
A group is the data from a single integration. Each correlator board
is written with a separate hdr and data array. The structure returned will
contain 1 to 4 elements depending on the number of boards being used.
b.b1
b.b2
b.b3
b.b4
each bN will have:
b.b3.h - complete header for this board
b.b3.p[2] int. 1-polA,2->polB, 0, no data this sbc
b.b3.accum double . accumulate scale factor (if used)
b.b3.d[nfreqchan,nsbc] the data. nsbc will be 1 or 2 depending on
how many sbc are using in this board.
use pol to determine what pol each sbc is. It will also tell you if
there is only 1 sbc pol[1] = 0. It will not compensate for
zeeman switching..
The header will contain:
.h.std - standard header
.h.cor - correlator portion of header
.h.pnt - pointing portion of header
.h.iflo- if,lo portion of header
.h.dop - doppler frequency/velocity portion of header
.h.proc- tcl procedure portion of header
The data is returned in increasing frequency order as floats.
If an i/o error occurs (hit eof) or the hdrid is incorrect (you are not
correctly positioned at the start of a header), then an error message
is output and the file is left positioned at the position on entry to the
routine.
EXAMPLE:
.. assume 2 boards used, pola,b per board (lagconfig 9)
istat=corget(lun,b)
b.b1.h - header first board
b.b2.d[*,0] - data from 2nd board, polA
SEE ALSO:
posscan,corgethdr
(See /pkg/rsi/local/libao/phil/Cor2/corget.pro)
NAME:
corgethdr - return the correlator header for the next group
SYNTAX: istat=corgethdr(lun,rethdr)
ARGS:
lun: int opened file to read
rethdr[nbrds]: {hdr} return array of headers brd 1.. thru brd n
(See /pkg/rsi/local/libao/phil/Cor2/corgethdr.pro)
NAME:
corgetm - input multiple correlator records to array.
SYNTAX: istat=corgetm(lun,ngrps,bb,scan=scan,sl=sl,han=han,noscale=noscale,$
check=check)
ARGS:
lun: int open file to read
KEYWORDS:
scan: long if provided, position to scan before reading.
sl[]: {sl} returned from getsl(). If provided then direct access
to scan is provided.
han: if set then hanning smooth
noscale: if set, then do not scale each sub correlator to the
9level corrected 0 lag.
check : if set, then verify that the records are compatible to store
in a single record.
RETURNS:
bb[]:{corget} array of corget structures that were read in
istat: 1 got all of the groups requested.
0 did not get all of the groups requested.
DESCRIPTION:
Read ngrp consecutive correlator records from disc and return them
in bb (an array of {corget} structures). The routine will read from the
current position on disc (after optional positioning to scan). It will
read across scan number boundaries if needed.
(See /pkg/rsi/local/libao/phil/Cor2/corgetm.pro)
NAME:
corhan - hanning smooth correlator data.
SYNTAX: corhan,b,bsmo
ARGS:
b[m]: {corget} data to hanning smooth
bsmo[m]: {corget} return smoothed data here. Note.. if only 1 argument
is passed to the routine then the data is smoothed in
place (it is returned in b).
DESCRIPTION:
corhan will hanning smooth the data in the structure b. If a single
argument is passed to the routine, then the smoothed data is returned
in place. If two arguments are passed (b,bsmo) then the data is returned
in the second argument.
EXAMPLE:
print,corget(lun,b)
corhan,b ; this smooths the data an returns it in b.
print,corget(lun,b)
corhan,b,bsmo ; this smooths the data an returns it in bsmo
(See /pkg/rsi/local/libao/phil/Cor2/corhan.pro)
NAME: corhcalrec- check if an input rec is a cal rec. SYNTAX: istat=corhcalrec(hdr) ARGS : hdr - header from 1 of the boards RETURNS: istat- 0 not a cal rec, 1-on, 2-off DESCRIPTION: Check if a record input is part of a cal record. EXAMPLE: corget(lun,b) istat=corhcalrec(b.b1.h)
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhcalval - return the pol A/B cal values for a sbc.
SYNTAX: stat=corhcalval(hdr,calval,date=date,swappol=swappol)
ARGS:
hdr: {hdr} header for board to check
KEYWORDS:
date[2]: intarray [year,dayNum] if provided, then compute the calvalues
that were valid at this epoch.
swappol: if set then swap the pola,polb cal values. This can
be used to correct for the 1320 hipass cable switch
problem or the use of a xfer switch in the iflo.
RETURNS:
calval[2]: float .. calValues in deg K for polA,polB
stat: int .. -1 error, 1 got the values ok.
DESCRIPTION:
Return the cal values in degrees K for the requested sbc. This
routine always returns 2 values (polA then polB) even if the header
is for a board that uses only one polarization.
The calvalues for the receiver in use are looked up and then the
values are interpolated to the observing frequency.
EXAMPLE:
input a correlator record and then get the calvalues for the
3rd correlator board:
print,corget(lun,b)
istat=corhcalval(b.b3.h,calval)
.. calval[2] now has the cal values in degrees K for polA and polB.
NOTE:
Some cals have measurements at a limited range of frequencies (in some
cases only 1 frequency). If the frequency is outside the range of measured
frequencies, then the closest measured calvalue is used (there is no
extrapolation).
The year daynum from the header is used to determine which set of
calvalue measurements to use (if the receiver has multiple timestamped
sets).
This routine computes the frequency of the sbc from hdr and then calls
calget().
SEE ALSO:
gen/calget gen/calval.pro, gen/calinpdata.pro
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhcfrrest - return the rest freq of band center.
SYNTAX: cfr=corhcfrrest(hdr)
ARGS:
hdr: {hdr} .. header for board to check
RETURNS:
cfr: double .. rest center freq of band in Mhz
DESCRIPTION:
Return the rf rest frequency for the band center of the requested board.
EXAMPLE:
get the rest freq of the 3nd board: cfrMhz=corhcfrrest(b.b3.h)
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhcfrtop - return the topocentric freq of band center.
SYNTAX: cfr=corhcfrtop(hdr)
ARGS:
hdr: {hdr} .. header for board to check
RETURNS:
cfr: double .. center freq of band in Mhz
DESCRIPTION:
Return the topocentric rf center of the band for the requested board.
hdr can be an array of hdrs but they must all come from the same
board. eg b[].b1.h
EXAMPLE:
get the freq of the 2nd board: cfrMhz=corhcfrtop(b.b2.h)
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhdnyquist - check if rec taken in double nyquist mode
SYNTAX: istat=corhdnyquist(hdr)
ARGS : hdr{} - header from 1 of the boards of the record
RETURNS: istat- 0 not taken in double nyquist mode.
1 taken in double nyquist mode.
DESCRIPTION:
Check if a record was taken in double nyquist mode.
EXAMPLE:
corget(lun,b)
istat=corhdnyquist(b.b1.h)
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhflipped - Obsolete..check if current data is flipped in freq.
SYNTAX: stat=corhflipped(corhdr)
ARGS:
corhdr[]: {corhdr} to check.
RETURNS: istat- 0 increasing frq, 1- decreasing freq order.
if corhdr[] is an array then istat will be an array of ints.
DESCRIPTION:
This routine has been replaced by corhflippedh(). corhflipped
was using the bit in the header that was not always being set correctly
for 1 ghz bandwidths.
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhflippedh - check if current data is flipped in freq.
SYNTAX: stat=corhflippedh(hdr,sbcnum)
ARGS:
hdr[]: {hdr} to check.
sbcnum:long sbc to check 1 thru 4.
RETURNS: istat- 0 increasing frq, 1- decreasing freq order.
if corhdr[] is an array then istat will be an array of ints.
DESCRIPTION:
Check if the correlator data for this sbc is stored in increasing
or decreasing frequency order (even or odd number of high side lo's).
This routine replaces Corhflipped. Corhflipped gave incorrect results for
1 ghz if's if the LO was below 1500 Mhz and it was being used as hi side, or
the lo was above 1500 Mhz and it was being used as a low side lo.
corhflippedh should work in all cases.
EXAMPLE:
sbcnum=1
check first board: istat=corhflippedh(b.b1.h,sbcnum)
(See /pkg/rsi/local/libao/phil/Cor2/corhflippedh.pro)
NAME:
corhgainget - return the gain given a header
SYNTAX: stat=corhgainget(hdr,gainval,date=date,az=az,za=za,onlyza=onlyza)
ARGS:
hdr[n]: {hdr} header for board to check
KEYWORDS:
date[2]: intarray [year,dayNum] if provided, then compute the gain value
at this epoch.
az[n]: fltarray If provided, use this for the azimuth value rather than
the header values
za[n]: fltarray If provided, use this for the zenith angle values rather
than the header values
onlyza: If set then return the za dependence (average of az)
RETURNS:
gainval: float .. gainvalue in K/Jy
stat: int -1 --> error, no data returned
0 --> requested freq is outside freq range of fits.
Return gain of the closed frequency.
1 --> frequency interpolated gain value returned.
DESCRIPTION:
Return the telescope gain value in K/Jy for the requested sbc.
The gain fits for the receiver in use are input and then the
values are interpolated to the az, za and observing frequency.
If hdr[] is an array then the following restrictions are:
1. each element must be from the same receiver and at the same
frequency (eg. all the records from a single scan).
2. If the az,za keywords are provided, they must be dimensioned the same
as hdr
EXAMPLE:
input a correlator record and then get the gain value for the
3rd correlator board:
print,corget(lun,b)
istat=corhgainget(b.b3.h,gain)
.. gain now has the gain value in K/Jy
input an entire scan and compute the gain for all
records of 1 sbc. assume 300 records..
print,corinpscan(lun,bar)
istat=corhgainget(bar.b3.h,gain)
gain is now a array of 300 records
NOTE:
Some receivers have measurements at a limited range of frequencies (in some
cases only 1 frequency). If the frequency is outside the range of measured
frequencies, then the closest measured gain is used (there is no
extrapolation in frequency).
The year daynum from the header is used to determine which set of
gain fits to use (if the receiver has multiple timestamped sets).
This routine takes the az,za, date, and frequency from the
header and then calls gainget().
If you input an array of corget recs , then they must all be from the
same sbc.
SEE ALSO:
gen/gainget gen/gaininpdata.pro
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhintrec - return integration time for a record
SYNTAX: secs=corhintrec(hdr)
ARGS : hdr{} - header from 1 of the boards of the record
RETURNS: secs:float seconds of integration for this record
DESCRIPTION:
Return the integration time for a record.
EXAMPLE:
corget(lun,b)
istat=corhintrec(b.b1.h)
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhsefdget - return the sefd given a header
SYNTAX: stat=corhsefdget(hdr,sefdval,date=date,az=az,za=za,onlyza=onlyza)
ARGS:
hdr[n]: {hdr} header for board to check
KEYWORDS:
date[2]: intarray [year,dayNum] if provided, then compute the gain value
at this epoch.
az[n]: fltarray If provided, use this for the azimuth value rather than
the header values
za[n]: fltarray If provided, use this for the zenith angle values rather
than the header values
onlyza: If set then return the za dependence (average of az)
RETURNS:
sefdval: float .. sefd value in Jy
stat: int -1 --> error, no data returned
0 --> requested freq is outside freq range of fits.
Return sefd of the closed frequency.
1 --> frequency interpolated sefd value returned.
DESCRIPTION:
Return the telescope sefd value in Jy for the requested sbc.
The sefd fits for the receiver in use are input and then the
values are interpolated to the az, za and observing frequency.
If hdr[] is an array then the following restrictions are:
1. each element must be from the same receiver and at the same
frequency (eg. all the records from a single scan).
2. If the az,za keywords are provided, they must be dimensioned the same
as hdr
EXAMPLE:
input a correlator record and then get the sefd value for the
3rd correlator board:
print,corget(lun,b)
istat=corhsefdget(b.b3.h,sefdval)
.. sefdval now has the sefd value in Jy
input an entire scan and compute the sefd for all
records of 1 sbc. assume 300 records..
print,corinpscan(lun,bar)
istat=corhsefdget(bar.b3.h,sefdAr)
sefdAr is now a array of 300 records
NOTE:
Some receivers have measurements at a limited range of frequencies (in some
cases only 1 frequency). If the frequency is outside the range of measured
frequencies, then the closest measured sefd is used (there is no
extrapolation in frequency).
The year daynum from the header is used to determine which set of
sefd fits to use (if the receiver has multiple timestamped sets).
This routine takes the az,za, date, and frequency from the
header and then calls sefdget().
If you input an array of corget recs , then they must all be from the
same sbc.
SEE ALSO:
gen/sefdget gen/sefdinpdata.pro
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhstate - decode status words for correlator header
SYNTAX: statInfo=corhstate(corhdr)
ARGS:
corhdr[]:{hdrcor}
RETURNS:
statInfo[]{corstate} decoded status structure
DESCRIPTION:
The correlator header contains various info that is encoded in bitmaps
(h.cor.state) This routine decodes this
bitmap and returns it in a structure. The input can be 1 or more
cor headers.
Example:
print,corget(lun,b)
corstate=corhstate(b.b1.h.cor)
help,corstate,/st ; to print it out
The returned structure format is:
a={corstate,$
dnyquist : 0 ,$; 1--> double nyquist
chiptest : 0 ,$; 1--> chip testmoded
blankOn : 0 ,$; 1--> radar blanking on
pwrcntInc : 0 ,$; 1--> power counter data included
relbitshift: 0 ,$; relative bit shift used
relbitshiftsgn: 0 ,$; 1-->relbitshift neg, 0--> relbitshift pos
rawacf : 0 ,$; 1--> raw acf (not bias removal
cmbacf : 0 ,$; 1--> corrected acfs
spc : 0 ,$; 1--> spectra
pack : 0 ,$; 1--> packed data (for decoder)
startImdSw: 0 ,$; 1--> start immediate software
startImdHw: 0 ,$; 1--> start immediate hardware
startTick1: 0 ,$; 1--> start 1 sec tick
startTick10: 0 ,$; 1--> start 10 sec tick
dmpDelTrgInt: 0 ,$; 1--> dump delay trig from interrupt
dmpDelTrgTick: 0 ,$; 1--> dump delay trig from tick
adjPwrStartScn: 0 ,$; 1--> adjpower start of scan (automatic)
adjPwrStartRec: 0 ,$; 1--> adjpower start of rece (automatic)
spcFlipped: 0 ,$; 1--> spectra is flipped on disc
isACor: 0 ,$; 1--> running as a correlator (0-->decoder)
totCntIncluded: 0 ,$; 1--> total counts included
calOff : 0 ,$; 1--> cal is off this rec
calOn : 0 ,$; 1--> cal is on this rec
complexDat: 0 ,$; 1--> complex data taken
level9: 0 ,$; 1--> 9 level data. 0--> 3 level
stokes: 0 ,$; 1--> stokes data
pwrCntI: 0 ,$; 1--> power count I included
pwrCntQ: 0 }; 1--> power count Q included
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corhstokes - check if record taken in stokes mode
SYNTAX: istat=corhstokes(hdr)
ARGS : hdr{} - header from 1 of the boards of the record
RETURNS: istat- 0 not taken in stokes mode. 1-->taken in stokes mode
DESCRIPTION:
Check if a record was taken stokes (polarization) mode.
EXAMPLE:
corget(lun,b)
istat=corhstokes(b.b1.h)
(See /pkg/rsi/local/libao/phil/Cor2/corhquery.pro)
NAME:
corimg - create an image of freq vs time for 1 sbc.
SYNTAX: d=corimg(b,brd,sbc,bscl=bscl,cmpscl=cmpscl,edgchn=edgchn)
ARGS:
b[n] : {corget} input data
brd : 1-4 which board to use
sbc : 1,2 in this board to use
KEYWORDS:
bscl : {corget} divide into each entry in b[n]. Usually comes from
coravgint().
cmpscl: int 1 compute mean ,use to flatten image (default)
2 compute median,use to flatten image
3 donot flatten image.
edgchn: int number of channels on each edge to set to unity
RETURNS:
d[lensbc,n]: floats .. scaled data
DESCRIPTION:
For the board and sbc of interest return a floating point array
d[nlags,nrecs]. It will be normalized by :
if bscl provided then
d=b/bscl for the particular sbc.
if cmpscl returns
0,1 d=b/avg(b)
2 d=b/median(b)
3 just return the data with no bandpass correction.
(See /pkg/rsi/local/libao/phil/Cor2/corimg.pro)
NAME:
corimgdisp - display a set of correlator records as an image
SYNTAX: img=corimgdisp(b,clip=clip,sbc=sbc,brdlist=brdlist,pol=pol,col=col,$
median=median,bpc=bpc,ravg=ravg,nobpc=nobpc,$
win=win,wxlen=wxlen,wylen=wylen,wxpos=wxpos,wypos=wypos,$
samewin=samewin,zx=zx,zy=zy,$
scan=scan,lun=lun,han=han,maxrecs=maxrecs,mytitle=mytitle,$
hlind=hlind,hlval=hlval,hldash=hldash,hlvlines=hlvlines,$
useind=useind,ln=ln,chn=chn)
ARGS:
b[nrecs]: {corget} correlator data to make image of
RETURNS:
img[nchns,nrecs]: float image of last sbc displayed (before clipping).
KEYWORDS:
clip[]: float value to clip the normalized data to. Default is
.02 (of Tsys). If clip has 1 value then
normalize to (img > (-clip)) < clip. If two value are
provided, then they will the [min,max].
sbc: int 1-4. sub correlator to plot. default is all 4.
brdlist: long A single number whose decimal digits specify the boards
to display. eg brds 1,2,3,4 would be: brdlist=1234
The boards are numbered 1 thru 8.
pol: int 1,2, (3,4 if stokes) polarization to plot. default is 1:polA
win: int window number to plot in. default is 1.
wxlen: int xlen of window. default 700
wylen: int ylen of window. default 870
wxpos: int xpos of lower left edge of window.default 445
wypos: int ypos of lower left edge of window.default:35
samewin: if set then use the current dimension for the image win.
If you are calling corimgdisp in a loop then setting this
(after the first call) lets the user dynamically adjust the
window size,position.
zx: int ..-3,-2,2,3,4 zoom factor x dimension. Negative numbers
shrink the image positive numbers expand. Negative number
must divide evenly into the number of channels.
zy: int ..-3,-2,2,3,4 zoom factor y dimension (same format as zx)
col[2]: int .. x columns to use to flatten the image in the time
direction. count 0..numberoflags-1. If multiple sbc
plotted then the same cols are used for all sbc. The
default is no flattening in the time direction.
chn: if set then plot vs channel number rather than freq
bpc:{corget} if supplied then this data will be used to do the
bandpass correction. The default is to average over
all of the nrecs.
nobpc: if set then no bandpass correction is done.
median: if set and bpc not provided, then bandpass correct using
the median of the nrecs rather than the average.
ravg: long bandpass correct with a running average of ravg spectra
scan : long if provided,then routine will input scans data into
b[] before making image. In this case you must also
supply the lun keyword to tell the routine where to
read from.
lun : int if scan keyword provided, then you must also supply
this keyword. It should contain the lun for the
corfile that you have previously opened.
sl[]:{scanlist} This array can be used for direct access when the
scan keyword is used. The sl[] (scanlist) array is
returned from the sl=getsl(lun) routine. The routine
scans the entire file recording where the scans start.
han: if set and scan keyword set, then hanning smooth the
data on input.
maxrecs: int if lun used then the max records of a scan to input.
default:300
mytitle:string user supplied tittle instead of scan,srcname,az,za
az,za at top of the plot.
hlind[]: ind index into img array (2nd dimension) to draw
horizontal lines.
hlval : float value to use for horizontal lines (in img units)
default is clip value.
hldash : int The dash lengths to used for the horizontal lines.
2*ldash must divide into x dimension.default is 4
hlvlines:int Number of engths to used for the horizontal lines.
default=1
useind[2]:int if provided then use these indices from data array
0 .. lengthsbc -1
default=1
ln :int linenumber for title..0..33 def:3
extra_=e this allows you to input keywords that will be
passed to the plotting routine. eg title=title..
DESCRIPTION:
Corimgdisp will display a set of correlator records (usually a scans
worth) as an image. By default it will make a separate image for each
subcorrelator (board). The sbc keyword lets you choose just 1 sbc to image.
The data for the image is taken from polA by default. Use the pol keyword
to make the image from the other polarization.
You can input the data outside of this routine (eg with corinpscan) or
use the scan,lun keywords to have corimgdisp input the data directly.
By default, the image is bandpass normalized by the average of all the
records (sbc/avg(sbc) - 1). If the median keyword is used then
avg(sbc) is replaced by median(sbc). The bpc keyword allows you to input
a separate correlator record to use as the normalization.
The col keyword lets you also flatten the image in the time (record)
dimension by specifying the first/last columns to average and then divide
into all of the other columns (the columns are counted from 0). By default
this is not done.
After bandpass correction and flattening in the record dimension, the
image is clipped (by default) to +/- .02 (tsys) and then scaled
from 0 to 256 for the image display. The clip keyword lets you change the
clipping value.
The zx,zy keywords let you scale the image in the x and y dimensions by
integral amounts. Negative numbers will shrink it by that amount (Note:
the negative numbers must divide evenly into the number of channels in
each sbc). -1,0,1 do no scaling. This scaling is only applied for the single
sbc displays. The multiple sbc displays are always scaled to fit
inside the window (700 by 870 pixels).
The scan keyword can be used to input the data directly from a file.
In this case the scannumber is the scan keyword value and the lun keyword
should be set the the logical unit number of the already open file. The
sl keyword can be used to do direct access to the file. It's argument
comes from a call to sl=getsl(lun) prior to calling this routine.
After displaying the image, use xloadct to manipulate the color table.
The routine returns the last image displayed (before clipping).
EXAMPLES:
input a scans worth of data.
print,corinpscan(lun,b,/han) .. input scan with hanning smoothing
1. display the image of all 4 sbc.
img=corimgdisp(b)
2. display only sbc 2, scale y by 2, and x by -2
img=corimgdisp(b,sbc=1,zx=-2,zy=2)
3. display all 4, clip to .015 Tsys , display polB and median filter the
bandpass correction:
img=corimgdisp(b,pol=2,/median)
4. assume scan 12730084 in an on position and 12730085 is the off position.
Suppose you want to display the on image normalized by the off.
scan=12730084L
print,corinpscan(lun,bon,scan=scan,/han)
print,corinpscan(lun,boff,scan=scan+1,/han)
bpc=coravgint(boff) ; compute the average of the off.
img=corimgdisp(bon,bpc=bpc)
5. Have the routine input and plot a scans worth of data:
openr,lun,'/share/olcor/corfile.27sep01.x101.1',/get_lun
scan=12730084L
img=corimgscan(b,lun=lun,scan=scan)
6. Have the routine input and plot a scans worth of data, use the sl
keyword for direct access.
openr,lun,'/share/olcor/corfile.27sep01.x101.1',/get_lun
sl=getsl(lun)
scan=12730084L
img=corimgscan(b,lun=lun,scan=scan,sl=sl)
This routine calls imgflat, and imgdisp for the image scaling and display.
NOTE:
For stokes channels q,u,v the bandpass flattening is done by offseting
the mean bp by 1. and dividing (should actually divide by sqrt(pola*polB)..:
SEE ALSO:
Imgdisp,imgflat
(See /pkg/rsi/local/libao/phil/Cor2/corimgdisp.pro)
NAME:
corimgonoff - make image of an on off pair
SYNTAX ball=corimgonoff(lun,scan,ver=ver,clip=clip,col=col,red=red,han=han,$
bon=bon,boff=boff,sl=sl,cont=cont,_extra=e,$
brdlist=brdlist,pol=pol)
ARGS:
lun :int logical unit number for file that is already open.
scan:long scannumber for the on position scan
KEYWORDS:
ver[2]: float vertical scale (min,max) for the line plot superimposed
on the image: (on/off -1)-median(on/off-1). The default
value is +/- .005 Tsys.
clip[2]: float min/max clipping value for the image (on/off-1).;
The default value is .02 Tsys
col[2]: int columns of image to use for flattening in the time direction
0 through nchannels-1.
han : if set,then hanning smooth the data on input
bon[] :{corget} pass in on scan rather than input it
boff[]:{corget} pass in off scan rather than input it
sl[] :{sl} scanlist array returned from sl=getsl(lun) that can be used
for direct access i/o.
red : if set then the on/off-1 graph will be red rather than
white
cont : don't remove the continuum info from on/off-1 plot
brdlist : sbc to use. 1,2,12,1234..default is all
pol : to use.1,2 12, default is all
_extra : e pass args to corimgdisp
RETURNS:
ball[]:{corget} all of the on,off records.
DESCRIPTION:
corimgonoff will read in an on/off position switch pair of scans from
the file pointed to by lun. scan should be the scan number for the on
position. It will then call corimgdisp passing in the on and off records
so an image of frequency versus records will be made for each board.
The off records will be used for the bandpass normalization for the entire
image. Time for each image goes from bottom to top so the on records start
at the bottom of each image.
The routine computes y=(avg(on)/avg(off)-1.) where avg() averages over
the individual records. If multiple polarizations are present, they will
also be averaged. The routine then plots:
y= y - median(y)
on top of the corresponding image.
The options sl,ver,col,han,clip are passed to corimgdisp. If the
col=[colmin,colmax] is used to flatten the image in time, then dashed
lines will be superimposed on the plot showing the range that was used.
The line plot will normally be done in white (the largest lookup table
value). If the /red keyword is set then the line plot will appear as red.
This looks nicer but it will probably be changed to black if you play with
the color lookup table (via xloadct) after making the image.
The routine pfcorimnonoff will loop through a file calling this
routine for every onoff pair in the file.
EXAMPLES:
1. default processing:
suppose on scan is 127300099L and lun has already been opened.
hanning smooth the data. Call getsl() once to get the scanlist for
random access.
sl=getsl(lun)
scan=127300099L
ball=corimgonoff(lun,scan,sl=sl,/han)
2. blowup the vertical scale of the line plot and make it red.
ball=corimgonoff(lun,scan,sl=sl,/red,ver=[-.002,.002],/han)
3. Add to 2., scaling of the image in the time dimension using columns
625 through 675 .
ball=corimgonoff(lun,scan,sl=sl,/red,ver=[-.002,.002],/han,$
col=[625,675])
SEE ALSO:
corimgdisp,imgdisp,imgflat,pfcorimgonoff
(See /pkg/rsi/local/libao/phil/Cor2/corimgonoff.pro)
NAME:
corinpscan - input a scans worth of data
SYNTAX: istat=corinpscan(lun,b,brecs,sum=sum,scan=scan,maxrecs=maxrecs,sl=sl,
han=han,df=df)
ARGS:
lun: int .. assigned to open file
b[ngrps]: {} corget structs.
brecs[ngrps]: {} corget structs. optional argument. If
/sum is set and brecs is supplied, then the summary
will be returned in b and the individual recs brecs.
..
KEYWORDS:
sum: int .. if not zero then return avgerage of the scan records
in b[]. The header is from the first record
(with modifications). If brecs is also supplied then
the individual recs will be returned in brecs;
maxrecs: int .. max number of records we can read in when not
computing summary. default is 500. If you have
more records in scan then set this arg.
scan: long position to scan before inputing.
sl[]: {sl} returned from getsl(). If provided then direct access
to scan is available.
han: if set then hanning smooth the data.
df : if set then remove the digital filter bandpass from each record.
istat: int .. returns 1 ok, 0 error.
DESCRIPTION:
corinpscan will input an entire scan.
It will return:
- all the records in b[] ( /sum not set )
- only the average of the scan in b( /sum set, brecs not specified)
- the average of the scan in b, all recs in brecs[](/sum set, brecs specified)
The user can position to the scan using the keyword scan. Input must
start at the first record of the scan (if no positioning is selected).
The returned data will be a single {corget} structure (if /sum set) or
and array b[ngrps] (or brecs[ngrsp] of {corget} structures
(one for each integration in the scan.
When /sum is set then the header returned will be from the first
group of the scan with the following modifications:
h.std.gr,az,chttd will be the average value.
h.std.grpnum will be the number of groups input.
h.cor.lag0pwrratio will be the average
mod history:
31jun00 updated to new corget format
04jul00 return average az,gr,ch position and avg lag0pwrratio
13aug00 added brecs optional arg
26dec02 fixed bug where last scan of file skipped if 1 record long and
we read sequentially with corinpscan().
(See /pkg/rsi/local/libao/phil/Cor2/corinpscan.pro)
NAME:
corinterprfi - interpolate spectra across rfi.
SYNTAX: binterp=corinterpolrfi(b,fsin=fsin,deg=deg,ndel=ndel,verb=verb,$
_extra=e)
ARGS:
b: {corget} Corget structure to interpolate.
KEYWORDS:
fsin: int degree of sin fit to use (default 8).see corblauto
deg : int degree of polynomial fit to use (default 1).see corblauto
ndel: int The number of adjacent points to exclude when a bad point
is found (default 5). see corblauto
verb: int if -1 then plot the baseline fitting while it's being done.
_extra: Any extra keywords will be passed to corblauto (see corblauto
for a list of keywords you can use).
RETURNS:
binterp:{corget} the {corget} data strucuture b with interpolation
over the excluded (outlying) data points.
DESCRIPTION:
This routine finds bandpass outliers (by fitting polynomials and sines
to the bandpass using corblauto). It then interpolates across the
outliers using the idl interpol() routine (with the default linear
interpolation).
Normally the "outliers" include the edges of the bandpass
(since they are hard to fit). This is not desired when fitting for
"just rfi". To correct this, the routine searches in from both edges
of the bandpass until it finds a "non outlier" point. All the outliers
from the edge to this point are not interpolated across.
EXAMPLE:
; input a scan, create a median bandpass, and then interpolate
; across any rfi. Use this average bandpass as the bandpass correction
; for the image display routine corimgdisp().
istat=corinpscan(b)
bavg=cormedian(bavg)
bpc=corinterprfi(b,verb=-1)
img=corimgdisp(b,bpc=bpc)
(See /pkg/rsi/local/libao/phil/Cor2/corinterprfi.pro)
NAME:
corlist - list contents of correlator data file
SYNTAX: corlist,lun,recperpage,scan=scan,sl=sl,ast=ast
ARGS:
lun: int assigned to open file
recperpage: lines per page before wait for response. def:30
KEYWORDS:
scan: long position to scan before listing. def:rewind then list
sl[]: {sl} returned from getsl(). If provided then direct
access is available.
ast: if set then print local time rather than lst for
start of each scan
(See /pkg/rsi/local/libao/phil/Cor2/corlist.pro)
NAME:
cormask - interactively create a mask for each correlator board
SYNTAX: cormask,b,cmask,edgefract=edgefract
ARGS:
b: {corget} correlator data
KEYWORDS:
edgefract[n]: float . if provided then create a mask ignoring
lagsPerSbc*edgefract lags on each edge of the bandpass.
In this case the routine returns immediately with no
user interaction.
If edgeFraction is a single number, then use it
on both sides of the bandpass. If n=2 then
use edgefraction[0] for the left side and
use edgefraction[1] for the right side.
RETURNS :
cmask : {cormask} structure holding a mask for each board
DESCRIPTION:
cormask lets the user interactively create a mask for each board of
a correlator dataset. The routine will pass back a cormask structure that
contains 1 mask for each board. If the user does not specify a mask for a
sbc, the last mask entered will be duplicated. If no masks are ever defined,
then masks with all 1's will be returned. The program menu is:
KEY ARGS FUNCTION
m .. define mask for current board
h h1 h2 .. change horizontal scale for plot to h1,h2
v v1 v2 .. change vertical scale for plot to v1,v2
b boardnum.. board number to process (1..nbrds)
q .. quit
By default the first board is displayed on entry. The user should adjust the
horizontal and vertical scale with h h1 h2 and v v1 v2. Each time one of
these is entered the plot will be redisplayed with the new limits
(do not use commas here..).
When the plot shows the correct limits, enter m to define the mask. The
leftmost mouse button is used to define the starting and ending portions
of the data that will be used for the mask. You can have as many of these
sections as you want. When you are done with the mask, click the right
mouse button. The routine will then autoincrement to the next board. You
can change the order of the next board or redo a board with the b key.
q is how you exit the routine.
Suppose you have 2 boards with 1024 and 2048 lags. Then the cormask structure
returned will be:
int cmask.b1[1024]
int cmask.b2[2048]
It is an anonymous structure whose elements will change depending on the
{corget} structure passed in.
(See /pkg/rsi/local/libao/phil/Cor2/cormask.pro)
NAME:
cormaskmk - create the cormask structure from the corget structure
SYNTAX: cmask=cormaskmk(b,edgefract=edgefract,ones=ones,raw=raw)
ARGS:
b[n]: {corget} correlator data
KEYWORDS:
edgefract[1or2]: float . if provided then create a mask ignoring
lagsPerSbc*edgefract lags on each edge of the bandpass.
if edgefract has two entries then use
edgefract[0] for the left and edgefract[1] for the right.
ones: if set and edgefraction not supplied, return the
mask filled with ones rather than zero.
raw: if the raw keyword is supplied, return a mask with
1 brd, 1 pol of raw pnts.
RETURNS :
cmask : {cormask} structure holding a zeroed mask for each board
if edgefract was supplied then the values inside
the edgefraction will be set to 1.
DESCRIPTION:
create a cormask structure from the corget data structure.
The mask structure contains float arrays:
cmask.b1[lensbc1,npol1]
cmask.b2[lensbc2,npol2]
cmask.b3[lensbc3,npol3]
cmask.b4[lensbc4,npol4]
The arrays will be zeroed. If the edgefract keyword is supplied, then
edgefract*lensbc on each edge will be set to zero and the reset will
be set to 1.
This routine is called by cormask and corcumfilter.
(See /pkg/rsi/local/libao/phil/Cor2/cormaskmk.pro)
NAME:
cormath - perform math on correlator data
SYNTAX: bout=cormath(b1,b2,sub=sub,add=add,div=div,mul=mul,$
sadd=sadd,ssub=ssub,sdiv=sdiv,smul=smul,$
norm=norm,mask=mask)
ARGS:
b1[n] : {corget} first math argument
b2[n or 1]: {corget} second math argument
KEYWORDS:
sub: if set then bout= b1-b2
add: if set then bout= b1+b2
div: if set then bout= b1/b2
mul: if set then bout= b1*b2
sadd: if set then bout= b1+sadd (scalar add)
ssub: if set then bout= b1-ssub (scalar subtract)
sdiv: if set then bout= b1/sdiv (scalar divide)
smul: if set then bout= b1*smul (scalar multiply)
norm: if set then normalize(b1) one argument, or norm only
normalize(b2) if two args and a math operation
If a mask defined, normalize within the mask.
If a math operation is also defined, perform the
math operation after normalization.
mask :{cormask} a mask to use with norm keyword (see cormask).
Note that cormath only uses the polA mask
of each board. If a 2nd mask is defined for polb
(eg corcumfilter), the pola mask is used for both pols
RETURNS:
bout[n]: {corget} result of arithmetic operation.
DESCRIPTION:
Perform arithmetic on correlator data structures. The operation performed
depends on the keyword supplied.
The keywords sub,add,div, or mul perform the operations on the two
structures b1,b2. There are two modes of combining b1,b2:
1.b1 and b2 are the same length:
bout[i]=b1[i] op b2[i]
2.b1[n] and b2[1] then
bout[i]= b1[i] op b2[0] ..
The 2nd form could be used when bandpass correcting a set of records
with a single bandpass.
The keywords sadd,ssub,sdiv,smul take a single scalar value and
perform the operation on all elements of b1[n] (b2 is ignored if provided).
The scalar value provided can be:
1. a single number. In this case all sbc,pol, freqbins are operated on
with this number.
eg: bout[i]=b[i] op scalarvavlue
2. an array of numbers with Numpol*Numsbc entries. In this case each
pol,sbc is operated on with its number.
eg. if b has 4sbc by 2 pols, sclarr[2,4] would have a separate
number for each sbc and the multilply would be:
b[n].(sbc).d[*,pol] op sclarr[pol,sbc]
The keyword norm will normalize b1 (so mean value is 1) if only b1 is
provided. If 2 args are provided then normalize b2 before performing the
operation. If a mask is provided, then compute the normalization over the
mask.
EXAMPLE:
print,corget(lun,b1)
print,corget(lun,b2)
badd=cormath(b1,b2,/add)
;
; subtract a smoothed bandpass from a single rec
;
corsmo,b1,bsmo,smo=15
bout=cormath(b1,bsmo,/sub)
;
; bandpass correct a smoothed bandpass from a single rec
;
corsmo,b1,bsmo,smo=15
bout=cormath(b1,bsmo,/div,/norm)
; multiply data by a scalar value 1.34
bout=cormath(b1,smul=1.34);
; input an entire scan, compute normalized band pass and then bandpass correct
; each record
print,corinpscan(lun,binp)
bbandpass=coravg(binp,/norm)
bout=cormath(binp,bbandpass,/div)
;
; input on,off data and normalize on to off record by record
; then average the spectra. It is probably simpler to add the records
; before dividing but you could do it this way to verify that the results
; are the same.
;
print,corinpscan(lun,bon1) ; bon1 is n recs
print,corinpscan(lun,boff1) ; boff1 is n recs
bonoff1=cormath(bon1,boff1,/div,/norm);bonoff1 is n recs
bonoff1=coravg(bonoff1) ; bonoff1 is 1 averaged rec
; The faster way would be
print,corinpscan(lun,bon2,/avg) ; bon2 is 1 averaged record
print,corinpscan(lun,boff2,/avg) ; boff2 is 1 averaged record
bonoff2=cormath(bon2,boff2,/div,/norm)
; Use a mask (see cormask() to normalize the bandpass where the mask is
; non-zero instead of the entire bandpass.
bonoff2=cormath(bon2,boff2,/div,/norm,mask=cmask)
SEE ALSO: cormask,coravg,corinpscan
(See /pkg/rsi/local/libao/phil/Cor2/cormath.pro)
NAME:
cormedian - median filter a set of integrations
SYNTAX: bfiltered=cormedian(b)
ARGS: b[] : {corget} array of corget structures
RETURNS:
bfiltered : {corget} after median filtering each sbc,brd of b
DESCRIPTION:
if b[] is an array of n corget structures, cormedian() will compute the
median (over the records). It will return a single averaged corget struct.
This is the same as coravg() except that it uses the median rather than
the mean.
(See /pkg/rsi/local/libao/phil/Cor2/cormedian.pro)
NAME:
cormon - monitor data from file.
SYNTAX: cormon,lun,b,m=pltmsk,han=han,vel=vel,pol=pol,scan=scan,
avgscan=avgscan,quiet=quiet,no1rec=no1rec,sl=sl,delay=delay
secperwait=secperwait
ARGS:
lun: int assigned to open file.
KEYWORDS:
m: which sbc to plot.. bitmask b0->b3 for brd1->4
han: if set then hanning smooth the data
pol: limit the polarizations to look at 1- polA, 2-polB
vel: if set then plot versus velocity. def:freq.
scan: long position to scan before starting the monitoring
avgscan: if set then only plot scan averages
quiet : if set and avgscan set, then do not print out rec numbers
no1rec : if set and avgscan set, do not plot scans that contain
only 1 record.
sl[] : {sl} if scan list is provided (see getsl()) then
only display the scans in sl. When done return
rather than wait for next rec.
delay : float seconds to wait between plots (useful if you are
scanning a file offline. default:0
secperwait: float When no data is available, the routine will
wait 1 sec before checking again. You can override
this (if you are dumping faster than once a second.
RETURNS:
b: {corget} data from last read
DESCRIPTION:
Monitor the data in a file. When the routine hits the end of file it
will continue plotting as new data becomes available. This routine
is normally used for online monitoring. To exit the routine use ctrl-c.
Use the avgscan keyword if you want to only plot the scan averages. This
is a good thing to do if you are observing remotely and the network
connection can not keep up with displaying every record.
An annoyance with the /avgscan setting is that it always displays the
scan average of the last completed scan. Suppose you are doing 5 minute
on/offs followed by 10 second cal/offs. You will wait 5 minutes for the
on scan, 5 minutes for the off scan, and then in 10 seconds the cal on
scan will be plotted (losing the off scan plot). The /no1rec keyword tells
cormon to not bother plotting scans that only have 1 record (the calon,
caloff records are separate scans of 1 record each). This will give
you the complete 5 minutes to look at the 5 minute off.
The sl option can be used to scan files offline. The routine will
return when done rather than wait. It also allows you to plot a subset
of a file.
EXAMPLES:
cormon,lun .. monitor every rec, all sbc,pol
cormon,lun,m=1 .. monitor every rec, sbc 1, both pol.
cormon,lun,m=2 .. monitor every rec, sbc 2, both pol.
cormon,lun,m=5 .. monitor every rec, sbc 1,3, both pol.
cormon,lun,m=8,pol=1 ..monitor every rec, sbc 4, pol A.
cormon,lun,/avgscan.. monitor scan averages only
cormon,lun,/avgscan,/no1rec.. monitor scan averages,skip scans with 1 rec
cormon,lun,/avgscan,/quiet.. monitor scan averages,don't output rec #'s
use the sl option to plot all records from receiver 6
sl=getsl(lun)
ind=where(sl.rcvnum eq 6,count)
if count gt 0 then cormon,lun,sl=sl[ind]
SEE ALSO:
corplot, getsl()
(See /pkg/rsi/local/libao/phil/Cor2/cormon.pro)
NAME:
cormonall - monitor correlator data, bandpass,rms,img
SYNTAX: cormonall,lun,vel=vel,pol=pol,brdlist=brdlist,maxrecs=maxrecs,$
spcavg=spcavg,norms=norms,nospectra=nospectra,noimg=noimg,$
base=base,edgefr=edgefr,vrms=vrms,vavgspc=vavgspc,$
imgpol=imgpol,ravg=ravg,col=col,clip=clip,imgmedian=imgmedian,$
imgbrdlist=imbbrdlist,vonoff=vonoff,cont=cont,$
sl=sl,step=step,scan=scan,nolut=nolut,pltmsk=pltmsk,samewin=samewin
ARGS:
lun: int file descriptor to input from. No file positioning done.
KEYWORDS for plotting:
vel: if set then plot x axis is velocity. default topocentric freq.
pol: int/char limit plots to just 1 polarization (if multiple pols per
board). 'a',1 for polA or 'b',2 for polB, 3,4=stokes,u,v
default is all pols.
brdlist: long specify the boards to use. Count 1..4 (or 1..8 with alfa).
To use brds 1,2,3,4 use brdlist=1234 . This replaces
pltmask=. See imgbrdlist for the images
df: if set then remove digital filter bandpass from average spectra
KEYWORDS for screen layout:
norms: if set then don't display rms plots
nospectra: if set then don't display individual spectra
noimg: if set then don't display the image.
maxrecs: int max number of records in a scan (for image or rms). By
default the routine allows up to 300 recs. If your
scans have more than 300 recs then set maxrecs to the
maximum number you will have (but not too much bigger since
it will allocate a buffer to hold an entire scan).
recsimg: int minimum numbers or records for an image to be made. Def:40
recsrms: int minimum numbers or records for rms to be made. Def:10
spcavg: int if non zero then also plot avg spectra. Update average plot
every spcavg recs. 1--> redisplay average every record.
2--> redisplay average every 2 recs etc..
base: int if non zero then remove baseline of order "base" from
average spectra. Ignore edgefr fraction of bandpass on
each side
edgefr[n]:flt fraction of bandpass on each side to ignore in fit. If
n=1 then use edgefr[0] for both edges. If n=2 then use [0]
for the left edge and [1] for the right edge. If edgefr is
not supplied, then use .15 if not /df and .1 if /df set.
samewin: if set then you are recalling cormonall with same parameters.
don't bother to recreate the windows. If you are looping thru
multiple files, this lets you look at the old pictures while
inputting the new ones.
KEYWORDS for image:
imgpol:int/strin polarization of image to display (if multiple pols
per board). a,1:polA, b,2:polB, 3,4 for cross spectra
default is polA
imgbrdlist:int Boards to use for image display. This has the same coding
as brdlist keyword above except that it applies to the image
display. This lets you display all of the board spectra but
make images of only a subset of them. This is handy if the
nrecs*nbrds is too large to be effectively displayed on the
screen.
ravg: int use a running average of (ravg) spectra for the bandpass
correction. The default is to average over the entire scan.
Note that for position switched spectra, the display
of the combined on,off image will always use the mean
of the off for the bandpass correction.
col[2]: int average frequency channels between col[0] and col[1] and then
divide the image by this. flattens the image in the time
direction.
clip[2]: flt By default the image max/min display is scaled to .02
of Tsys. The clip[2] will clip the display to the min,max of
clip[0],clip[1].
imgmedian: If set (/imgmedian) then use the median rather than the mean
when computing the bandpass correction (for flattening the
image). This helps if there is strong continum in the image.
vonoff[2]:flt images of on/off position switching will include a lineplot
of [(on/off-1) - median(on/off-1)]. The default vertical scale
for this plot is +/- .005 Tsys (25Mhz, 1024 channels, and
300 secs on gives a noise sigma of .0005). veronoff[2] lets
you change the min,max to veronoff[0],veronoff[1].
vrms[2]: flt min, max value for rms plot. default is autoscale
cont: if set then don't remove the continuum from the on/off line
plot that is drawn on top of the image. Normally the plot
is (on/off-1) - median(on/off-1). Setting this keyword plots
(on/off-1). You will probably have to also use veronoff[] to
keep the drawing within the yaxis range.
vavgspc[2]:flt min, max value for plot of average spectra. If you have
baslined it then it should be about zero.
nolut: by default the routine will load a linear grey scale ramp
into the color lookup table (lut) for the image. If you have
already loaded the lut (eg xloadct) then /nolut will tell
this routine to not modify the lut.
KEYWORDS misc:
sl[]: {scanlist} The scanlist returned by getsl() can be used to
speed the processing of the file.
step: if set then prompt the user to continue after each scan.
han: if set then hanning smooth the data on input
scan: long position to scan before monitoring
pltmsk: int replaced by brdlist=
bitmask to display a subset of the boards.
1 brd1, 2 brd2, 3 brd1&brd2, 4 brd3, 5 brd1&brd3,6 brd2&brd3,
7 brd1&brd2&brd3, 8 brd4,9 brd1&brd4,10 brd4&brd2...
DESCRIPTION:
cormonall processes a file displaying the individual spectral records,
the rms by channel of each scan, an image of frequency by time for
each scan, and a running average and final average of the spectra.
The data is read from the file descriptor lun. It must be opened
prior to calling this routine. The program starts reading from the current
position in the file. Use rew,lun or print,posscan(lun,scan) if you want
to position the file before processing.
By default all spectral records are displayed, all scans with more than
10 records will have the rms by channel computed, and all scans with 40
records or more will have an image made. The norms, nospectra, noimg keywords
can be used to not display one or more of these options.
To include the running and final average of the spectra use the keyword
spcavg. It will display the running average every spcavg records. Two
windows will be generated: one for the running average of a scan and a
second for the complete scan average. Scans with 1 record are not averaged.
These two window occupy the same screen area as the image display.
If images are to be plotted the color lookup table (lut) will be
loaded with a greyscale linear ramp. If you have already loaded the
lut (eg with xloadct) then /nolut will tell the program to not touch
the current values.
The routine will processes all of the records in a file and then wait
for the file to grow so the next record can be displayed. You can hit
any key and the program will prompt you to enter return to continue or
q to quit. This allows you to stare at a plot or image for awhile before
it gets overwritten.
You can limit the number of boards displayed with the brdlist= keyword.
To display only boards 1,2, and 3 use brdlist=123
SPECTRA/RMS/AVG
The spectra, rms, and average windows are plotted versus topocentric
velocity. The vel keyword can switch this to velocity (this keyword does not
affect the image x axis). The vrms and vavgspc can be used to change the
vertical scale of the rms plot and the average spectra plot.
IMAGES
An image of each scan with more than 40 records is made (unless /noimg
is selected). If two polarizations per board are taken then the image will
default to pol A. The imgpol keyword lets you choose the other polarization.
The image will be bandpas corrected by the mean spectra of the entire scan
(or median if /median is supplied). The display range is scaled to
+/- .02 Tsys. You can change this range with the clip keyword. You can
flatten the image in the time direction with the col (columns) keyword.
It will average over these columns and then divide this into the image.
If on/off position switch data is found then the images work a little
differently. The first image to appear is the on. The second image has the
on records in the bottom half of the image and the off records in the
top half of the image. The average off scan is used for the bandpass
correction for both the entire image (on and off scans)
A line plot is overplotted in the center of the image. It is:
(on/off-1) - (median(on/off-1))
The default vertical scale for this plot is +/-.005 Tsys. You can change
the range with the vonoff keyword. By default any continuum in the on or
off is removed (via the -median(on/off-1). The /cont keyword will leave the
continuum in the line plot. You will probably have to use the vonoff[]
keyword to readust the vertical scale.
EXAMPLES:
; load the grey scale lookup table
loadct,0
; open a file and display it..
openr,lun,'/share/olcor/corfile.29may02.p1525.1',/get_lun
; default call
cormonall,lun
; at the end of the file, hit a key to get the prompt back, then enter
; q to exit the routine
; change the dynamic range of the image to +/-.01 Tsys. blowup any
; on/off-1 plots to .004 Tsys
; only plot sbc 1 and 3
rew,lun
cormonall,lun,clip=[-.01,.01],vonoff[-.004,.004],brdlist=13
SEE ALSO:
The routine calls corplot,corrms,corimgdisp,and corimgonoff
(See /pkg/rsi/local/libao/phil/Cor2/cormonall.pro)
NAME:
cormovie1d - make a movie of 1d spectra
SYNTAX: cormovie1d,b,loop=loop,win=win,xdim=xdim,ydim=ydim,$
delay=delay,yauto=yauto,$
corplotargs..
ARGS:
b[]: {corget} array of corget data to plot
KEYWORDS:
loop: If set then continue looping through the data array.
you can stop by hitting a key and then q for quit.
By default the routine exits after making 1 pass through
the array.
win : int The window to use for the plotting. By default it uses
the first free window above 0.
xdim: int if supplied, then the x dimension for the plot window
in pixel units. the default is 640
ydim: int if supplied, then the y dimension for the plot window
in pixel units. the default if 510
delay: float The number of seconds (or fractions of seconds) to wait
between each plot. The default is 0.
yauto: If set then automatically set the vertical scale by
finding the min,max of all the data.
corplotargs: You can pass in any keywords that are accepted by
corplot (eg m=2,/vel, etc..)
DESCRIPTION:
Make a movie of an array of corget spectra. It will make a call
to corplot for each spectra in the array. Pixwins are used to
eliminate flashing between plots. By default the routine exits after
make 1 pass through the array. The /loop keyword will cause
the routine to continue looping through the array until the user
asks to quit.
While plotting, the user can resize the window with the mouse. The
routine will change to the new window size.
If the user hits any key while plotting, the routine will pause
and print :
idl> return to continue,s step,c exit step mode, q to quit
The user can enter:
return - the routine will continue whatever it was doing.
s - the routine will go to single step mode. The user will
have to hit return after each plot.
c - go back to continuous plotting. This exits step mode.
q - quit. The routine will exit.
When using the routine, you should set the vertical scale before calling
the routine or use the /yauto keyword. When the vertical scale is not
set then each spectra displayed will autoscale to the min,max.
Any corplot keywords will be passed through to the corplot call.
EXAMPLES:
1. input a scan
IDL> istat=corinpscan(lun,b)
2. make a movie of all sbc in b, use auto scaling.
IDL> cormovie1d,b,/yauto
3. set the vertical scale and loop forever. Only plot sbc 1.
IDL> ver,0,2
cormovie1d,b,/loop,m=1
SEE ALSO: corplot
(See /pkg/rsi/local/libao/phil/Cor2/cormovie1d.pro)
NAME:
cornext - input and plot next rec from disc
SYNTAX: cornext,lun,b,m=pltmsk
ARGS:
lun: int assigned to open file.
b: {corget} data from last read
KEYWORDS:
m: which sbc to plot.. bitmask b0->b3 for brd1->4
(See /pkg/rsi/local/libao/phil/Cor2/cornext.pro)
NAME: coronl - open the online datafile SYNTAX: - lun=coronl(calfile=calfile) ARGS: calfile if set then use the calibrate file rather than corfile DESCRIPTION: open the online datafile /share/olcor/corfile. Return the lun to access the file. This should only be used by the people that are actively taking data on the telescope.
(See /pkg/rsi/local/libao/phil/Cor2/coronl.pro)
NAME:
coroutput - output a cor data structure to disc
SYNTAX: coroutput,lunout,b
ARGS:
lunout : int lun for file to write to.
b[]: {corget} a corget structure or an array of corget structs.
RETURNS:
nothing
DESCRIPTION:
Output a corget structure to disc. You can then read this back in
with corget at a later time.
WARNINGS:
1. The {corget} structure has an extra element after the hdr structure. It
is stripped off before writing it out so that data can be read in with
corget().
2. On input corget() will flip the input data order if it is in
decreasing frequency order on disc (determined by a bit in h.cor.state).
This output routine assumes that the data has already been converted to
increasing frequency order and it will set the bit in h.cor.state
so that it will not be reflipped when input with corget().
3. corget() will scale the data on input using the total power (this is
the nine level total power correction(). It is ok to apply the
scaling multiple times since it is:
scale*(data/Mean).. but if you are outputing data that is the
rms/mean rather than spectral data, then be sure you use the /noscale
option to corget().
(See /pkg/rsi/local/libao/phil/Cor2/coroutput.pro)
NAME:
corplot - plot the correlator data.
SYNTAX : corplot,b,brdlist=brdlist,vel=vel,rest=rest,off=ploff,pol=pol,
over=over,nolab=nolab,extitle=extitle,newtitle=newtitle,
col=col,sym=sym,chn=chn,cmask=cmask,$
xtitle=xtitle,ytitle=ytitle,m=pltmsk,flag=flag,$
lnsflag=lnsflag,font=font,cs=cs,ncs=ncs
ARGS:
b: {corget} data to plot
KEYWORDS:
brdlist: int number containing the boards to plot. To plot boards
1,3,5,7 use brdlist=1357.By default all boards are plotted.
This replaces the m= keyword.
vel: if set, then plot versus velocity
rest: if set, then plot versus rest frequency
off: float if plotting multiple integrations then this is the
offset to add between each plot of a single sbc. def 0.
pol: int if set, then plot only one of the pol(1,2 ..3,4 for stokes)
over: if set then overplot with whatever was previously plotted.
This only works if you plot 1 sbc at a time (eg m=1).
nolab: if set then don't plot the power level labels in the middle
of the plot for stokes data.
extitle: string add to title line
newtitle: string replace title line
xtitle: string title for x axis
ytitle: string title for y axis
col[4]: int Change the order of the default colors. The
order is PolA,polB,stokesU,stokesV colors.
values are:1=white(on black),2-red,3=green,4=blue,5=yellow
sym: int symbol to plot at each point. The symbols are labeled
1 to 9. Positive numbers only plot the symbol (no
connecting lines). Negative numbers (-1 to -9) plots the
symbol and a connecting line. sym=10 will plot in
histogram mode.
chn: if set then plot vs channel
cmask:{cmask} if provided then overplot the mask on each plot
m: int bitmask 0-15 to say which boards to plot. This has been
deprecated (see brdlist for the replacement).
flag[]: float if present then flag each of these channels with a vertical
line. The units should be that of the x axis (Mhz,chan, or
km/sec).
lnsflag: int linestyle for flag. 0=solid 1-dashed,2=dashed,3=dots
font : int font to use. def: reg. 1 =truetypelinestyle for flag. 0=solid 1-dashed,2=dashed,3=dots
cs : float charsize. def set by number of rows
ncs : float charsize for notes (tit).
DESCRIPTION:
corplot will plot the correlator data in b. It can be from a single
record (corget) or it can be an array from corinpscan. You can plot:
- any combination of sbc using m=pltmsk and pol=pol
- by topocentric frequency (defaul), rest frequency, or by velocity
- make a strip plot if b is an array by using off=ploff.
If the data passed in has been accumulated by coraccum() then the
average will be computed before plotting.
EXAMPLES:
corget,lun,b
corplot,b plot all sbc
corplot,b,brd=1 plot first sbc only
corplot,b,brd=2 plot 2nd sbc
corplot,b,brd=3,pol=2 plot 3rd sbc,polB only
corplot,b,brd=4,pol=1 plot 4th sbc, pola only
corplot,b,brd=13,/vel plot 1st and 3rd sbc by velocity
istat=corinpscan(lun,b) .. input entire scan
corplot,b,brd=1,pol=1 plot sbc 1, pola, overplotting
corplot,b,brd=2,off=.01 plot sbc 2, with .01 between records
.. overplot sbc 1 from b[0] (red) and b[1] yellow
corplot,b[0],brd=1
corplot,b[1],brd=1,/over,col=[5,5]
NOTE:
use hor,ver to set the plotting scale.
oveplotting only works if you display 1 sbc at a time.
SEE ALSO:
corget,coraccum. hor,ver in the generic idl routines.
(See /pkg/rsi/local/libao/phil/Cor2/corplot.pro)
NAME:
corplotrl - plot correlator data flagging recomb lines.
SYNTAX : corplotrl,b,recombI,m=pltmsk,pol=pol,extitle=extitle,rest=rest,$
across=across,down=down,font=font,cs=cs,ncs=ncs
ARGS :
b: {corget} corget structure to plot
recombI[]: {} recombination lines to flag (see recombsearch()).
KEYWORDS :
m: int bitmask 0-15 to say which boards to plot
1=brd 1 , 2=brd2 , 4=brd3, 8=brd4,
3=brd12 , 5=brd13 , 7=brd123, 9=brd14..
pol : int The polarization to plot (1 or 2). By default
all polarizations are plotted.
extitle:string add to title line
across : int if supplied then the number of plots across the page
down : int if supplied then the number of plots down the page
font : int font to use (1=truetype)
cs : float chars=cs value to use for plot
ncs : float chars=ncs value to use for note labels &$
DESCRPIPTION:
Plot correlator structure flagging the recombination line freqeuncies.
The input data b should be a single corget structure. The recombination
line information (recombI) is from the routine recombsearch().
The data is plotted in topocentric frequencies. The rest frequencies
of each recombination line will be doppler shifted by the doppler
shift used for datataking (the doppler shift for the rf band center or the
doppler shift for the center of each sub correlator). The routine
uses color to differentiate the atoms. The linestyle (dots, dashes, etc..)
is separates the transitions (alpha,beta,..). Up to 8 atoms and 5
transition steps can be displayed at once.
You can limit the subbcorrelators that are plotted with the m=pltmsk
keyword. m is a bitmask where a 1 in a bit means plot the sbc. bit0(1)=sbc1,
bit1(2)=sbc2, bit2(4)=sbc3,bit3(8)=sbc4. Any combination of these bits
can be used.
EXAMPLES:
; Flag any recombination lines between 1664 and 1680 for transitions 1
;thru 5.
; .. input correlator data into b
; .. look for recomb lines in this range
n=recombsearch(1664,1680,recombI,lineStepR=[1,5])
corplotrl,b,recombI
; search for transitions 1 thru 10, plot the first 5 then the second five
n=recombsearch(1664,1680,recombI,lineStepR=[1,10])
ind1=where(recombi.lineStep le 5,count)
ind2=where(recombi.lineStep gt 5,count)
; steps 1.. 5
corplotrl,b,recombI[ind1]
; steps 6..10
corplotrl,b,recombI[ind2]
SEE ALSO: recombsearch(), recombfreq() in idl/gen
NOTE:
The recombination line flagged will be offset from the measured value
if the doppler shift used was different than that of the emission region.
(See /pkg/rsi/local/libao/phil/Cor2/corplotrl.pro)
NAME:
corplttpza - plot total power vs za for each scan
SYNTAX: corplttpza,pwr,horza=horza,ver1=ver1,ver2=ver2,printscan=printscan,
title=tile,sbc=sbc,scanpos=scanpos,avgpol=avgpol
ARGS :
pwr[] :{corpwr} power data returned from corpwr
KEYWORDS:
horza[2]: float za range to plot.def.. that of data
ver1[2]: vertical range top plot. default:range of data
ver2[2]: vertical range bottom plot. default:range of data
printscan: if set then print scannumber on bottom plot
title : string title to use for top plot.
sbc : int 0 or 1 to plot 1st or 2nd sbc of board.
For setups with 2 polarizations per board this
lets you choose polA or polb. The default is 0=polA
scanpos: float yposition for the start of the scan number list (if
printscanis requested). Units are fration of the
window (measured from the bottom). The default value
is .15
avgpol: if set then average polarizations
DESCRIPTION:
Display the total power versus zenith angle for all scans in a file.
You must first input the total power information from the file using:
nrecs=corpwrfile(fileanme,p,maxrecs=20000)
Corplottpza will then make two plots:
1. the total power versus zenith angle.
2. totalpower/median(totalpowerInscan) -1
By default it will plot the total power from the first sbc of a board. The
sbc keyword can be used to display the total power from the second
sbc (if there were two sbc/board). Scans with only 1 record will be
ignored. Each scan will be plotted in a different color (cycling through
10 colors). If the colors do not show up try ldcolph to load the color
lookup table.
This routine is handy to look at onoff position switching data since
the on and off position will follow the same za track. The color order
printed on the left helps you determine which was the on and which was
the off (the on is taken before the off so its color will appear first
in the list).
The bottom plot shows the fractional change in power with za. Ideally this
change should be the same for the on,off positions (if there is no
continuum present). Rapid changes may be interference or the gain
changing with temperature.
The keyword horza lets you adjust the horizontal scale.
The vertical scale can be adjusted with ver1[] (top plot) and ver2[]
bottom plot.
The keyword printscan will print the scannumbers (with the same color coding
as the plots) at the bottom of the page.
EXAMPLE:
nrecs=corpwrfile(fileanme,p,maxrecs=20000);
ldcolph .. make sure colortable is loaded.
corplttpza,p
corplttpaz,p,ver1=[.8,1.5],horza[0,20] ..modify vertical/horizontal scale
corplttpza,p,sbc=1,/printscan ..plot other polarization
print scan numbers at bottom.
(See /pkg/rsi/local/libao/phil/Cor2/corplttpza.pro)
NAME:
corposonoff - process a position switch scan
SYNTAX: istat=corposonoff(lun,b,t,cals,bonrecs,boffrecs,sclcal=sclcal,
sclJy=sclJy,scan=scan,maxrecs=maxrecs,sl=sl,dbg=dbg,han=han,
swappol=swappol,median=median)
ARGS :
lun: int .. assigned to open file.
b: {corget} return (on-off)*X here..
X is:
sclCal true : (kelvins/corcnt)/.. units are then K
sclcal false: off .. units are TsysOff
t: {temps} total power temp info contains..
t.src[2,8] .. on/off-1 (same as lag0pwrratio)
t.on[2,8] .. temp on src
t.off[2,8] .. temp off
units are Kelvins if /sclcal else Tsys Off
cals[] {cal} return cal info. see coronoffcal..
this is an optional argument
bonrecs[] {corget} return individual on records. This is an
optional argument. The units will be K (if /sclcal),
Jy if /sclJy, else correlator units.
boffrecs[] {corget} return individual off records. This is an
optional argument. The units will be K (if /sclcal),
Jy if /sclJy, else correlator units.
KEYWORDS:
sclcal: if not zero then scale to kelvins using the cal on/off
that follows the scan
scljy: if set then convert the on/off-1 spectra to Janskies.
It first scales to kelvins using the cals and then applies
the gain curve.
scan: position to first record of scan before reading.
sl[]: {sl} returned from getsl(). If provided then direct
access is available for scan.
maxrecs: maximum records in a scan. default 300.
han : if set then hanning smooth the data.
swappol: if set then swap polA,polB cal values. This can be
used to correct for the 1320 hipass cable reversal or
the use of a xfer switch in the iflo.
median : if set then median filter rather than average
each scan
_extra : calval=caltouse[2,nbrds] .. for corcalonoff..
RETURNS:
istat int 1 gotit, 0 didnot complete
cal value to use polA,B by board instead of
value in file
DESCRIPTION:
Process a position switch on,off pair. Return a single spectrum
per sbc of (posOn-posOff)/posOff. The units will be:
Kelvins if /sclcal is set
Janskies if /scljy is set .
TsysUnits if none of the above is set.
The t structure holds the total power info for (on/off -1), on source,
and off source. These units are either Kelvins (/usecal or /usejy set) or
TsysOff. The cal record info is also returned in the cals structure.
The header will be that from the first record of the on with the
following modifications:
h.std.gr,az,chttd is the average for the on scan
h.cor.lag0pwratio is the average of (on-off)/off the units will be
TsysPosOff or kelvins depending on /sclcal.
If bonrecs is specified then the individual on records will also
be returned. If boffrecs is specified then the individual off records
will be returned. The units for these spectra will be Kelvins
if /sclcal, Jy if /scljy is set or correlator units (linear in power).
structures:
b - holds the spectral data and headers: has N elements:
b.b1 - board 1 thru
b.bn - board n
each b.bN has:
b.bN.h - header
b.bN.d[nlags,nsbc] - spectral data
t - is a structure that holds temperature/power data .
For each element the first dimension holds the polarizations.
If there are two polarizations per board then polA is in t.aaa[0,x] and
and polB is in t.aaa[1,x]. If there is only 1 polarization per
board then the data is always in t.aaa[0,x].
The second dimension is the board.
t.src[2,8] (on-off)*X power
t.on[2,8] on power
t.off[2,8] off power
The power units will be Kelvins if /corscl is set. If not,
the units will be TsysOffPos.
cals[nbrds] holds the cal info. it is an array of {cal} structures:
cals[i].h - header from first cal on.
cals[i].calval[2] - cal value for each sbc this board
cals[i].calscl[2] - to convert corunits to kelvins sbc,1,2
cals[i].h.cor.calOn - has the total power calOn
cals[i].h.cor.calOff - has the total power calOff .. units are
correlator units.
NOTES:
1. If the individual records are returned and /sclJy is returned then
the values will be the SEFD (system equivalent flux density).
No bandpass correction is done on these individual records.
values for the ons,offs with be Tsys/gain.
2. There is a difference between < on/off - 1> and
/ namely the bandpass shape .. let f be the bandpass shape:
then < f*on/(f*off) -1> does not have the bandpass
but / does the basic problem is that
/ is not equal to
(See /pkg/rsi/local/libao/phil/Cor2/corposonoff.pro)
NAME:
corposonoffm - process a on/off scan using a mask
SYNTAX: istat=corposonoffm(lun,m,b,t,cals,sclcal=sclcal,scan=scan,sl=sl,$
maxrecs=maxrecs,dbg=dbg,swappol=swappol)
ARGS :
lun: int .. assigned to open file.
m: {cormask} holds the mask in the bi.data portion.
should be 1 or 0.
Note: if m.b1[1024,2] has 2 masks per board, this
routine uses the first maskfor both pols..
KEYWORDS:
sclcal: if not zero then scale to kelvins using the cal on/off
that follows the scan
scan: position to first record of scan before reading.
sl[]: {sl} returned from getsl(). If provided then direct access
is available to scan.
maxrecs: maximum records in a scan. default 300.
swappol: if set then swap polA,polB cal values. This can be
used to correct for the 1320 hipass cable reversal or
the use of a xfer switch in the iflo.
_extra : calval=caltouse[2,nbrds] .. for corcalonoff..
RETURNS:
b: {corget} return (on-off)*X here..
X is:
sclCal true : (kelvins/corcnt)/.. units are then K
sclcal false: off .. units are TsysOff
t: {temps} total power temp info contains..
t.src .. on/off-1 (same as lag0pwrratio)
t.on .. temp on src
t.off .. temp off
units are Kelvins if /sclcal else Tsys Off
cals[]: {cal} return cal info. see coronoffcal..
this is an optional argument
istat: int 1 gotit, 0 didnot complete
cal value to use polA,B by board instead of
value in file
DESCRIPTION:
Process a position switch on,off pair. Return a single spectrum
per sbc of (posOn-posOff)/posOff. The units will be Tsys posOff or
Kelvins (if /sclcal) is set. The t structure holds the total power info
for on/off -1, on source, and off source. The units are either
Kelvins (/usecal set) or TsysOff.
The cal record info is also returned in the cals structure.
The header will be that from the first record of the on with the
following modifications:
h.std.gr,az,chttd is the average for the on scan
h.cor.lag0pwratio is the average of (on-off)/off the units will be
TsysPosOff or kelvins depending on /sclcal.
structures:
b - holds the spectral data and headers: has N elements:
b.b1 - board 1 thru
b.bn - board n
each b.bN has:
b.bN.h - header
b.bN.d[nlags,nsbc] - spectral data
t - is a structure that holds temperature/power data .
For each element the first dimension holds the polarizations.
If there are two polarizations per board then polA is in t.aaa[0,x] and
and polB is in t.aaa[1,x]. If there is only 1 polarization per
board then the data is always in t.aaa[0,x].
The second dimension is the board.
t.src[2,8] (on-off)*X power
t.on[2,8] on power
t.off[2,8] off power
The power units will be Kelvins if /corscl is set. If not,
the units will be TsysOffPos.
cals[nbrds] holds the cal info. it is an array of {cal} structures:
cals[i].h - header from first cal on.
cals[i].calval[2] - cal value for each sbc this board
cals[i].calscl[2] - to convert corunits to kelvins sbc,1,2
cals[i].h.cor.calOn - has the total power calOn
cals[i].h.cor.calOff - has the total power calOff .. units are
correlator units.
Note.. there is a difference between < on/off - 1> and
/ namely the bandpass shape .. let f be the bandpass shape:
then < f*on/(f*off) -1> doesnot have the bandpass
but / does the basic problem is that
/ is not equal to
(See /pkg/rsi/local/libao/phil/Cor2/corposonoffm.pro)
NAME:
corposonoffrfi - process a position switch pair with rfi excision.
SYNTAX: istat=corposonoffrfi(lun,bonoff,calScl,bonrecs,boffrecs,$
scan=scan,maxrecs=maxrecs,sl=sl,$
han=han,sclcal=sclcal,sclJy=sclJy,$
dataPassedIn: boninp=boninp,boffinp=boffinp,bcalinp=bcalinp,$
smorfi=smorfi,
rfiExcParms: frqChnNsig=frqChnNsig,tpNSig=tpNsig,$
flatTsys=flatTsys,adjKer=adjKern,$
ReturnRfiInfo: bsmpPerChn=bsmpPerChn,bresdat=bresdat,bmask=bmask,$
bsecPerChn=bsecPerChn,brmsByChn=brmsByChn,$
ReturnRfiTpInfo: tpdat=tpdat,tpmask=tpmask,$
verbtp=verbtp,verbwt=verbwt,verball=verball
INPUT ARGS:
lun: int file lun if program is to read the data from disc.
INPUT KEYWORDS:
POSITIONING:
scan: long scan number of on scan. If not supplied then start
reading from the current position.
maxrecs: long If then number of records in a scan in more than
300, then you need to specify this value using the
maxrecs keyword. If not, the routine will only
read the first 300 records.
sl[]: {getsl} If you first scan the file with getsl() then passing
the sl[] array in will allow the routine to do
random access (rather than sequentially searching
for the scan on interest).
DATA SCALING:
han: If set then hanning smooth the data on input.
sclCal: if set then scale the data to kelvins using the
cal scan that follows the position on/off.
By default the units are Tsys (since we divide
by to off source).
sclJy: If set then scale the data to Janskies. This is done
by first scaling to kelvins using the cal on/off
scans and then using the gain curve to scale from
kelvins to Janskies. This only works for those
receivers that have gain curves.
cal scan that follows the position on/off.
useMask: If true then use the bmask passed in as a starting
point. This lets you call the routine once with one set of
parameters. A second call with a different set of
parameters can then be made with the first calls
mask included in the final mask.
USERINPUTS DATA:
boninp[]: {corget} The user passes in the on source records rather than
reading them from disc.
boffinp[]:{corget} The user passes in the off source records rather than
reading them from disc.
bcalinp[]:{corget} The user passes in the cal on/off records rather than
reading them from disc.
RFI EXCISION:
smorfi: long When searching for rfi, smooth the data by
smorfi channels before searching. smorfi should
be an odd number. This is only used for searching
for the rfi, the returned data is not smoothed.
flatTsys: If set, then divide each spectra by
its median value before searching each freq channel
along the time direction. This will allow records with
total power fluctuations to be included.
adjKer[i,j]: float convol this array with the mask array for each
sbc. If the convolved value is not equal to the
sum of adjKer then also exclude this point. This
array lets you exclude points adjacent to bad
points (xdim is frequency, ydim is time). The
array dimensions should be odd.
frqChnNsig: float The number of sigmas to use when excluding a point.
linear fits are done by channel along the time axis.
Points beyond frqChnNsig will be excluded from the
averaging.
tpNsig: float The total power is computed vs time for the points
that pass the freqchannel test. A linear fit is
then done versus time. Any time points whose
residuals are greater than tpNSig will be ignored.
The default is 0 (this computation is not done)
rowThrLev: float Row threshhold. If the fraction of good data points
(non zero mask value) in a row is less than rowThrLev
then the entire row will not be used in the average.
eg rowthrlev=.8 means that any rows with fewer than
80% good points will not be included in the average.
The default is to use all good points in the mask.
rowThrDat[] will return each rows fraction of good pnts
RETURNED ARGS:
bonoff: {corget} return processed (on-off)/off here. The units will
by Tsys,K, or K/Jy depending on the sclCal and
sclJy keywords.
calScl[2,nbrds]:float values used to scale correlator counts to kelvins.
bonrecs[]:{corget} return individual on records. This is an
optional argument. The units will be K (if /sclcal),
Jy if /sclJy else correlator units
boffrecs[] : {corget} return then individual off records. This is an
optional argument. The units will be K (if /sclcal),
Jy if /sclJy else correlator units
RETURNED KEYWORDS:
bsmpPerChn: {corget} holds the number of samples averaged for each
freq channel.
bsecPerChn: {corget} holds the number of secs averaged for each
freq channel. same as bsmpperchn but scaled to time
brmsByChn: {corget} rms/mean for each channel. This is before
adjKer,rowThr, or exclusion by total power.
bresdat[nsmp]:{corget} holds the residuals (in units of sigma along a
frequency channel) for all the input data.
bmask[nsmp]:{corget} hold mask for points used. 0 excluded, 1. included.
tpdat[ntmsmp,2,nbrds]:float Total power vs time for each subband used
for total power clipping. This is after the
bad frequency channels are removed. If only one
polarization per board, then the 2nd pol index will
have 0's.
tpmask[ntmsmp,2,nbrds]:float Total power mask vs time for each subband used
will contain 1 if time sample used, 0 if it
was not used.
rowThrDat[ntmsmp,2,nbrds]:float Fraction of good points in each row.
You can look at this to see where the rowthr
should be set. This array is returned if the
keyword rowThrLev is used.
DESCRIPTION:
Process a position switch on,off pair. Return a single spectrum
per sbc of (posOn-posOff)/posOff. Try to excise rfi before averaging.
The units will be:
Kelvins if /sclcal is set
Janskies if /scljy is set .
TsysUnits if none of the above is set.
If bonrecs is specified then the individual on records will also
be returned. If boffrecs is specified then the individual off records
will be returned. The units for these spectra will be Kelvins
(if /sclcal,/scljy is set) or correlator units (linear in power).
The user can optionally pass in the data rather than reading it
from disc (using the keywords boninp,boffinp,bcalinp). If you pass
in data, you must pass in all the data (ons,offs, and cals).
The processing is:
1. Input all of the on records. Hanning smooth them if requested.
2. Input all of the off records. Hanning smooth them if requested.
3. Input the two calon, caloff records if we are scaling to kelvins, Jy,
or the user supplied the calscl keyword.
4. Compute onoffAr[nsmp]= on/off for each record.
5. For each sbc of each correlator board call bpccmpbychan(onoffAr).
bpcbychan() does the following:
0. If usemask is set then initialize bmask starts with the values
passed in. If usemask is not set then bmask is initialized to all 1's.
a.For each freq chan iteratively do a linear fit along the time direction
throwing out any points that have fit residuals greater than freqChnNsig
(the default is 3 sigma). Continue fitting the same channel until no
new points are thrown out (or an iteration limit is hit).
b.Flag all of the good and bad points found for a channel in bmask:
good values get 1, bad values get 0. Also store the fit residuals in
bresdat[].
6. If tpNsig is provided, then compute the total power in the onoffar[]
for each time sample using the "good" data points. Do an interative
linear fit throwing out complete time points if the fit residuals for
a time point is greater than tpNsig. Any time points that get flagged
as bad will have all of their frequency channels flagged as bad.
This processing is only done if the keyword tpNsig is provided.
7. If rowThrLev is supplied then for each time sample (row) compute the
number of good points in the bmask If this average is
less than rowThrLev then mark every point in that row of the
mask as bad (set it to 0). It will not be included in the average.
7. AFter flagging all of the good and bad points, average onoffAr[] over
the good points and then subtract 1 (on-off)/off = (on/off - 1). Store
the number of good points found in each frequency channel in the array
bSmpPerChn.
8. If we use the cals for scaling, also compute the average off source
value over the "good" data points. This is needed since we've been
working the whole time with on/off.
avgOffSrc
9. If scaling to Kelvins,or Jy call corcalcmp(calrecs, /blda) with the cal
on,off records. The routine will:
a. compute calRatio=(calOn-calOff)/caloff.
b. call corblauto(calRatio,deg=1,mask=calmask). It will interatively
fit a linear function to (calratio) throwing out any points that
are above 3 sigma. The cal mask will be all of the frequency channels
that fall below 3 sigma.
c. compute corToK= calInK/avg(calon-caloff) averaging over the "good"
cal points.
-We've compute on/off which is in units of Tsys. The conversion factor
computed in c. converts from correlator counts to Kelvins. We need to
rescale the onoff computed in 7. back to correlator counts. This should
be done using the same frequency channels that were used to compute the
cal scaling factor.
meanTpoff = mean(avgOffSrc[over channels where calmask is 1])
- compute onoff*meanTpOff * corToK
10. If scaling to Jy then get the gain for this frequency and az,za in
K/Jy. Divide each frequency point in onoff by this value to get Jy.
HOW WELL DOES THE EXCISION WORK:
1. For a signal to be identified as rfi it must vary in time relative
to the median value in the channel. If the signal is stable or slowly
drifts by a single channel over the integration time, then it will
not be excised.
2. We typically take 1 second records for 300 seconds on and off source.
This routine fits the 1 second data and uses the 1 second sigmas to try
and find bad data points.
It then averages over these good data points by up to 300 which ;
reduces the sigma by sqrt(300). So a weak signal will pass the test.
the smorfi keyword will smooth over frequency channels before searching
for the bad points. This may help (depending on the width of the rfi).
EXAMPLES:
Suppose we've opened a data (lun points to it) and it has the
following data:
W49A 417600054 300 onoff on 19:26:12 11
W49A 417600055 300 onoff off 19:32:14 11
W49A 417600056 1 calonoff on 19:37:17 11
W49A 417600057 1 calonoff off 19:37:28 11
istat=corposonoffrfi(lun,bonoff,calscl,scan=417600054L,/han,/sclcal,$
bsmpPerchn=bsmpPerchn,bmask=bmask)
bonoff[1] will then have on/off-1 in Kelvins
bsmpPerchn[1] will have the number of samples used in each frq channel
bmask[300] will have a 1,0 for the good and bad data points.
; plot on/off -1 in kelvins
corplot,bonof
;plot the good points per channel
corplot,bsmpperchn
; make an image of freq vs time of the points that were used (bright)
; and not used (dark).
; (first do xloadct to load a linear ramp in the color table)..
img=corimgdisp(bmask,/nobpc,corplot)
12nov04 - bonrecs,boffrecs were always being returned in correlator
units... changed so they are returned in Cor,Tsys,K , or Jy.
21mar05 - added secsPerChn
28jul06 - added rowThrLev, rowThrDat
added usemask
removed badNsig
(See /pkg/rsi/local/libao/phil/Cor2/corposonoffrfi.pro)
NAME:
corposonoffrfisep - process a position switch pair with rfi excision.
this version searchs on,off separately rather than
on/off -1
SYNTAX: istat=corposonoffrfisep(lun,bonoff,calscl,bonrecs,boffrecs,$
scan=scan,maxrecs=maxrecs,sl=sl,$
han=han,sclcal=sclcal,sclJy=sclJy,$
dataPassedIn: boninp=boninp,boffinp=boffinp,bcalinp=bcalinp,$
smorfi=smorfi,
rfiExcParms: frqChnNsig=frqChnNsig,tpNSig=tpNsig,badNsig=badNsig,$
flatTsys=flatTsys,adjKer=adjKern,$
ReturnRfiInfo: bsmpPerChn=bsmpPerChn,bresdat=bresdat,bmask=bmask,$
ReturnRfiTpInfo: tpdat=tpdat,tpmask=tpmask,$
verbtp=verbtp,verbwt=verbwt,verball=verball
INPUT ARGS:
lun: int file lun if program is to read the data from disc.
INPUT KEYWORDS:
POSITIONING:
scan: long scan number of on scan. If not supplied then start
reading from the current position.
maxrecs: long If then number of records in a scan in more than
300, then you need to specify this value using the
maxrecs keyword. If not, the routine will only
read the first 300 records.
sl[]: {getsl} If you first scan the file with getsl() then passing
the sl[] array in will allow the routine to do
random access (rather than sequentially searching
for the scan on interest).
DATA SCALING:
han: If set then hanning smooth the data on input.
sclCal: if set then scale the data to kelvins using the
cal scan that follows the position on/off.
By default the units are Tsys (since we divide
by to off source).
sclJy: if set then scale the data to Janskies. This is done
by first scaling to kelvins using the cal on/off
scans and then using the gain curve to scale from
kelvins to Janskies. This only works for those
receivers that have gain curves.
cal scan that follows the position on/off.
useMask: If true then use the bmask passed in as a starting
point. This lets you call the routine once with one set of
parameters. A second call with a different set of
parameters can then be made with the first calls
mask included in the final mask.
USERINPUTS DATA:
boninp[]: {corget} The user passes in the on source records rather than
reading them from disc.
boffinp[]: {corget} The user passes in the off source records rather than
reading them from disc.
bcalinp[]: {corget} The user passes in the cal on/off records rather than
reading them from disc.
RFI EXCISION:
smorfi: long When searching for rfi, smooth the data by
smorfi channels before searching. smorfi should
be an odd number. This is only used for searching
for the rfi, the returned data is not smoothed.
flatTsys: If flatTsys is set, then divide each spectra by
its median value before searching each freq channel
along the time direction. This will allow records with
total power fluctuations to be included.
adjKer[i,j]: float convol this array with the mask array for each
sbc. If the convolved value is not equal to the
sum of adjKer then also exclude this point. This
array lets you exclude points adjacent to bad
points (xdim is frequency, ydim is time). The
array dimensions should be odd.
frqChnNsig: float The number of sigmas to use when excluding a point.
linear fits are done by channel along the time axis.
Points beyond frqChnNsig will be excluded from the
averaging.
tpNsig: float The total power is computed vs time for the points
that pass the freqchannel test. A linear fit is
then done versus time. Any time points whose
residuals are greater than tpNSig will be ignored.
The default is 0 (this computation is not done)
rowThrLev: float Row threshhold. If the fraction of good data points
(non zero mask value) in a row is less than rowThrLev
then the entire row will not be used in the average.
eg rowthrlev=.8 means that any rows with fewer than
80% good points will not be included in the average.
The default is to use all good points in the mask.
rowThrDat[] will return each rows fraction of good pnts
RETURNED ARGS:
bonoff: {corget} return processed (on-off)/off here. The units will
by Tsys,K, or K/Jy depending on the sclCal and
sclJy keywords.
calScl[2,nbrds]:float values used to scale correlator counts to kelvins.
bonrecs[] : {corget} return individual on records. This is an
optional argument. The units will be K (if /sclcal) or
correlator units.
boffrecs[] : {corget} return then individual off records. This is an
optional argument. The units will be K (if /sclcal) or
correlator units.
RETURNED KEYWORDS:
bsmpPerChn: {corget} holds the number of samples averaged for each
freq channel.
bresdat[nsmp]:{corget} holds the residuals (in units of sigma along a
frequency channel) for all the input data.
bmask[nsmp]:{corget} hold mask for points used. 0 excluded, 1. included.
it has anded the on and off source masks.
tpdat[ntmsmp,2,nbrds]:float Total power vs time for each subband used
for total power clipping. This is after the
bad frequency channels are removed. If only one
polarization per board, then the 2nd pol index will
have 0's.
tpmask[ntmsmp,2,nbrds]:float Total power mask vs time for each subband used
will contain 1 if time sample used, 0 if it
was not used.
rowThrDat[ntmsmp,2,nbrds]:float Fraction of good points in each row.
You can look at this to see where the rowthr
should be set. This array is returned if the
keyword rowThrLev is used.
DESCRIPTION:
Process a position switch on,off pair. Return a single spectrum
per sbc of (posOn-posOff)/posOff. Try to excise rfi before averaging.
The units will be:
Kelvins if /sclcal is set
Janskies if /scljy is set .
TsysUnits if none of the above is set.
If bonrecs is specified then the individual on records will also
be returned. If boffrecs is specified then the individual off records
will be returned. The units for these spectra will be Kelvins
(if /sclcal,/scljy is set) or correlator units (linear in power).
.. note the current units for onrecs, offrecs is Tsys.. i need to
fix this..
The user can optionally pass in the data rather than reading it
from disc (using the keywords boninp,boffinp,bcalinp). If you pass
in data, you must pass in all the data (ons,offs, and cals).
The processing is:
1. Input all of the on records. Hanning smooth them if requested.
2. Input all of the off records. Hanning smooth them if requested.
3. Input the two calon, caloff records if we are scaling to kelvins, Jy,
or the user supplied the calscl keyword.
5. For each sbc of each correlator board call bpccmpbychan(onrecs) and
bpccmpbychan(offrecs)
bpcbychan() does the following:
a.For each freq chan iteratively do a linear fit along the time direction
throwing out any points that have fit residuals greater than freqChnNsig
(the default is 3 sigma). Continue fitting the same channel until no
new points are thrown out (or an iteration limit is hit).
b.Flag all of the good and bad points found for a channel in bmask:
good values get 1, bad values get 0. Also store the fit residuals in
bresdat[].
6. And the masks created for the on recs and the off recs. Use this
mask for the combined dataset. This means that if an on or an off
recs i bad, then that data point will not be used.
7. If tpNsig is provided, then compute the total power in the onrecs,offrecs
arrays for each time sample using the "good" data points. Do an interative
linear fit throwing out complete time points if the fit residuals for
a time point is greater than tpNsig. Any time points that get flagged
as bad in the on or the off will have all of their frequency channels
flagged as bad.
This processing is only done if the keyword tpNsig is provided.
8. Compute bonavg,boffavg averaging over the good data points.
Then compute bonoff=(bonavg/boffavg - 1).
9. If scaling to Kelvins,or Jy call corcalcmp(calrecs, /blda) with the cal
on,off records. The routine will:
a. compute calRatio=(calOn-calOff)/caloff.
b. call corblauto(calRatio,deg=1,mask=calmask). It will interatively
fit a linear function to (calratio) throwing out any points that
are above 3 sigma. The cal mask will be all of the frequency channels
that fall below 3 sigma.
c. compute corToK= calInK/avg(calon-caloff) averaging over the "good"
cal points.
-We've compute on/off which is in units of Tsys. The conversion factor
computed in c. converts from correlator counts to Kelvins. We need to
rescale the onoff computed in 8. back to correlator counts. This should
be done using the same frequency channels that were used to compute the
cal scaling factor.
meanTpoff = mean(boffavg[over channels where calmask is 1])
- compute bonoff*meanTpOff * corToK
10. If scaling to Jy then get the gain for this frequency and az,za in
K/Jy. Divide each frequency point in onoff by this value to get Jy.
SOME OF THE DIFFERENT WAYS TO USE THE ROUTINE:
The simplest way is to call it with all of the defaults. It will use
3 sigma as the clipping level.
Since final onoff spectra will be averaged over the individual records,
its rms will be sqrt(nrecs) better than the individual records. We are
using the rms on the individual records to try and find rfi. There can
be rfi that slips under the individual rfi rms but sticks out in the
averaged onoff spectra. The smorfi keyword can be used to improve this
situation (depending on the signature of the rfi). It will smooth smorfi
adjacent frequency channels before searching along the time axis for
rfi. This will improve the individual record rms's by sqrt(smorfi).
Another way to catch weak rfi is to use the tpnsig option. This will compute
the total power at each time sample using only the "good" data points
(found by searching along the time axis for each frequency channel). Any
total power points that stick up more than tpNsig from a linear fit to
tp will have the entire time sample discarded.
Often times rfi will be strong in a few time or frequency channels and
weaker in some adjacent ones. The adjker (Adjacent Kernel) keyword lets
you flag data points next to known bad points as bad. You construct
a 2-d array (1st dimension is freq, 2nd dimension is time). Fill the
array with 1's and 0's. The mask of good data points is computed by
searching along the time axis for each frequency channel and placing a
1 at the good points and a 0 at the bad points. The adjKer is then
convolved with the mask. Any points after convolution that do not equal
the sum of the adjker will be marked as bad. If you center the adjker
on a data point, then all data points under the non-zero elements of the
adjker must be 1. If not the center point is marked bad.
To mark frequency points next to bad points as bad use:
adjker=fltarr[3,1]+1
To mark adjacent time points as bad use:
adjker=fltarr[1,3]+1
To mark adjacent time and freq points as bad (not including diagonols use:
adjker=fltarr(3,3)+1
ind=[0,2,6,8]
adjker[ind]=0
EXAMPLES:
Suppose we've opened a data (lun points to it) and it has the
following data:
W49A 417600054 300 onoff on 19:26:12 11
W49A 417600055 300 onoff off 19:32:14 11
W49A 417600056 1 calonoff on 19:37:17 11
W49A 417600057 1 calonoff off 19:37:28 11
istat=corposonoffrfsep(lun,bonoff,calscl,scan=417600054L,/han,/sclcal,$
bsmpPerchn=bsmpPerchn,bmask=bmask)
bonoff[1] will then have on/off-1 in Kelvins
bsmpPerchn[1] will have the number of samples used in each frq channel
bmask[300] will have a 1,0 for the good and bad data points.
; plot on/off -1 in kelvins
corplot,bonof
;plot the good points per channel
corplot,bsmpperchn
; make an image of freq vs time of the points that were used (bright)
; and not used (dark).
; (first do xloadct to load a linear ramp in the color table)..
img=corimgdisp(bmask,/nobpc)
(See /pkg/rsi/local/libao/phil/Cor2/corposonoffrfisep.pro)
NAME:
corpwr - return power information for a number of recs
SYNTAX: nrecs=corpwr(lun,reqrecs,pwra,lasthdr,pwrcnt=pwrcnt)
ARGS:
lun - file open to
reqrecs -requested records to return
KEYWORDS:
pwrcnt if set the return power counter data (50Mhz) rather than
0 lag. This is measured before the digitizer with a
vTof converter (it is also not blanked).
RETURNS:
pwra - returns an array pwra[nrecs] of {corpwr} struct
nrecs - number of recs found, 0 if at eof, -1 if hdr alignment/io error
DESCRIPTION:
Return the total power information for the requested number of
records. The data is returned in the array pwra. Each element of the
array contains:
pwra[i].scan - scan number
pwra[i].rec - record number
pwra[i].time - seconds from midnight end of record.
pwra[i].nbrds- number of boards in use
pwra[i].az - az (deg) end of this record.
pwra[i].za - za (deg) end of this record.
pwra[i].azErr- az (asecs great circle) end of this record.
pwra[i].zaErr- za (asecs great circle) end of this record.
pwra[i].pwr[2,4] - total power info. first index in pol1,pol2
2nd index are 4 correlator boards.
There will only be valid data in the first pwra[i].nbrds entries of
pwra.
In pwra[i].[i,j], i=0,1 will be pola, polb if two polarizations were
recorded in the board. If only one polarization was recorded on the
board, then i=0 holds the data (either pola,polB) and i=1 has no data.
(See /pkg/rsi/local/libao/phil/Cor2/corpwr.pro)
NAME:
corpwrfile - input the power information from the entire file.
SYNTAX: nrecs=corpwrfile(filename,pwra,maxrecs=maxrec)
ARGS
filename: string.. filename to read from.
KEYWORDS
maxrecs : long max number of recs to read in. def:20000
RETURNS
pwra - returns an array pwra[nrecs] of {corpwr} struct
nrecs - number of recs found, 0 if at eof, -1 if hdr alignment/io error
DESCRIPTION:
This routine opens the file, calls corpwr() and then closes the
file. See corpwr() for a description of the power structure.
You only need to use the maxrecs keyword if the file has more than
20,000 records.
(See /pkg/rsi/local/libao/phil/Cor2/corpwrfile.pro)
NAME:
corradecj - get ra,dec J2000 positions
SYNTAX: npnts=corradecj(bscan,raHrJ,decDegJ)
ARGS:
bscan[n]: {corget} get ra,dec positions for this scans data
KEYWORDS:
RETURNS:
RAHRJ[N]: double ra J2000 in hours for records in scan
decDegI[N]: double dec J2000 in degrees for records in scan
DESCRIPTION:
Take the requested ra/dec positions from the header and interpolate
them to the center of the data samples. If bscan only has a single
record then no interpolation is done.
NOTE: If the data was taken with a mapping routine ;(cormap,cordrift,
etc) you should probably use cormapinp().
(See /pkg/rsi/local/libao/phil/Cor2/corradecj.pro)
NAME:
correcinfo - print record info
SYNTAX: correcinfo,b,inphdr=inphdr,outfd=outfd
ARGS :
b: {corget} record to print info from
KEYWORDS:
inphdr: if set, then b is a header rather than a correlator
record: eg b.b1.h rather than b
outfd : int if provided, then write the information to this
file descriptor (logical unit number) instead of
the terminal. You should open the desctripor with
openw,outfd,filename,/get_lun before calling this
routine.
DESCRIPTION:
Print Info about the correlator rec passed in.
An example output is:
IDL> correcinfo,b
scan:102000270 ra:16:45:42 dec:00:20:58 rcv: 5 azErr: -3.9 zaE: 0.4Asec
yymmdd: 10120 dayno: 20 astTm:09:42:35 az,gr,ch: 386.024 17.800 8.834
recNum: 1 sec/Rec:60
brd cfr bw nchn lvl npol lagCf CorDb lag0
1 1612.2 0.39 1024 9 2 9 11 11 1.27 1.52
2 1666.0 6.25 1024 9 2 9 10 12 1.49 1.42
3 1665.4 0.39 1024 9 2 9 8 9 1.31 1.33
4 1667.3 0.39 1024 9 2 9 11 9 1.37 1.41
src:J1645+021 pattern:onoff
on nrecs: 1 cal:1 movTmSec: 60
(See /pkg/rsi/local/libao/phil/Cor2/correcinfo.pro)
NAME:
corrms - compute rms/Mean by channel
SYNTAX: brms=corrms(bin,rembase=rembase)
ARGS:
bin[] : {corget} array of corget structures
KEYWORDS:
rembase : If set then remove a linear baseline (by channel) before
computing the rms.
RETURNS:
brms : {corget} where brms.bN.d[] is now the rms/mean by channel.
DESCRIPTION
The input bin should be an array of {corget} structures. corrms will
compute the rms/Mean by channel for each frequency bin.
if bin[] has 300 records with 2048 freq channels in brd 1 then:
brms.b1.d[i] is rms(bin[0:299].b1.d[i])/mean(bin[0:299].b1.d[i]
where the mean and rms is computed over the 300 records for each
frequency channel.
(See /pkg/rsi/local/libao/phil/Cor2/corrms.pro)
NAME:
corsavbysrc - create arrays by source.
SYNTAX - istat=corsavbysrc(savelistinp,savefileout,noavg=noavg,$
nosingle=nosingle,/filelist)
ARGS:
savelistinp[]: string list of save files to use as input
These files should contain:
bf[n]: an array of {corget} structures.
savefileout: string name of save file to store the arrays by source.
KEYWORDS:
noavg: if set then do not bother to save averages. This may be
needed if you end up with too many variables to save at once.
nosingle: if set then do not bother to save the arrays of individual sources.
(it will save b_all and bavg_ if noavg is not set)
filelist: if set then savelistinp is the data filelist used to input the
data. In this case this routine will generate the standard
save filenames to read:
corfile.ddmonyy.projid.n -> ddmonyy.projid.n.sav
It will assume that the .sav files are in the current directory.
DESCRIPTION:
Process a set of idl save files each containing an array bf[n] of
corget structures. This routine will input all of the save files into a single
b_all[m] array and then create separate arrays by source name.
The new array names will be the source name with the following
modifications:
1. Each source name will be prepended with
b_ for the regular data.
bavg_ for the averaged data.
2. The following replacements are made within the src name:
+ -> p
- -> m
This is needed to make the source names legal variable names in idl.
For each source array also create an average array that averages
over all of the entries for a single source and over polarization
(there is no weighting for gain or tsys in the averaging).
Save the new set of arrays in idl save format to the file specified by
the savefileout parameter. This is the simplest way to pass a variable
number of arguments back to the main routine.
This routine was written to combine datasets output by
multiple calls to pfposonoff(). pfposonoff() creates a save file bf[n]
containing all of the processed on/off-1 position switch scans in a datafile.
corsavbysrc() can then be used to combine these multiple files (days)
so there is a single set of arrays by source. You can then edit/average
each source array giving the final result.
The /noavg keyword will not save the average of each source.
The /nosingle will not store the array for each individual source
(before averaging).
The array b_all() containing all of the sources (before averaging) is
always saved. If you have many sources (gt say 150) then the save
command may fail. In that case you should consider not saving the
averages or the single arrays.
EXAMPLE:
Suppose project a1721 generated the following corfiles (either in
/share/olcor or /proj/a1721). You could process all the data with the
following code:
; The pfposonoff processing:
flist=['corfile.08may03.a1721.2',$
'corfile.08may03.a1721.3',$
'corfile.09may03.a1721.1']
dir=['/share/olcor/','/proj/a1721/'] ; directories to look for the files
voff=.03 ; plotting increment
han=1 ; hanning smooth
hor
ver
for i=0,n_elements(flist)-1 do begin &$
npairs=pfposonoff(flist[i],bf,tf,calsf,han=han,dir=dir,/scljy,/sav) &$
; plot the data so we can see what's up. offset each strip by voff
ver,-voff,(npairs+1)*voff &$
corplot,bf,off=voff &$
endfor
; The processed on/off-1 are now in the save files (one for each corfile)
; Use corsavbysrc() to:
; 1. input all the data from the daily save files.
; 2. create two arrays for each source:
; b_srcname, bavg_srcname which have the individual on/off-1
; and the average for the source. Some characters in the source
; name will be changed: + -> p - -> m
; 3. save all of these arrays in the savbysrc filename.
; 4. the /filelist keyword flags that the string array has the
; original datafilenames and this routine should generate the
; input save file names.
;
savefileout='savebysrc.sav' ; place to store the new data
istat=corsavbysrc(flist,savefileout,/filelist) ;this does all the work
;
; to load the data by source into memory:
restore,savefileout,/verbose
; the averaged data:
BAVG_NGC3384 STRUCT = -> Array[1]
BAVG_NGC3389 STRUCT = -> Array[1]
BAVG_NGC3412 STRUCT = -> Array[1]
BAVG_NGC3489 STRUCT = -> Array[1]
BAVG_NGC3607 STRUCT = -> Array[1]
BAVG_U4385 STRUCT = -> Array[1]
BAVG_U5326 STRUCT = -> Array[1]
BAVG_U6018 STRUCT = -> Array[1]
; the data before averaging. each entry is 1 processed on,off scan pair
B_NGC3384 STRUCT = -> Array[6]
B_NGC3389 STRUCT = -> Array[1]
B_NGC3412 STRUCT = -> Array[6]
B_NGC3489 STRUCT = -> Array[6]
B_NGC3607 STRUCT = -> Array[3]
B_U4385 STRUCT = -> Array[2]
B_U5326 STRUCT = -> Array[1]
B_U6018 STRUCT = -> Array[1]
; to plot an averaged source:
ver,-.01.01 ; since the galaxy blows up the scale
corplot,BAVG_NGC3384,/vel
WARNING:
The data is save by executing a command:
istat=execute(cmd) where cmd is a string holding the save command and all
of the sources to save. If the length of cmd gets above about 2800 bytes
i've seen the save fail. If you have so many sources, try using the
/noavg keyword to not save the averages.
(See /pkg/rsi/local/libao/phil/Cor2/corsavbysrc.pro)
NAME:
corsclcal - scale spectrum to K using the cals.
SYNTAX: corsclcal,b,calI
ARGS: b[];{corget} data input from corget.. single integration
or an array of integrations
calI[nbrds]:{corcal} structure returned from corcalonoff.
or calI:{corcalpf} struct returned from pfcalonoff
DESCRIPTION
scale spectra in b to Kelvins using the cal scale factors that are
stored in the calI structure.
The calI can come from two sources:
1. you called corcalonoff()
In this case calI[nbrd] will be an array with one entry per board.
2. You called pfcalonoff() . This routine returns a single structure that
has the calI for all the boards
This routine will check to see which type is passed to it.
EXAMPLES:
;1. use corcalonoff.
calonscan=315200009L
istat=corcalonoff(lun,calI,scan=calonscan)
; input position on data
poson=315200007L
istat=corinpscan(lun,bposon,scan=poson)
; scale this data to kelvins using our cals info
corsclcal,bposon,calI
;2. use pfcalonoff
calonscan=315200009L
istat=pfcalonoff(lun,junk,pfcalI,scanOnAr=calonscan,/blda)
; input position on data
poson=315200007L
istat=corinpscan(lun,bposon,scan=poson)
; scale this data to kelvins using our cals info
corsclcal,bposon,pfcalI
;
(See /pkg/rsi/local/libao/phil/Cor2/corsclcal.pro)
NAME:
corsmo - smooth correlator data.
SYNTAX: corsmo,b,bsmo,smo=nchn,savgol=degree
ARGS:
b[n]: {corget} data to boxcar smooth
smo[n]: {corget} return smoothed data here. Note.. if only 1 argument
is passed to the routine then the data is smoothed in
place (it is returned in b).
savgol: int Perform savitzky-golay smoothing filter (see idl savgol
routine). Fit a polynomial of order degree to nchn
at a time (degree must be less than nchn).
KEYWORDS:
smo: int number of channels to smooth. default is 3
DESCRIPTION:
corsmo will smooth the data in the structure b. By default boxcar
smoothing is performed. The savgol keyword will smooth by fitting a polynomial
of order "degree" to "nchn"'s at a time (see idl savgol documentation with
order=0).
If a single argument is passed to the routine, then the smoothed data is
returned in place. If two arguments are passed in (b,bsmo) then the data is
returned in the second argument. If b is an array of records then each one
will be smoothed individually.
The edge points will only be smoothed by the number of available datapoints
on each side (see the idl /edge_truncate keyword in the smooth() or savgol
routine).
EXAMPLE:
print,corget(lun,b)
corsmo,b ; this smooths the data and returns it in b.
print,corget(lun,b)
corsmo,b,bsmo ; this smooths the data and returns it in bsmo
; input an entire scan then smooth each record by 9 channels.
print,corinpscan(lun,b)
corsmo,b,bsmo,smo=9
;
; use polynomial smoothing of order 3 on 31 points at a time.
corsmo,b,bsmo,smo=31,savgol=3
(See /pkg/rsi/local/libao/phil/Cor2/corsmo.pro)
NAME:
corstat - compute mean,rms by sbc
SYNTAX : corstat=corstat(b,mask=mask,/print,/median)
ARGS :
b[n]: {corget} correlator data
KEYWORDS :
mask: {cormask} Compute mean, rms within mask (see cormask).
print: If set then output info to stdout.
median: If set then use median rather than the mean.
RETURNS :
corstat[n]: {corstat} stat info
DESCRIPTION:
corstat computes the mean and rms by sbc for a correlator dataset. The
input data b can be a single {corget} structure or an array of {corget}
structures. The mean,avg will be computed for each record of b[n]. If
a mask is provided then the mean and average will be computed within the
non-zero elements of the mask (see cormask()). The same mask will be used
for all records in b[n].
The returned data structure corstat consists of:
corstat.avg[2,8] the averages
corstat.rms[2,8] the rms's
corstat.fracMask[8] fraction of the bandpass that the mask covered.
a single mask is used for each board.
corstat.p[2,8] This will contain a 1 if pola, 2 if polB and 0 if this
entry is not used.
EXAMPLES:
Process a position switch scan then compute the rms and mean using a mask
that the user defines.
istat=corposonoff(lun,b,t,cals,/sclcal,scan=scan)
cormask,b,mask
cstat=corstat(b,mask=mask,/print)
(See /pkg/rsi/local/libao/phil/Cor2/corstat.pro)
NAME:
corstokes - input and intensity calibrate stokes data.
SYNTAX: istat=corstokes(lun,bdat,bcal,scan=scan,calscan=calscan,$
calinp=calinp,datinp=datinp,avg=avg,maxrecs=maxrecs,$
han=han,sl=sl,edgefract=edgefract,mask=mask,phase=phase
bpc=bpc,fitbpc=fitbpc,smobpc=smobpc,phsmo=phsmo,$
nochk=nochk,mmcor=mmcor,mmret=mmret
ARGS:
lun : int file descriptor to read from
KEYWORDS:
scan : long scan number for data. default is current position
calscan : long scan number for cal on scan. def: scan following
calinp[2]:{corget} pass in the caldata rather than read it from lun
datinp[n]:{corget} pass in the data rather than read it from lun
avg: if set then return the averaged source record
maxrecs: long maximum number of recs to read in. default is 300.
han: if set then hanning smooth the data.
sl[]: {sl} array used to do direct access to scan.
edgefract: float fraction of bandpass on each side to not use during
calibration. default .1
mask:{cormask} mask structure created for each brd via cormask routine
Note: if mask.b1[1024,2] has two masks per board,
this routine uses the first mask for all sbc.
phase: if set then phase calibrate the data
bpc: int 1 band pass correct with cal off
2 band pass correct with calon-caloff
3 band pass correct with mean(data)
4 band pass correct with the median(data)
Last two good for time varying measurements:
mapping,flares, etc..
fitbpc: int fit a polynomial of order fitbpc to the masked
version of the band pass correction and use the
polynomial rather than the data for the "interior
portions of the mask (0 portions of the mask excluding
the outside edges where the filter falls off.
This is only used if bpc =1,or 2.
Note: This ends up only averaging the "interior"
portion of the bandpass( where the line is). Outside
of this, the bandpass has the integration time of the
bpc (say calOn-calOff). 24sep12 i added
polybpc. this fits then entire non-zero portion
of the mask
polybpc: int Do a polynomial fit to chosen bpc (see bpc) over the
non-zero part of the mask. Use this for the entire
bpc. this differs from fitbpc by:
fitbpc - does polynomial fit over non zero part of mask
- then replaces "interior 0 mask" with the
fit. The idea was to mask out the line in the
center of the bp and then replace
just this masked portion with the fit.
The problem is that outside this area,
the measured bpc is used. This normally has
small integration time. Added polybpc
for things like freqsw where you need to
bandpass correct over a larger area.
smobpc: int smooth the bandpass correction by smobpc channels
It should be an odd number of channels. This is only
valid if bpc =1 or 2.
phsmo: int number of channels to smooth the sin,cos of the
phase angle before computing the arctan. default 11.
nochk: if set then don't bother to check if these are valid
cal records. Good to use if data from a non standardprg.
mmcor: if set then apply the mueller matrix correction
to the data (if it exists for the receiver)
mmret[4,4,nbrds]: if mmcor set and mmret is provided, return the mueller
matrix for each board.
RETURNS:
bdat: {corget} intensity calibrated data spectra
bcal: {corget} intensity calibrated cal spectra
istat: 1 ok
: 0 hiteof
:-1 no cal onoff recs
:-2 could not get cal value
:-3 cal,data scans different configs
:-4 at least 1 board did not have stokes data
:-5 sbc length does not match mask length
:-6 illegal bandpass correction requested
DESCRIPTION:
corstokes will intensity calibrate (and optionally phase calibrate)
stokes data given a data scan and cal on,off scans. The data can be read
from disc or input directly to this routine. On output bdat and bcal will be
in units of Kelvins. The /avg keyword will return the average of the
records in the scan after calibration.
If the data is input from disc then lun should be the file descriptor
from the open command. By default it will start reading the data scan from
the current file position. The scan=scan keyword lets you position to the
data scan before reading. By default the calscans will be the two scans
following the data scan. The calscan=calscan keyword lets you position
to the cal on scan before reading them. If the scans on disc have more than
300 records you need to use the maxrecs= keywords so the entire scan will
be input.
By default 10% of the bandpass on each edge is not used for the calibration.
You can increase or decrease this with the edgefract keyword. The mask
keyword allows you to create a mask for each sbc. The calibration will only
use the channels within the mask when computing the gain and phase calibration
factors. You can use this to exclude rfi or spectral lines.
Bandpass correction can be done with the cal off scan or with the
calon-caloff difference spectrum. Since the integration time for the cal is
usually much less than the data integration time, you need to do some
type of averaging to the bandpass so the signal to noise does not increase.
The program lets you fit an N order polynomial to the bandpass with the
fitbpc keyword. An alternative would be to use the smobpc= keyword to
smooth the bandpass.
Phase calibration can be included by using the /phase keyword.
You can pass in the data and/or calscans directly by using the
datinp, calinp keywords.
THE PROCESSING:
Let X and Y be the two polarizations, Px,Py be the total power, and
MN be the correlation of M and N. cosft is a cosine transmform and cacf
is a complex acf from YX and XY. Then the raw data is stored as:
sbc1: I (Px*cosft(XX) + Py*cosft(YY))/(Px+Py)
sbc2: Q (Px*cosft(XX) - Py*cosft(YY))/(Px+Py)
sbc3: U (real(fft(cacf)*Px*Py)/(Px+Py)
sbc4: V (-img(fft(cacf)*Px*Py)/(Px+Py)
(the Q,U,V naming assumes linear polariztion).
The intensity calibration consists of:
1. Scale all of the spectra by (Px+Py) to get unnormalized data.
2. Convert, I,Q back into XX,YY spectra
XX=(I+Q)*.5,YY=(I-Q)*.5
3. Compute the average <calon>,<calOff> over the specified channels.
The specified channels are determined by:
a. The mask for each board from the mask keyword
b. Use edgefract to throw out this fraction of channels at each edge.
c. Use an edgefraction of .1 (10%)
4. The conversion factor Tcal/(<calOn> - <calOff>) is computed for
the two polarizations: scaleXX, scaleYY
5. If band pass correction is done, multiply
scaleXX=scaleXX/normalized(bandpassXX)
scaleYY=scaleYY/normalized(bandpassYY)
bandpassXX can be taken from the calon or calon-caloff. You can
smooth the bpc (smobpc=) or you fit a polynomial to it (fitbpc=). If
fitting is selected then the channels specified in 3. above are used
for the fit.
The normalization of the bandpass is also computed over the channels
selected in 3. above.
6. For the cals and the data compute
sbc1:I = (XX*scaleXX + YY*scaleYY)
sbc2:Q = (XX*scaleXX - YY*scaleYY)
sbc3:U = U*scaleXX*scaleYY
sbc4:V = V*scaleXX*scaleYY
7. The phase correction will do a linear fit to the phase using the
channels selected in 3. above.
When deciding on the mask or edge fraction to use, you should have
a region where the calon-calOff is relatively flat (no filter rolloff and
no rfi).
EXAMPLE:
Suppose we have the following scans:
corlist,lun
SOURCE SCAN GRPS PROCEDURE STEP LST RCV
B1641+173 210200235 5 on on 17:42:08 5
B1641+173 210200236 1 calonoff on 17:47:10 5
B1641+173 210200237 1 calonoff off 17:47:21 5
40.42+0.70 210200238 5 on on 18:02:53 5
40.42+0.70 210200239 1 calonoff on 18:07:56 5
40.42+0.70 210200240 1 calonoff off 18:08:07 5
To process the first two sets:
--rew,lun
--print,corstokes(lun,bdat,bcal); will process the first set
--print,corstokes(lun,bdat,bcal,/han); will process the 2nd set with hanning
To process the 2nd set directly with an edgefraction=.12:
--print,corstokes(lun,bdat,bcal,scan=210200238L,edgefract=.12)
To input the data first, interactively create a mask, and then process
the data with a mask
--print,corinpscan(lun,bdatinp,scan=210200238L,/han)
--print,corgetm(lun,2 ,bcalinp,/han) ; reads the next two records
--cormask,bcalinp,mask ; interactively creates the mask
--print,corstokes(lun,bdat,bcal,calinp=bcalinp,datinp=bdatinp,mask=mask)
Use the same cal for multiple data scans:
--print,corgetm(lun,2 ,bcalinp,scan=210200236L/han);
--print,corstokes(lun,bdat1,bcal1,calinp=bcalinp,scan=210200235L)
--print,corstokes(lun,bdat2,bcal2,calinp=bcalinp,scan=210200238L)
Do amplitude and phase calibration. Use the cal off for the bandpass
correction. Use a 3rd order polynomial fit to the cal off for the bandpass
correction.
--print,corstokes(lun,bdat,bcal,scan=210200238L,/phase,bpc=1,fitbpc=3)
The bandpass correction is a bit tricky and depends on what type of
observations you are doing. The integration time for the off is usually
a lot less than the on positions so you need to use either the bandpass
fit or smoothing. It would probably be a good idea to add an option for
the user to input a bandpass to use for the correction (from an off src
position).
SEE ALSO:
cormask,corlist
(See /pkg/rsi/local/libao/phil/Cor2/corstokes.pro)
NAME:
corstostr - move a structure to an array of structs.
SYNTAX: corstostr,b,ind,barr
ARGS: b[m]: {corget} structure to move
ind: long index (0..n-1) in barr to store b
barr[n]: {corget} array of structures. move to here.
DESCRIPTION:
The data returned by corget is an anonymous structure. This allows
it to be dynamically generated. A drawback is that you can't combine
different anonymous structures into an array. This routine allows
you to place corget data from multiple calls into an array. The
data placed in barr must all be the same correlator configuration
(number of boards, polarizations, and lags).
The input structure can be a single structure or an array of structures.
You must pre-allocate barr with corallocstr().
EXAMPLE:
..storing 1 structure at a time. Assume you want space for numrecs.
numrecs=200L
for i=1,numrecs-1 do begin
istat=corget(lun,b) .. get the next rec
if i eq 0 then bb=corallocstr(b,numrecs) ;;.. if first, allocate bb
corstostr,b,i,bb .. store this rec
endfor
..combine 4 scans separated by 10 scan numbers each. assume 60 recs per
scan. bb will hold 4*60 records when done.
scan=112500481L
for i=0,3 do begin
istat=corinpscan(lun,b) ;;.. get the scan
if i eq 0 then bb=corallocstr(b[0],4*60) ;;.. if first, allocate bb
corstostr,b,i*60,bb ;;.. store in correct slot
endfor
SEE ALSO: corallocstr.
(See /pkg/rsi/local/libao/phil/Cor2/corstostr.pro)
NAME:
corstripch - stripchart recording of total power
SYNTAX: corstripch,lun,maxpnt=maxpnt,v1=v1,v2=v2,mean=mean,$
scan=scan,rec=rec,rdDelay=rdDelay,pltDelay=pltDelay,$
rdStep=rdStep,win=win,pwrcnt=pwrcnt
ARGS:
lun: int to read power from
KEYWORDS:
maxpnt: long max number of points to display at once. default: 1000
v1[2] : float min,max value for top plot
v2[2] : float min,max value for bottom plot (polA-polB)
mean : long if set then remove mean from both plots.
scan : long position to scan before starting. The
default is the current position
rec : long record of scan to position to before starting.
default: if scan present then:1 else current rec.
rdDelay: float number of seconds to wait between reads if no new
data available.
pltDelay: float number of seconds to wait after each plot. Use this
to slowly scan a file that already exists.
rdStep: long max number of points to try and read at a time. Plot
will increment by this amount when you read a file.
pwrcnt: if set then display power counter data rather than
0 lag.
win : long window number to plot in. Default is 0
DESCRIPTION:
The routine makes a stripchart recording of the total power and power
difference polA-polB. It will step through the data file displaying
up to MAXPNT pnts on the screen. When this value is reached, the plot will
start scolling to the left. When the end of file is hit, the routine
will wait rdDelay seconds (default of 1 second) and then try to read
any new data in the file. This allows you to monitor data as it is
being taken. If you are scrolling through a file offline, use PLTDELAY
to slow down the plotting rate (if it is too fast). At the end of the file
hit any key and the enter q to quit (see below).
The top plot is the 0lag total power. The units are linear in power and the
definition is measured/optimum power (for statistics of 3/9level sampling).
You can change the vertical scale with the v1[min,max] keyword or from the
menu displayed you hit any key. The line colors correspond to the 8
subcorrelators available on the interim correlator.
The bottom plot is the power difference PolA-PolB for each correlator
board (NOTE: currently this only works if polA,polB are on the same board).
The vertical scale can be changed with the v2=v2[min,max] keyword or from
the menu displayed when you hit any key.
You can stop the plotting by touching any key on the keyboard.
A menu will be presented that lets you modify some parameters and then
continue. The menu is:
command function
q to quit
r rewind file
v1 min max new min,max for top window
v2 min max new min,max for bottom window
blank line ..continue
You can quit, rewind the file and continue, change the vertical scale of
the top plot with v1 min max, or change the vertical scale of the bottom
plot. Any other inputline will let you continue from where you are
(unlike cormonall, you have to enter return for the program to read the
inputline.
EXAMPLES:
1. monitor the online datataking:
lun=coronl()
corstripch,lun
2. set fewer points per plot for better resolution. Set the top vertical
scale to .8,2. and the bottom plot to -.4,.4.
corstripch,lun,maxpnt=600,v1=[.8,2],v2=[-.4,.4]
3. If you want to restart from the begining of the file:
hit any character
r
and it will rewind an continue
4. If you want to monitor a file offline, with no wait between
updates, and 500 points plotted:
openr,lun,'/share/olcor/corfile.30apr03.x102.1',/get_lun
corstripch,lun,maxpnt=500,v1=[.5,3]
5. Do the same thing but wait 1 second between plots and read 200 points at
a time:
corstripch,lun,maxpnt=500,v1=[.5,3],pltDelay=1,rdstep=200
NOTE:
You can change the size of the plot by expanding or contracting the
plot window. The plot will rescale itself.
(See /pkg/rsi/local/libao/phil/Cor2/corstripch.pro)
NAME:
corsubset - make a copy of data keeping only specified boards
SYNTAX: bret=corsubset(b,brd,pol=pol)
ARGS: b[m]: {corget} original structure
brd[]: int brd's to keep.. count 1,2,3,4
RETURNS:
bret[m] : {corget} return subset data here
if an illegal brd is requested then '' is returned.
You can check this with keyword_set(bret) eq 0.
KEYWORDS:
pol : if set then just return first polarization.
This is used by coravg...
DESCRIPTION:
corsubset will create a subset of a correlator data structure
keeping only the specified boards. It will also update
some header locations so the header will reflect the data.
The input structure can be a single structure or an array of structures.
EXAMPLE:
..Keep boards 1 and 3 of all the records of the data structure.
print,corinpscan(lun,bsum,b,/sum)
b13=corsubset(b,[1,3]) ; all the records
b13sum=corsubset(bsum,[1,3]); the summary record
(See /pkg/rsi/local/libao/phil/Cor2/corsubset.pro)
NAME:
cortblradius - get indices for all positions within radius.
SYNTAX: num=cortlbradius(ra,dec,radius,slar,indlist,hms=hms,radDeg=radDeg,$
encErr=encErr,pattype=pattype,dist=dist)
ARGS :
ra: float/double ra Hours J2000
dec: float/double dec Degrees J2000
radius: float/double radius to search for source. Default Arcminutes
(keyword raddeg switches to degress).
slAr[]: {slcor} array returned from arch_gettbl(../cor)
KEYWORDS:
hms : if set then ra is in hms.s dec is in dms.s rather than hour,degrees
radDeg: if set then radius is in degrees. default is arcminutes.
encErr: float. Average encoder error (Asecs) of each scan selected
must be less than this value. The default is 30asec.
pattype: int restrict the scans to a particular pattern type.
1 - on/off position switch with cal on/off
2 - on/off position switch whether or not cal there
3 - on followed by cal on ,off
4 - heiles calibrate scan two crosses
5 - heiles calibrate scan 1 or more crosses
6 - cal on,off
7 - x111auto with calonoff
8 - x111auto with or without cal
If a pattern type is specified then the returned indlist array
will contain pointers to the first scan of each pattern (
not every scan within the pattern). As an example
pattype=1 would return the indices for the on position scans.
The default is to return all scans independant of pattern
type.
RETURNS:
indList[num]: long indices into slar for positions that match the request.
num: long number of positions found
DESCRIPTION:
Find all of the positions in the SLAR that lie within the specified
RADIUS of the postion RA,DEC. Return the indicies into SLAR in the array
INDLIST. The radius is measured as a true angle on the sky.
The slar positions are the requested ra,dec (rather than the actual).
The routine requires that the average encoder error of each selected
scan be less than 30 asecs (you can change this with the encerr keyword).
EXAMPLE:
The most efficient way to use this routine is to first extract a
subset of the archive into an slar, and then call this routine multiple
times with different radii to see how many sources are available. You
can then view the sources and extract the actual data records from the archive.
.. get indices of all lband data between 1370 and 1420 Mhz taken during
.. 2002.
n=arch_gettbl(010101,021231,slar,slfilear,freq=[1200.,1420],/cor)
... search for observations of 3C286 (J2000 position).
ra =133108.3
dec=303033.0
... search within 5 arcminutes
radius=60.
num=cortblradius(ra,dec,radius,slar,indlist,/hms)
... print the project id's that took this data.
... print the record types (see arch_gettbl for a list).
print,slar[indlist].projid
print,slar[indlist].rectype
... Look for all of the on/off patterns withing 5amin of Virgo A
radius=20.
ra=123049
dec=122328
num=cortblradius(ra,dec,radius,slar,indlist,/hms)
... get the actual data scans processing on/off-1...
nfound=arch_getonoff(slar,slfilear,indlist,bout,infoAr,$
incompat=incompat,/sclJy,/verbose)
... look at an individual on/off-1
corplot,bout[0]
(See /pkg/rsi/local/libao/phil/Cor2/cortblradius.pro)
NAME:
cortpcaltog - compute total power for toggling cals
SYNTAX:nrecs=cortpcaltog(lun,tpon,tpoff,scanst=scanst,frq=frq,hdr=hdr,$
verb=verb,calonInd=calOnind,median=median)
ARGS:
lun int for file to read
KEYWORDS:
scanst: long scan to position to. 0, -1 --> no position
verb : pass to corblauto, -1--> plot fits
ndeg : pass to corblauto. poly fit deg. def=1
fsin : pass to corblauto. sin fit order. def=6
median : if set then use median when averaging
RETURNS:
ncals: long number of cal toggles found (nrecs/2)
tpon[ncals,nsbc,npol] : float total power cal on
tpoff[ncals,nsbc,npol]: float total power cal off
frq[nsbc]: float freq Mhz for center freq each sbc.
hdr[nsbc]: {hdr} return headers from first rec
calOnInd: long 0--> first rec calon, 1--> 2nd rec cal on
DESCRIPTION:
Process a scan where the cal is toggling on/off each record.
Steps:
- input the records of scan with hanning smoothing
- compute bavg=calOnAvg/calOffAvg.. this removes bandpass
- fit 6 order harmonic, 1 order poly to bavg using corblauto.
create mask of all points with residuals < 3 sigma in fit.
- for each record:
compute totalPwr over mask. place in tpon or tpOff
(See /pkg/rsi/local/libao/phil/Cor2/cortpcaltog.pro)
NAME:
corwriteascii - output an ascii dump of the correlator data
SYNTAX: istat=corwriteascii(b,outfilename)
ARGS:
b[n]: {corget} data to write.
outfilename: string name of output file
DESCRIPTION:
corwriteascii will output a correlator dataset in ascii to the specified
filename. Data is printed by scan. Within each scan the data is printed by
board. Each channel of spectral data will appear on a separate line with
channelTopoFreq channelVel pol1Data pol2Data.
Keywords (which start the line with #) are used to separate the scans and
the boards of a scan. The keywords are:
# scan: scannum srcname
# brd: brdnum numlag numpols
# com: other useful comments
The idea is that you could use a text processor (awk, perl) to scan
the file and find out where the scans, records, and boards start.
The only header information is that in the keyword lines:
scannumber, srcname, number of boards in the scan,number of lags
in the board, and number of polarizations in the board.
An example file is listed below. The lines with --> are not part of the
file format (I added them to put a few comments in this documentation).
--> start of file.. list the keys so people remember the format
# com: keys: scan: scanNum nrecs srcName
# com: keys: rec: recnumber (1..nrecs)
# com: keys: brd: brdNumber(1-4) numlags numpol
# com: data: freq vel pol1 pol2
--> scanNum numRecs srcName
# scan: 224801156 1 MAPS_P463_498572
# rec: 1
--> brd nlags 2pol
# brd: 1 1024 2
--> Freqchn1 velChn1 pol1Data pol2Data
1.3300587E+03 2.0371273E+04 4.3640506E-01 2.3575836E-01
--> channel 2 data
1.3300831E+03 2.0365396E+04 3.8925239E-01 2.0697139E-01
--> .... till end of this board
# brd: 2 1024 2
1.3500587E+03 1.5628420E+04 3.6222970E-01 1.3802753E-01
1.3500831E+03 1.5622716E+04 3.0191767E-01 1.6439275E-01
--> ... for the rest of the boards in this rec
--> ... If there are multiple records in the scan then rec 2 starts..
--> ... else skip to the next scan
# rec: 2
# brd: 1 1024 2
--> ... When all the recs of this scan are done, start the next scan
# scan: 224801160 1 MAPS_P463_790971
# rec: 1
# brd: 1 1024 2
1.3300579E+03 2.0371298E+04 -2.0002886E-03 4.1539539E-04
1.3300823E+03 2.0365421E+04 -1.2376260E-03 1.9774768E-03
EXAMPLE:
process onoff position switch data, convert to janskies then output
ascii file holding just this processed data.
scan=212300123
corposonoff(lun,b,t,cals,scan=scan,/han,/scljy,/sum)
nrecs=corwriteascii(b,'outfile.ascii')
NOTES:
A scan with 1 record of 4 boards, each board having 2 polarizations of
1024 lags each took .25 megabytes. The binary data was 32Kbytes so the file
size gets blownup by about a factor of 8.
(See /pkg/rsi/local/libao/phil/Cor2/corwriteascii.pro)
NAME:
mm0proc - do mueller 0 processing of data in a file.
--> now use spider/pro/mm0/mmproc
SYNTAX: npat=mm0proc(filename,mm,skycor=skycor,astro=astro,$
noplot=noplot,noprint=noprint,keywait=keywait,cumcorr=cumcorr,$
board=board,rcvnum=rcvnum,slinp=slinp,no_mcorr=no_mcorr)
ARGS:
filename: string. name of file with correlator data.
RETURNS:
mm[npat]: {mueller} array of structs with data (1 per pattern)
npat : long number of patterns returned
KEYWORDS:
skycor : int default is to do the sky correction. skycor=0
will not do the sky correction.
astro : int default is to do the astronomical correction.
astro=0 will not do it.
noplot : int if set then do not plot out the spectra. default
is to plot.
noprint : int if set then do not plot out the spectra. default
noplot : int if set then do not print out the results as you
to along. default is to print them out.
keywait : int if set then wait for a keystroke after every pattern.
cumcorr : int if set then do the cum correction to excise rfi
board : int 0 to 3. If set then just process this correlator
board. The default is to process 4 boards.
tcalx[4]: float polA cal values to use for the 4 boards.
If not supplied then look them up.
tcaly[4]: float polB cal values to use for the 4 boards.
If not supplied then look them up.
rcvnum : int rcvr number to process.Default is process all
receivers.
slinp[] : {sl} scan list to use if user has already scaned file
with getsl. If slinp=null and rcvnum != null then
the scanned sl will be returned in slinp.
no_mcorr: if set then no mueller correction is done
(independent of skycor,astro,etc).
DESCRIPTION:
mm0proc will do the mueller0 processing of all of the calibration patterns
in a file. It will return the results (mm) as an array of {mueller} structures
(see mmrestore for a description of whats in the structure).
If the user specifies the board keyword (0 to 3) then only that
board will be processed. By default all 4 boards are processed (a feature/bug
is that integrations using less then 4 boards still loop through all
the board numbers.. the data is ok, it just takes a little longer).
When you are done with this routine you can save the data as
a save file in idl..
save,mm,filenamae=' xxxx'
When starting idl you need to source carl's startup routines:
@~heiles/allcal/idlprocs/xxx/start_cross3.idl
or make your IDL_STARTUP environment variable point there.
SEE ALSO:
mmrestore, mm0proclist
(See /pkg/rsi/local/libao/phil/Cor2/mm0proc.pro)
NAME:
mm0proclist - do mueller 0 processing for a set of data files
--> now use ../spider/pro/mm0/mmproclist
SYNTAX: npat=mm0proclist(listfile,mm,maxpat=maxpat,skycor=skycor,astro=astro,$
noplot=noplot,noprint=noprint,keywait=keywait,$
cumcorr=cumcorr,board=board,rcvnum=rcvnum,fnameind=fnameind)
ARGS:
listfile: string. filename holding names of correlator files
(one per line). semi-colon as the first char is a comment.
RETURNS:
mm[npat]: {mueller} array of structs with data (1 per pattern)
npat : long number of patterns returned
KEYWORDS:
maxpat : long maximum number of patterns to process. The default
is 1000. It will pre-allocate an array of structures
of this size to hold the data.
fnameind: int index of word in line for filename (count from 0)
default is the first word of line. words are separated by
whitespace.
extra : for other keywords see mm0proc
DESCRIPTION:
Call mm0proc for every filename in listfile. Return all of the
data in one array of structures. An example of the listfile is:
; a comment
/share/olcor/calfile.19mar01.a1400.1
/share/olcor/calfile.20apr01.a1446.1
/share/olcor/calfile.20apr01.a1489.1
/share/olcor/calfile.20mar01.a1389.1
It would process all 4 files and return the data in mm.
(See /pkg/rsi/local/libao/phil/Cor2/mm0proclist.pro)
NAME:
mmbeammap - generate beam map from mm data structure.
SYNTAX: mmbeammap,mm,beammap,show=show,mappixels=mappixels,$
mbmap=mbmap,slmap=slmap,axisrange,$
pntsperstrip=pntsperstrip, hpbw=hpbw,nterms=nterms
ARGS:
mm[0]: struct returned from mmgetarchive() holds calibration
fit info.
KEYWORDS:
show : int if true then display images of:
window 0: sidelobe map
window 1: main beam + sidelobe map
the map will be "mappixels square"
the maps are linear scales
mappixels:int the size of the image maps to return. The default is 200
(so bmmap is 200x200 pixels. The fit is evaluated on the
pntsperstripxpntsperstrip grid and this it is interpolated
onto mappixels x mappixels using the idl routine congrid().
Keywords not normally used:
pntsperstrip: int number of points per strip in original data.
The default is 60. This and hpbw deterimine the az,za grid
that is used to evaluate the fit on. Probably best to just
use the default of 60 (since that's what the datataking uses).
hpbw: float half power beamm width used (in arcminutes).
By default it uses what is in the mm structure fit info.
nterms:int number of coef's in the original fft fit to use when creating
the sidelobes. Fit has 8 terms and is hermetian. The default
value is 6. Whatever you enter it must be even and <=8.
RETURNS:
beammap[n,n]:float image of main beam and 1st sidelobes.
n= mappixels which defaults to 200
It is a linear scale with the peak value normalized to 1.
mbmap[n,n] :float image of main beam by itself
slmap[n,n] :float image of first sidelobe by itself.
axisrange[n]:float map angular offsets in arcminutes
(eg beammap[*,0] is evaluated at axisrange[*])
DESCRIPTION:
The spider scan calibration of carl heiles does a 2d fit to the main
beam and the first sidelobe for each spider scan taken (two crosses rotated
by 45 degrees). This information is archived at AO and is accesssible
using the mmgetarchive() routine (in Cor2 lib).
This routine will take the info from 1 pattern (mm) and compute an
image of the main beam and first sidelobe using the original fit results.
The fits are done in great circle az,za. beammap[az,za]. The axisrange[n]
keyword returns the angular offsets where these points were evaluated.
The show=1 keyword will have the routine display the images before
returning.
The maxpixels=N lets the user change the size of the output image.
The fit is always evaluated on a 60x60 grid (the original size of the input
data set). The routine then interpolates this onto the maxpixels x maxpixels
grid. This defaults to 200.
There are a few keywords that allow the user to do non standard
processing: pntsperstrip,hpbw,nterms. You should normally let the routine
use the default values.
EXAMPLE:
- get cband data for 2010, sort by frequency, then build a beammap
for one of the patterns.
idl
@corinit
addpath,'spider/pro/mm0'
.compile mmbeammap
yymmdd1=100101
yymmdd2=101231
istat=mmgetarchive(yymmdd,yymmdd2,mm,rcvnum=9)
; get the unique frequencies data measured at
ufrq=mm[uniq(mm.cfr,sort(mm.cfr))].cfr
print,ufrq
4500.00 4860.00 5000.00 5400.00 5690.00 5900.0
; get 4500 Mhz data
ii=where(mm.cfr eq ufrq[0],cnt)
; make beammap first of these
mmbeammap,mm[ii[0]],beammap,/show,mbmap=mbmap,slmap=slmap,$
axisrange=axisrange
; display the main beam map
imgdisp,mbmap,zx=2,zy=2,nsigclip=3
NOTE:
- This routine uses the ffteval() routine in spider/pro/mm0. Before calling
this routine you need to :
addpath,'spider/pro/mm0'
so the routine can find it. You may have to then say
.compile ffteval if idl doesn't autocompile it.
- You will have to run this routine at AO since the database that
mmgetarchive() accesses is only at ao.
(See /pkg/rsi/local/libao/phil/Cor2/mmbeammap.pro)
NAME:
mmcmpmatrix - compute the mueller matrix from the parameters
SYNTAX:
mm=mmcmpmatrix(mmparams)
ARGS:
mmparams:{mmparams} mm params structure with the parameter
data already loaded: alpha,deltag,epsilon,
phi,psi,chi
RETURNS:
mm[4,4]: float. The computed mueller matrix without the astronomical
correction included.
DESCRIPTION:
Compute the mueller matrix from the parameters that define it. The
parameter values should have already been loaded into the mmparams
structure before calling this routine (see mmgetparams). The 4 by 4 matrix
is returned in the function call. It is not loaded into the mmmparams.mm
element in the structure.
correction to get to sky coordinates.
SEE ALSO:
AO technical memo 2000-05 (The Mueller matrix parameters for Arecibo'S
receiver systems. http://www.naic.edu/aomenu.htm
(See /pkg/rsi/local/libao/phil/Cor2/mmcmpmatrix.pro)
NAME:
mmcmpsefd - compute sefd for each entry in mm structure
SYNTAX: sefd=mmcmpsefd(mm)
ARGS:
mm[n]: {mueller} mueller structure from mmgetarchive
RETURNS:
sefd[n]: float the Sefd System equivalent flux density
DESCRIPTION:
Compute the SEFD for each entry in the mm structure.
(See /pkg/rsi/local/libao/phil/Cor2/mmcmpsefd.pro)
NAME:
mmdoit - input allcal sav files, convert to struct and sav
SYNTAX:
mm=mmdoit(dir=dir)
ARGS: none
KEYWORDS:
dir : string. path for mmfiles.dat input and .sav outputs
The default is '/proj/x102/cal/'
RETURNS:
mm[npats]: {mueller} return the array of mueller structs
one for each pattern used.
DESCRPTION:
The x102 calibration normally runs mm_muelller0.idl to perform the
calibration. It writes the output to a number of save files. These
save files contain the information in arrays. For easier access to
this data, mmdoit will convert these arrays to an array of structures.
The typical sequence is:
1. run mm_mueller0.idl in a particular directory
eg. /share/megs/phil/x101/x102/allcal/test
2. create a directory where you want to store the .sav files as
structures. eg. /share/megs/phil/x101/010104
3. cp /proj/x102/cal/mmgenfiles.sc to the directory in 2.
4. edit mmgenfiles.sc
- set dirlist=( ..) list the directories that have outputs from
mm_mueller0.idl (you can do more than 1 directory at a time).
- set rcvlist=( "12" "9" ) list of receivers you want to process.
5. mmgenfiles.sc in the directory in 2 (in the shell). This will create
the files mmfiles.dat
6. in idl
@corinit
@mminit
mm=mmdoit(dir='/share/megs/phil/x101/010104/')
don't forget the trailling /
7. you can now use mmrestore(dir=dir) to read them back in.
also look at mmrestore documentation for a description of the data
structure.
NOTES:
You need write access to the directory in step 2 above.
(See /pkg/rsi/local/libao/phil/Cor2/mmdoit.pro)
NAME:
mmfindpattern - get the indices for the start of the patterns
SYNTAX: nfound=mmfindpattern(sl,indar,stripspat=stripspat)
ARGS:
sl[]: {getsl} scan list array from getsl
KEYWORDS:
stripspat: int minimum number of strips in pattern to count
at ok default is 4.
RETURNS:
ind[npat]: long indices into sl[] for cal rec at start of each pattern
npat : long number of patterns found
DESCRIPTION:
The getsl() routine will routine a scanlist array holding info about
every scan in a file. This routine will then search the sl[] array and
find the indices where the heiles scans start. For it to be included
there must be at least stripspat strips in the pattern (the default is 4).
EXAMPLE:
openr,lun,'/share/olcor/corfile.23aug02.x101.1',/get_lun
sl=getsl(lun)
npat=mmfindpattern(sl,indfound)
for i=0,npat-1 do begin
scan=sl[indfound[i]].scan
.. process this scan
(See /pkg/rsi/local/libao/phil/Cor2/mmfindpattern.pro)
NAME:
mmget - extract a subset of mueller array by key
SYNTAX: mnew=mmget(mm,count,freq=freq,rcvnum=rcvnum,rcvnam=rcvnam,brd=brd,$
src=src)
ARGS:
mm[] :{mueller} array of mueller structs read in via mmrestore
(or a subset of any of these)
KEYWORDS:
freq : float .. freq in Mhz to extract
rcvnum : long .. receiver number to extract:
1=327,2=430,3=610,5=lbw,6=lbn,7=sbw,9=cb,12=sbn
100=430ch
brd : long correllator board number to return 0-3
src : string source name to return
rcvnam : string rcvname to get:327,430,610,lbw,lbn,sbn,sbc,cband,430ch
RETURNS:
mmnew[count]: {mueller} subset meeting the requested critiera
count : long number of patterns found
DESCRIPTION:
The x102 calibration data can be input via mmrestore. It will input
arrays of {mueller} structures containing reduces data for all the receivers
and an individual array for each receiver:
eg: mm[] holds all the data
mm12[] holds only the rcvnumber 12: sbn data.
You can display the structure format with:
help,mm,/st and
help,mm.fit,/st total power fit info
help,mm.fitq,/st q fit info
help,mm.fitu,/st u fit info
help,mm.fitv,/st v fit info
mmget allows you to select a subset of one of these arrays. It uses
the specified keyword to determine what to return. It uses the
first keyword found.
EXAMPLES
return all of the brd 0 data from lbw
aa5=mmget(mm5,count,brd=0)
return the source B0106+130 for cband data
aa9=mmget(mm9,count,src='B0106+130')
return all the board 1 data for the above cband source.
aa9b0=mmget(aa9,count,brd=0)
NOTES:
beware that you donot accidentally wipe out the input array
by using that name as the output..
eg: mm5=mmget(mm5,brd=0) changes mm5 to only have the brd 0 data.
(See /pkg/rsi/local/libao/phil/Cor2/mmget.pro)
NAME:
mmgetarchive - restore all or part of calibration archive
SYNTAX: count=mmgetarchive(yymmdd1,yymmdd2,mm,rcvnum=rcvnum,dir=dir)
ARGS:
yymmdd1 : long year,month,day of first day to get (ast)
yymmdd2 : long year,month,day of last day to get (ast)
KEYWORDS:
rcvnum : long .. receiver number to extract:
1=327,2=430,3=610,5=lbw,6=lbn,7=sbw,8=sbw,9=cb,$
11=xb,12=sbn,100=430ch
dir : string If provided, then this is the path to acces the
archive save files. The default is:
/share/megs/phil/x101/x102/runs/ .
RETURNS:
mm[count]: {mueller} data found
count : long number of patterns found
DESCRIPTION:
This routine will restore the x102 calibration data stored in
/share/megs/phil/x101/x102/runs. It is updated monthly. You specify
the start and end dates of the data to extract. You can optionally
specify a receiver with keyword rcvnum. Once the data has been input
you can created subsets of the data with the mmget() routine.
See mmrestore for a description of the structure format.
EXAMPLES
;get all data for jan02->apr02
nrecs=mmgetarchive(020101,020430,mm)
;get the cband data for apr02
nrecs=mmgetarchive(020101,020430,mm,rcvnum=9)
;from the cband data extract the 5000 Mhz data
mm5=mmget(mm,count,freq=5000.)
(See /pkg/rsi/local/libao/phil/Cor2/mmgetarchive.pro)
NAME:
mmgetparams - input mueller matrix for rcvr.
SYNTAX:
istat=mmgetparams(rcvNum,cfr,mmparams,fname=fname,date=date)
ARGS:
rcvNum: 1 thru 16. receiver to use (same as hdr.iflo.stat1.rfnum).
cfr: float freq in Mhz where the values should be computed.
KEYWORDS:
fname: to specify an alternate data file with mueller matrix values.
The default file is aodefdir() + 'data/mm.datR{rcvNum}
date : [year,daynum] .. if specified the data when you want the
mueller matrix data for..default is most recent.
RETURNS:
istat: 1 ok
-1 couldn't open file (probably no mueller data.
mmparams:{mmparam} return mueller matrix info.
DESCRIPTION:
Input the mueller matrix parameters and for the receiver number rcvnum.
By default the data is read from the file:
aodefdir() + 'data/mm.datR{rcvNum}
aodefdir() is a function that returns the root of the aoroutines.
The keyword fname lets you use an alternate file. The date keyword will
search for data valid on or after this date (the default is to take the most
recent data).
There is 1 file per receiver that holds the mueller matrix info for
this receiver. Data blocks are defined by the date that they became
valid. The datafile format is:
- col 1 # : this is a comment line. They can be anywhere in the file
- !yyyy daynum : this is the start of the data set for data valid on or
after yyyy daynum. !0 1 will work for any time in the past.
- end in col 1: This marks the end of the dataset that started with the
previous !yyyy daynum
- The lines between the !yyyy daynum and the "end" line are idl
executable statements that define the parameters.
The input routine will scan the file looking for the date record that
the person wants to use. It will then input 1 line at a time and execute
them as idl statements (lines starting with # are ignored) until the
line with "end" in cols 1-3 is found. These lines must define the following
variables: alpha,epsilon,phi,chi,deltag,angle_astron,m_astron[4,4],
circular, and corcal.
An example file would look like:
# following block valid starting on daynumber 145 of 2002 (25may02).
#
!2002 145
alpha = .25*!dtor
epsilon= .0015
phi = -148.*!dtor
psi = -175.4*!dtor
chi = 90.*!dtor
cfr20=cfr - 1420.
deltag= 0.100 + 0.015* cos( 2.* !pi* cfr20/300.)
angle_astron=-45.
angle=angle_astron*!dtor
m_astron=fltarr(4,4)
m_astron[0,0] =1.
m_astron[3,3] =-1.
m_astron[ 1,1]= cos( 2.* angle)
m_astron[ 2,1]= sin( 2.* angle)
m_astron[ 2,2]= m_astron[ 1,1]
m_astron[ 1,2]= -m_astron[ 2,1]
corcal=1
circular=0
end
# following block valid after daynumber 130 of 2001
!2001 130
alpha = .25*!dtor
epsilon= .0015
phi = -148.*!dtor
psi = -175.4*!dtor
chi = 90.*!dtor
..
..
end
# valid for any time before daynum 130 of 2001
!0 1
alpha = .25*!dtor
epsilon= .0015
..
..
end
These variables are then loaded into the mmparams (mueller matrix
parameter structure. The structure format for {mmparams} is:
mmp.rcvNum int receiver number
mmp.year int year when this data becaume valid
mmp.day int daynumber of year when this data became valid
mmp.cirular int 1 if receiver is native circular
mmp.corcal int 1 if receiver has correlated cal
mmp.cfr float frequency in Mhz where the parameters are computed.
mmp.alpha float alpha parameter (in radians)
mmp.epsilon float epsilon parameter
mmp.phi float phi angle in radians
mmp.psi float psi angle in radians
mmp.chi float chi angle in radians
mmp.deltag float difference in cal values
mmp.astronAngle float astronomical angle (degrees)
mmp.m_astron[4,4] float matrix to apply after the mueller matrix
correction to get to sky coordinates.
SEE ALSO:
AO technical memo 2000-05 (The Mueller matrix parameters for Arecibo's
receiver systems. http://www.naic.edu/aomenu.htm
(See /pkg/rsi/local/libao/phil/Cor2/mmgetparams.pro)
NAME:
mmplot - x,y plot using different colors for each receiver
SYNTAX: mmplot,x,y,mm,xtitle=xtitle,ytitle=ytitle,title=title,over=over
ARGS: x[npts] : float xarray to plot
y[npts] : float yarray to plot
mm[npts] : {mueller} mueller info
KEYWORDS:
xtitle : string label for x axis
ytitle : string label for y axis
title : string label top of plot
sym : if set then use symbols rather than color
DESCRIPTION:
mmplot will plot the x,y data with a different color
for each reciever.
EXAMPLE:
mmplot,mm.za,mm.fit.tsys,mm,xtitle='za',ytitle='tsys[K]',title='tsys vs za'
(See /pkg/rsi/local/libao/phil/Cor2/mmplot.pro)
NAME:
mmplotazza - plot az,za coverage of dataset.
SYNTAX: mmplotazza,mm,sym=sym,dx=dx,over=over,color=color,title=title,fig=fig
ARGS:
mm[npts]: {mueller} mueller info
KEYWORDS:
sym: int symbol to plot at each position.Default is *.
over: if set then overplot this data with what is there.
dx: The step size in feet along the x,y axis. default is 10 feet.
color: int color index to use : 1..10 (see ldcolph)
title: string label for title of plot
fig : int figure number
DESCRIPTION:
Plot the azimuth, za positions as a cartesian x,y plot. The axes are
feet from the center of the dish (projected onto z=0). This routine can
give an idea of how well a set of sources has covered the dish.
EXAMPLE:
nrecs=mmgetarchive(020826,021015,mm,rcv=9)
mmplotazza,mm,title='cband coverage after focus change till 15oct02'
(See /pkg/rsi/local/libao/phil/Cor2/mmplotazza.pro)
NAME: mmplotbydate - plot calibration observations by receiver and date SYNTAX: mmplotbydate,date=date KEYWORDS: date[2]: long start, end date for plot. format is [yymmdd1,yymmdd2] flagmon: if set then flag the start of each month. DESCRIPTION: Make a plot of the cumulative calibration observations done vs date. Each receiver is plotted separately. By default the observations done in the previous 12 months are used. The keyword date lets you specify a different date range. NOTE: This routine uses the calibration database which is updated at the end of each month. The current months data will not be available until the start of the next month.
(See /pkg/rsi/local/libao/phil/Cor2/mmplotbydate.pro)
NAME:
mmplotcsme - plot coma,SideLobeHght, and beam efficiencies.
SYNTAX: mmplotcsme,mm,col=col,vc=vc,vs=vs,vm=vm,ve=ve,tit=tit,lnf,$
lnf=lnf,xpfrq=xpfrq,fig=fig,fgln=fgln,xpfig=xpfig,lns=lns,xps=xps,$
srccol=srccol,msrc=msrc,bypix=bypix,cs=cs,font=font
ARGS:
mm[n]: {mueller} data structure to plot
KEYWORDS:
col[m]: int color array to use
vc[2] : float vertical range for coma [min,max]
vs[2] : float vertical range for sideLobeHght [min,max]
vm[2] : float vertical range for main beam efficiency [min,max]
ve[2] : float vertical range for main beam +1st sidelobe efficiency
lnf : float line number to print frquencies
xpfrq : float 0.,1. horizontal position for frequency
fig : int figure number
fgln : float line number for figure number
xpfig : float 0,1. horizontal position for figure number
lns : float line number to start source labels
xps : float 0,1. horizontal position for source labels
srccol: if set then use colors for sources rather than freq
msrc : int max number of source names to print on left side of page
before moving to the right. (def=30)
bypix : if set then colors are for pixels not frequency
cs : float chars size for plots. def=1.5
font : int font to use.
DESCRIPTION:
Plot the coma parameter, 1st sidelobe height (in db), the main beam
efficiency , and the main beam + 1st sidelobe beam efficiency for all
of the sources in the mm array. Plot each source with a separate symbol.
Plot each frequency as a separate color. It is probably a good idea to
plot 1 receivers worth of data at a time.
You can change the colors used via the col keyword (use number 1 through
10 for colors..1=black,2-red,3-green,4-blue,5..). The vX array allows
you to set the vertical scale for each plot. The default is to autoscale
to the max, min of the data. The lnX keywords let you position the
source names on the plot.
The mm array is normally generated from the mm0proc routine.
EXAMPLES:
plot the xband info for jan02.
restore,'/share/megs/phil/x101/x102/runs/c020101_020131.sav'
ind=where(mm.rcvname eq 'xb') ; just get the xband data
mm=mm[ind]
fig=1
tit='02jan02 Xband'
mmplotcsme,mm,vc=[0,.3],vs=[-15,0],vm=[0,.5],ve=[0,.5],fig=fig,tit=tit
SEE ALSO:
mmplotgtsb, mm0proc.
(See /pkg/rsi/local/libao/phil/Cor2/mmplotcsme.pro)
NAME:
mmplotgtsb - plot gain, Tsys, sefd, and beamWidth for sources
SYNTAX: mmplotgtsb,mm,col=col,vg=vg,vt=vt,vs=vs,vb=vb,tit=tit,lns=lns,$
lnf=lnf,xpfrq=xpfrq,fig=fig,fgln=fgln,xpfig=xpfig,xps=xps,$
srccol=srccol,bwamin=bwamin,msrc=msrc,bypix=bypix
charth=charth,thick=thick,cs=cs ,font=font
ARGS:
mm[n]: {mueller} data structure to plot
KEYWORDS:
col[m]: int color array to use
vg[2] : float vertical range for gain [K/Jy] ,[min,max]
vt[2] : float vertical range for Tsys [K] ,[min,max]
vs[2] : float vertical range for sefd [Jy/Tsys] ,[min,max]
vb[2] : float vertical range for beamWidth[asecs],[min,max]
tit : string title for top of plot.
lns : float line number to start plotting source names.def=2
lnf : float line number to print frquencies/pixels. def=2
xpfrq : float 0.,1. hor position for frequency/pixel. def=.25
fig : int figure number
fgln : float line number for figure number
xpfig : float 0,1. horizontal position for figure number
xps : float 0,1. horizontal position for source labels
srccol: if set then use colors for sources rather than freq/pix
bwamin: if set then plot beam width in arc min rather than asecs
msrc : int max number of source names to print on left side of page
before moving to the right. (def=30)
bypix : if set then colors are for pixels not frequency
byaz : if set then plot vs azimuth. default is za
DESCRIPTION:
Plot the gain, Tsys, Sefd, and average beam width for all
of the sources in the mm array. Plot each source with a separate symbol.
Plot each frequency as a separate color. It is probably a good idea to
plot 1 receivers worth of data at a time.
You can change the colors used via the col keyword (use number 1 through
10 for colors..1=black,2-red,3-green,4-blue,5..). The vX array allows
you to set the vertical scale for each plot. The default is to autoscale
to the max, min of the data. The lns keyword lets you position the
source names on the plot.
The mm array is normally generated from the mm0proc routine.
Before calling this routine you should execute @corinit.
EXAMPLES:
plot the xband info for jan02.
restore,'/share/megs/phil/x101/x102/runs/c020101_020131.sav'
ind=where(mm.rcvname eq 'xb') ; just get the xband data
mm=mm[ind]
fig=1
tit='02jan02 Xband'
mmplotgtsb,mm,vg=[0,5.],vt=[40,60],vs=[0,30],vb=[30,45],fig=fig,tit=tit
SEE ALSO:
mmplotcsme, mmplotpnterr,mm0proc.
(See /pkg/rsi/local/libao/phil/Cor2/mmplotgtsb.pro)
NAME:
mmplotpnterr - plot pointing error.
SYNTAX: mmplotpnterr,mm,col=col,vza=vza,vaz=vaz,vtot=vtot,hza=hza,haz=haz,$
htot=htot,lns=lns,xps=xps,lnm=lnm,xpm=xpm,fig=fig,fgln=fgln,$
xpfig=xpfig,rise=rise,tit=tit ,notot=notot,msrc=msrc,bypix=bypix,$
byfrq=byfrq,cs=cs,fln=fln,font=font,fxp=fxp
ARGS:
mm[n]: {mueller} data structure to plot
KEYWORDS:
col[m]: int color array to use
vza[2] : float vertical range for zaErr in asecs ,[min,max]
vaz[2] : float vertical range for azErr in asecs ,[min,max]
vtot[2] : float vertical range for total pntErr in asecs ,[min,max]
hza[2] : float horizontal range for za in deg [min,max]
haz[2] : float horizontal range for az in deg [min,max]
lns : float line number to start source labels
xps : float 0.,1. horizontal position for source labels
lnm : float line number to start mean,rms labels
xpm : float 0.,1. horizontal position for mean,rms labels
fig : int figure number
fgln : float line number for figure number
xpfig : float 0,1. horizontal position for figure number
rise : if set then use solid, dash lines to differentiate between
rise set. default is to just plot symbols with no connecting
lines
tit : string title for top of page.
notot : if set then don't bother plotting the total errors.
msrc : int max number of source names to print on left side of page
before moving to the right. (def=30)
bypix : if set then colors are for pixels not frequency
byfrq : use colors to separate frequency (default is to use
colors to separate sources.
cs=cs : float if supplied then character size
fln : 0,31 line to start writing the frequencies
fxp : 0..1 xpos to start writing the freq. def=.25
DESCRIPTION:
Plot the za pointing error, az pointing error, and the total pointing
error (sqrt(zaErr^2+azErr^2)) versus zenith angle and azimuth.
of the sources in the mm array. Plot each source with a separate symbol
and color. The pointing error at each frequency will be overplotted.
You can change the colors used via the col keyword (use number 1 through
10 for colors..1=black,2-red,3-green,4-blue,5..). The vX array allows
you to set the vertical scale for each plot. The default is to autoscale
to the max, min of the data. The lnX keywords let you position the
source names on the plot.
The mm array is normally generated from the mm0proc routine.
Before calling this routine you should execute @corinit.
EXAMPLES:
plot the xband info for jan02.
restore,'/share/megs/phil/x101/x102/runs/c020101_020131.sav'
ind=where(mm.rcvname eq 'xb') ; just get the xband data
mm=mm[ind]
fig=1
tit='02jan02 Xband'
mmplotpnterr,mm,vg=[0,5.],vt=[40,60],vs=[0,30],vb=[30,45],fig=fig,tit=tit
SEE ALSO:
mmplotcsme, mmplotgtsb,mm0proc.
(See /pkg/rsi/local/libao/phil/Cor2/mmplotpnterr.pro)
NAME:
mmplotsrc - x,y plot using different colors for each source
SYNTAX: mmplotsrc,x,y,mm,xtitle=xtitle,ytitle=ytitle,title=title,over=over
xp=xp,ln=ln,nolab=nolab,col=col,sym=sym,sclln=sclln,
decord=decord,rise=rise,fln=fln,_extra=e,srcl=srcl,msrc=msrc
xlp=xlp,thick=thick,charth=charth,font=font,connect=connect
ARGS: x[npts] : float xarray to plot
y[npts] : float yarray to plot
mm[npts] : {mueller} mueller info
KEYWORDS:
xtitle : string label for x axis
ytitle : string label for y axis
title : string label top of plot
over : if set then overplot the data with what is there.
xp : float xposition for srcnames 0..1
ln : int line to start source list 0..33
sclln : float fractional step between labels (default 1.)
rise : if set then plot,rise,set as separate linestyles
sym : if set then use symbols rather than color to distinguish
sources. use 1st elm of col array for all colors
col[nsrc]:use the values for the color (1..10) rather than the defaults
nolab : if set then donot print the source names
decord : if set then plot src in dec order, rather then ra.
useful if x is just lindgen(npts) and you want to plot
the sources by increasing ra or dec.
byfrq : Use symbols for sources, color for freq. (does not work
with /rise option
bypix : if set then use symbols for sources color for pixels.
only for alfa, does not work with /rise option.
fln : float line to start writing the frequencies
fxp : float xpos (0...1) to start writing the frequencies
srcl[] : string user passes in source list rather than searching
through mm. Use this if you are making multiple
plots per page and all the plots may not have all the
same sources (eg. lbw 1415 Mhz is shared by the high
and low lbw measurements).
msrc : int max sources to list on left before move to right col
xlp : float xpos (0..1) for 2nd col. def:1.
connect : int if true then plot with solid lines connecting points of a source.
Useful if plotting vs az.
_extra : e .. passed to plot routine.
DESCRIPTION:
mmplotsrc will plot the x,y data with a different color
for each source.
EXAMPLE:
mmplotsrc,mm.za,mm.fit.tsys,mm,xtitle='za',ytitle='tsys[K]',
title='tsys vs za'
(See /pkg/rsi/local/libao/phil/Cor2/mmplotsrc.pro)
NAME:
mmrestore - input the muellar structure arrays.
SYNTAX: @mmrestore
ARGS: none
RETURNS:
mm[] : {mueller} all the data in a single array. sorted by receiver and
then source
mm1[] : {mueller} the 327 data
mm2[] : {mueller} the 430 data
mm3[] : {mueller} the 610 data
mm5[] : {mueller} the lbw data
mm6[] : {mueller} the lbn data
mm7[] : {mueller} the sbw data
mm9[] : {mueller} the cband data
mm12[] : {mueller} the sbn data
mm100[] : {mueller} the 430ch data
DESCRIPTION:
mmrestore will input the mueller structure arrays for all of the
receivers that have been reduced in the sep00 calibration run. All it
does is :
restore,'/proj/x102/cal/mmdata.sav',/verbose
If you have run mmdoit() and created a save file of a different
set of data, you can input that data with:
restore,'filename',/verbose
(you do not need to use this routine).
The structure format returned is described below. Those variables with
(I) in the comments fiels refer to total power I and not a single
polarization (eg Tsys ...... (I)).
IDL> help,mm,/st
** Structure MUELLER, 20 tags, length=372:
SRCNAME STRING 'B0300+162'
SRCFLUX FLOAT -1.00000 Jy (1 pol)
SCAN LONG 26500134
RA1950 FLOAT 3.00748 hours
DEC1950 FLOAT 16.2511 degrees
RCVNUM LONG 3
RCVNAM STRING '610' .. receiver name. valid names are:
327,430,610,lbw,lbn,sbn,sbw,cband,
430ch
UTSEC LONG 6 ut second from midnight
JULDAY FLOAT 51808.8 (MJD + .5) or (JD-2400000)
BANDWD FLOAT 1.56250 Mhz
CFR FLOAT 612.000 Mhz
BRD INT 0 corr brd #. 0-3
CALTEMP FLOAT Array[2] deg K
LST FLOAT 1.88105 hours
AZ FLOAT -86.4037 degrees
ZA FLOAT 16.8460 degrees
PARANGLE FLOAT -56.4730 degrees (paralactic angle)
ASTRONANGLE FLOAT 0.0000 degrees feed to astromical system
PASRC FLOAT 74.1880 deg position angle source on sky
PASRC_ERR FLOAT .2780 deg ..error
POLSRC FLOAT .0391 source fractional linear pol
POLSRC_ERR FLOAT .0004 err
BMWIDSCAN FLOAT 8.00000 Amin
MMCOR INT 2 muelcorrection.0=no,1-toaz/za,2=sky
FIT STRUCT -> MUELLERFITI
FITQ STRUCT -> MUELLERFITPOL
FITU STRUCT -> MUELLERFITPOL
FITV STRUCT -> MUELLERFITPOL
MMPARM STRUCT -> MUELLERPARAMS params used for mueller matrix
the total power fit info is:
IDL> help,mm.fit,/st
** Structure MUELLERFITI, 29 tags, length=132:
TSYS FLOAT 268.714 degK (I)
TSYS_ERR FLOAT 0.121488 degK (I)
GAIN FLOAT -41.4766 K/Jy (1 pol)
DTSYSDZA FLOAT 0.673788 K/Degza
DTSYSDZA_ERR FLOAT 0.248539 K/Degza
TSRC FLOAT 82.9533 K (I)
TSRC_ERR FLOAT 0.358256 K (I)
SIGMAPNTS FLOAT 1.47459 K
AZERR FLOAT -0.507198 Amin
AZERR_ERR FLOAT 0.0310421 Amin
ZAERR FLOAT 0.0396541 Amin
ZAERR_ERR FLOAT 0.0375799 Amin
BMWIDAVG FLOAT 7.65640 Amin (max+min)/2
BMWIDAVG_ERR FLOAT 0.0395753 Amin
BMWIDDELTA FLOAT 0.556200 Amin (max-min)/2
BMWIDDELTA_ERR FLOAT 0.0418182 Amin
BMPHI FLOAT 92.4938 deg
BMPHI_ERR FLOAT 2.17635 deg
COMA FLOAT 0.00683132 hpbw
COMA_ERR FLOAT 0.0153134 hpbw
COMAPHI FLOAT -44.1549 deg
COMAPHI_ERR FLOAT 129.217 deg
SLHGT FLOAT .02 fraction of MainBeam
SLCOEF COMPLEX Array[8] coef from fft of sdlbHgt/mainBeam
ETAMB FLOAT -2.69011 main beam efficiency
ETASL FLOAT -0.243481 sidelobe efficiency
CALPHASE FLOAT Array[2] a+b*(freq-cfr) Rd/Mhz
CALPHASE_ERR FLOAT Array[2]
SRCPHASE FLOAT Array[2] a+b*(freq-cfr) Rd/Mhz
SRCPHASE_ERR FLOAT Array[2]
The polarized fit info is:
IDL> help,mm.fitq,/st
** Structure MUELLERFITPOL, 14 tags, length=56:
OFFSET FLOAT 0.0201661 zero offset degK
OFFSET_ERR FLOAT 0.0593501
DOFFDZA FLOAT 0.623669 degK/degza
DOFFDZA_ERR FLOAT 0.134090
SRC FLOAT 5.66668 src deflection. fraction of I
SRC_ERR FLOAT 0.179385 fraction of I
SQUINTAMP FLOAT 0.0571257 Amin
SQUINTAMP_ERR FLOAT 0.0115344
SQUINTPA FLOAT 25.4974 position angle deg
SQUINTPA_ERR FLOAT 11.8336
SQUASHAMP FLOAT 0.169447 arcmin hpbw
SQUASHAMP_ERR FLOAT 0.0293442
SQUASHPA FLOAT -82.4178 position angle deg
SQUASHPA_ERR FLOAT 4.51572
there are 3 arrays of polarized data:
mm.fitq
mm.fitu
mm.fitv
mm.parameters used for mueller matrix correction
help,mm.mmparm,/st
** Structure MUELLERPARAMS, 6 tags, length=24:
DELTAG FLOAT -0.0230000
EPSILON FLOAT 0.00600000
ALPHA FLOAT 0.00872665
PHI FLOAT 1.37183
CHI FLOAT 1.57080
PSI FLOAT 2.61799
SEE ALSO
mmget
(See /pkg/rsi/local/libao/phil/Cor2/mmrestore.pro)
NAME:
mmselectgood - select good points from mm array
SYNTAX: n =mmselectgood(mm,mmGood,ltI,tdbad=tdbad,lrTdI=lrTdI,ldate=ldate,hard=hard,$
gain=gain,bmwid=bmwid,vg=vg,vb=vb,retall=retall,verb=verb)
ARGS:
mm[m]:{} calibration data to use as input
KEYWORDS:
tdbad[] : int specify tiedowns that have bad tensions
(count 1..6 -> td12-1,td12-2,td4-1,td4-2,td8-1,td8-1
the routine will copy the tension of the "other" cable
for that td. If both cables on a td are bad, you are
out of luck.
ldate : string for plots
hard : int if true then make hardcopy of result
psname : int if hard then psname to use. def: idl.ps
gain : int if set then include gain plot
bmwid : int if set then include beamwid
vg [2] : float vertical range for gain plot
vb [2] : float vertical range for beamwidth
ret_all : if set then return all, no user input
verb : if set then print progress
RETURNS:
n : long number of good points in mmGood
mmGood[n] :{} good data selected from mm
ltI[n] :{} tdkips and platform height for each good pnt
lrTdI[] :{} laser range and td info input that overlaps the
the mm time range.
DESCRIPTION:
Let the user select calibration scans that are good,
The routine will input the platform height and td tension data that
corresponds to the mm[m] time range.
It will then plot the azErr, zaErr, td tensions, and platform height.
The user can then interactively exclude points that have a bad platform height
or low td tension.
The array mmGood[n] are the resulting good calibration scans.
ltI[n] holds the interplated platform height, tdKips for each mmGood entry.
The lrTdI[] array holds the original td, platform height data input.
Before calling this routine, you must
@lrinit,@tdinit,@corinit
unction mmselectgood,mm,mmG,ltI,tdbad=tdbad,lrtdI=lrTdI,ldate=ldate,hard=hard,psname=psname,$
gain=gain,bmwid=bmwid,vb=vb,vg=vg,ret_all=ret_all,verb=verb
common colph,decomposedph,colph
use3=1 ; if missing distomat, use 3 dist model
verbl=keyword_set(verb)
dohard=keyword_set(hard)
if n_elements(ret_all) eq 0 then ret_all=0
if (dohard) then begin
if n_elements(psname) eq 0 then psname='idl.ps'
endif
if n_elements(vb) ne 2 then begin
vbl=[0,220]
endif else begin
vbl=vb
endelse
if n_elements(vg) ne 2 then begin
vgl=[0,12]
endif else begin
vgl=vg
endelse
minKipsOk=1.
ldate=(n_elements(ldate) gt 0) ? ldate:''
get the lr,td data for this time range
if verbl then print,"start mmgetlrtd"
n=mmgetlrtd(mm,lrTdI,tdbad=tdbad,verb=verb,use3=use3)
if n eq 0 then begin
print,"No lr,td info found for this time range"
return,0
endif
if verbl then print,"done mmgetlrtd"
for every mm entry record min,max,avg value for kips,platHght
a={ min:0.,$
max:0.,$
avg:0.}
aa={ mjd:0d,$ ; start
nsmp:0l,$
kips:replicate(a,3),$ ; use total kips td
plh :a}
pattern length, fraction of a day
patDura=321d/86400d
eps=10/86400.
nmm=n_elements(mm)
ltIinp=replicate(aa,n)
loop through mm finding min,max,avg lr,td info for each entry.
for i=0,nmm-1 do begin &$
ii=where((lrtdi.mjd ge (mm[i].julday-eps)) and (lrtdi.mjd le (mm[i].julday+patDura+eps)),cnt) &$
default to mm start. update ot lr time if we find lr samples
ltIinp[i].mjd=mm[i].julday &$
if cnt gt 0 then begin &$
for j=0,2 do begin &$
ltIinp[i].kips[j].min=min(lrtdI[ii].tdkips[*,j]) &$
ltIinp[i].kips[j].max=max(lrtdI[ii].tdkips[*,j]) &$
ltIinp[i].kips[j].avg=mean(lrtdI[ii].tdkips[*,j]) &$
endfor &$
; for avg plathght, exclude samples with no readings
avghPl=lrtdI[ii].avghpl &$
jj=where(avghPl gt 1.,cnt) &$
if (cnt ne 0) then begin &$
ltIinp[i].plh.min=min(avghPl[jj]) &$
ltIinp[i].plh.max=max(avghPl[jj]) &$
ltIinp[i].plh.avg=mean(avghPl[jj]) &$
endif else begin &$
ltIinp[i].plh.min=0. &$
ltIinp[i].plh.max=0. &$
ltIinp[i].plh.avg=0. &$
endelse &$
ltIinp[i].mjd=lrtdI[ii[0]].mjd &$
endif &$
ltIinp[i].nsmp=cnt &$
endfor
now loop twice
once to let them specify bad points
a second time to plot the results..
plot vs start time of mm array
edo:
selectBadDone=0 ; see if in the select loop
if ret_all then begin
ok=intarr(n)+1
iiOk=where(ok eq 1,cntiiok)
selectBadDone=1
endif
nbad=0l
while 1 do begin
hor
ver
x0=mm[0].julday ; start of x axis
xmm=(mm.julday - x0)*24. ; put in hours from start
xltInp=(ltIinp.mjd - x0)*24.
find points with any td tension less than minkipsOk
ii=where((reform(ltIinp.kips[0].min) lt minkipsok) or $
(reform(ltIinp.kips[1].min) lt minkipsok) or $
(reform(ltIinp.kips[2].min) lt minkipsok),cnt)
az=mm.az
za=mm.za
azErr=mm.fit.azerr*60
zaErr=mm.fit.zaerr*60
gainM=mm.fit.gain
bmWidM=mm.fit.bmwidavg*60.
smp=lindgen(nmm)
symAz=-1
symZa=-2
col=[1,3,4,5,6,7,8,9,10]
cs=2
if dohard and selectBadDone then pscol,psname,/full
ll=4
if keyword_set(gain) then ll++
if keyword_set(bmwid) then ll++
!p.multi=[0,1,ll]
ver,-30,30
hor,-10,nmm+5
xind=lindgen(nmm)
x=xind
sclln=.7
plot vs index
ln=2
ver,-40,40
mmplotsrc,x,zaerr,mm,/rise,col=col,chars=cs,sclln=sclln,ln=ln,$
xtit='pattern #',yit='Asec',tit=ldate + ' Za error vs pattern index'
if cnt gt 0 then oplot,x[ii],zaerr[ii],psym=1,col=colph[2]
flag bad points we actually used, only after selectBadDone
if selectBadDone then begin
for i=0,nbad-1 do begin &$
jj=where(mm.scan eq uscanbad[i]) &$
flag,jj,linestyle=1 &$
endfor &$
endif
(See /pkg/rsi/local/libao/phil/Cor2/mmselectgood.pro)
NAME: mmtostr - move mueller processed arrays to a structure
(See /pkg/rsi/local/libao/phil/Cor2/mmtostr.pro)
NAME:
mm_chkpattern - check if a mueller pattern is complete
SYNTAX: stat=mm_chkpattern(scan,sl)
ARGS:
scan : long scan number for cal on at start of pattern
sl[] : {getsl} scan list array
RETURNS:
istat : long >=0 pointer into sl array where pattern starts.
<0 not a complete pattern
-1 scan not in file
-2 pattern not complete , hit eof
-3 not a complete pattern, wrong rec types
-4 not a complete pattern, wrong rec nums
DESCRIPTION:
Check that a mueller pattern is complete. You must first input a
scanlist array for the file. The program will check for the correct record
types and number of records. If the pattern is complete then the routine
will return the index into the sl array that points to the scan that starts
the pattern. If is not complete then istat will be less than zero.
EXAMPLES:
file='/proj/a1489/calfile.17may01.a1489.1'
openr,lun,file,/get_lun
;
; scan then entire file, then search for scan 113700159L
; if it's ok, read all the records in
;
sl=getsl(lun)
scan=113700159L
istat=mm_chkpattern(scan,sl)
if istat ge 0 then istat=corgetm(lun,242,bb,scan=scan,sl=sl)
free_lun,lun
; scan only the 6 records about the requested scan (this will
; go faster if you only are interested in a single pattern)
; if it's ok, read all the records in
sl=getsl(lun,scan=scan,maxscans=6)
if istat ge 0 then istat=corgetm(lun,242,bb,scan=scan,sl=sl)
SEE ALSO:
mm_findpattern
(See /pkg/rsi/local/libao/phil/Cor2/mm_chkpattern.pro)
NAME:
mm_mon - monitor online datataking calibration scans
SYNTAX: @mm_mon
ARGS:
If the following variables are setn they will be used rather thean
the defaults:
startscan=scannumber if you don't want to start at the beginnning.
usecal=1 By default the gui uses /share/olcor/calfile to write calibrate
scans. When you are doing long tracks on sources the line oriented
code used /share/olcor/corfile for the output file. By default
the code used corfile. To use calfile set usecal=1 before running
the script.
usebrd=n By default board 0 is used. setting this varible will proces
a different board (0..3)
useextfile If set then the user is supplying the external filename in the
variable extfile. I you want to switch back to the online
files, be sure and set this variable to zero.
You need to use carls startup script if you are using this routine:
setenv IDL_STARTUP ~heiles/allcal/idlprocs/xxx/start_cross3.idl
before you start idl.
(See /pkg/rsi/local/libao/phil/Cor2/mm_mon.pro)
NAME:
pfcalib - intensity calibrate an entire file
SYNTAX: istat=pfcalib(filename,bf,bcalf,han=han,
scan=scan,dir=dir,sl=sl,maxrecs=maxrecs,
edgefract=edgefract,mask=mask,
bpc=bpc,smobpc=smobpc,blrem=blrem,svd=svd,$
rcv=rcv,sav=sav,fsav=fsav,rsl=rsl)
ARGS:
filename: string .. file to process
KEYWORDS
han: if set then hanning smooth the data
scan: long if set then start at this scan number rather than the
beginning of the file.
dir[]: string search for filename in these directories.
in this case filename should not have a path.
sl[]: {slist} Scan list to use when searching file. This saves having
to scan the file. Use this to speed processing if you
have already scanned the file (see the Rsl keyword).
maxrecs : long The maximum number of recs to allow in a scan
the default is 300.
edgefract[1/2]: float fraction of bandpass on each side to not use during
calibration. default .1
mask:{cormask} mask structure created for each brd via cormask routine
use this rather than edgefract.
bpc: int 1 band pass correct with cal off
2 band pass correct with calon-caloff
3 band pass correct (smooth or fit) with data spectra
The default is no bandpass correction
fitbpc: int fit a polynomial of order fitbpc to the masked
version of the band pass correction and use the
polynomial rather than the data to do the bandpass
correction. This is only done if bpc is specified.
smobpc: int smooth the bandpass correction spectra by smobpc channels
before using it to correct the data. The number should be
an odd number of channels. This is only done if bpc is
specified.
blrem: int Remove a polynomial baseline of order blrem. Fit to the
masked portion of the spectra. This is done after
any bandpass correction or averaging.
svd: If baselining is done (blrem) then the default fitting
routine is poly_fit (matrix inversion). Setting svd
will use svdfit (single value decomposition) which is
more robust but slower.
rcv: int If supplied then only process receiver number rcv.
(1-327,2-430,3-610,5-lbw,6-lbn,7-sbw,8-sbh,9-cb,11=xb,
12-sbn)
sav: if keyword set then save the processed data to
ddmmyy.projid.N.sav . The keyword fsav will let you
switch this to a different save file.
fsav: string name for save file. Only used if /sav is set.
file will be forced to end in .sav
RETURNS:
istat: int >=0 number of patterns found
<0 could not process file.
bdat[]: {corget} intensity calibrated data spectra
bcal[]: {corget} intensity calibrated cal spectra
Rsl[nscans]: {sl} return scan list that was used.
nBadIndAr: int number of scans not used because different type.
the indices into sl[] are returned in badindar.
badIndAr[nBadIndAr]: int .. hold indices into sl[] of on/off position
switch data that were not processed because they had
a different data structure. The number of elements will
be in nbadindar
DESCRIPTION:
pfcalib will do the corcalib() processing on all of the
on/calon/caloff patterns in filename. The keyword dir[] will search for
filename in all of the directories in dir.
The processing averages the data scan, converts the data and cal scans
to kelvins using the cals, and optionally bandpass corrects the data.
pfcalib requires that the cal scans follow immediately after the data scans.
See corcalib() for a discussion of the edgefraction,mask and bandpass
correction using fitting or smoothing.
pfcalib() will use the same edgefraction or mask, and bandpass correction
on all of the scans in the file. It also always averages the data scans
to 1 record.
The routine will first scan the file and create a scan list array (see
getsl()). If you have already scanned the file, then using the sl=sl keyword
lets you pass in the scan list array and save the time to scan the file
(the scanlist array used is returned in the rsl=rsl keyword).
The processing then proceeds one pattern at a time. If there are more than
300 records in a scan, then you must use the maxscans= keyword set to the
maximum number of records. If an incomplete on/calon/caloff pattern is found,
the routine will skip forward and continue searching. The processing can be
limited to a subset of the entire file by :
1. scan= keyword. Start processing on this scan.
2. rcv=N keyword. Only process data from receiver number N.
The data and cals are returned in arrays of {corget} structures.
For this to work, all of the data must meet the following criteria:
1. All scans use the same number of boards.
2. For a single board (say board n)
a. All scans should have the same number of polarizations.
b. all scans should have the same number of lags.
(It is ok if different boards have different number of polarizations or
lags).
Any scans whose data structure is different than the first pattern will be
skipped (with a message output to the terminal). The keywords nbadIndAr
and badIndAr will return the number and indices into rsl[] of the patterns
not included because the data type differed from the first pattern found.
The /sav keyword will store the processed data in an idl save file. The
default name is ddmmmyy.projid.n.sav where ddmmmyy.projid.n is taken
from the correlator datafile name. You can change the savefile name
with the fsav= keyword (it still must end in .sav). The variable names
stored in the file will be:
bf[npairs], bcalf[2,npairs]
srcnamesf[npairs]:string An array holding the source names
The save file lets you recall this data into a later session of idl using
restore:
restore,'02may02.x101.1.sav'
SEE ALSO: corcalib()
(See /pkg/rsi/local/libao/phil/Cor2/pfcalib.pro)
NAME:
pfcalonoff - extract calon/off scans and compute cal size and tsys.
SYNTAX: ncals=pfcalonoff(lun,sl,pfcalI,han=han,scanOnAr=scanOnAr,slAr=slAr,$
bla=bla,blda=blda,verb=verb)
ARGS:
lun : int open to file to process
sl[]: {slar} scan list of file from getsl()
KEYWORDS:
han: if set then hanning smooth the data on input
scanOnAr[n]:long if supplied then this is the list of calOn scans to process.
sl is ignored.
slAr[n]:long Output from slar=getsl(lun). When scanOnAr is used sl is
ignored. You can still pass in slAr to speed the file access.
bla: if set then use corblauto and fit to calOn-calOff to remove
outliers
blda: if set then use corblauto and fit to (calOn-calOff)/caloff
to remove outliers
verb : if set to -1 then have corblauto plot the fits and residuals
as it goes along.
otherParms any extra keywords are passed to corcalcmp. These mainly
deal itin will be passed to corcalcmp
they can be used to control the selection of then channels
in each spectra to use.
RETURNS:
ncals : long number of cal pairs found.
pfcalI[ncals]: {} return cal info (defined below)
DESCRIPTION:
Find all of the on/off pairs in the scan list and then process them.
The processing uses corcalcmp() to compute the cal value. It returns
the processed data in pfcalI[].
pfcalI contains:
pfcalI.scan scan number of cal on
pfcalI.az azimuth of cal on
pfcalI.za azimuth of cal on
pfcalI.nbrds number of boards this scan, entries in brdI[]
pfcalI.cfr[8] center frequency for each board (topocentric)
pfcalI.bw[8] bandwidth for each board (Mhz)
pfcalI.pols[2,8] [pols,brds] 0==> not use, 1==>polA, 2--> polB
pfcalI.calScl[2,8] cal scale factore (K/correaltor count)
pfcalI.calval[2,8] cal in kelvins
pfcalI.tsysOn[2,8] Tsys for cal on
pfcalI.tsysOff[2,8] Tsys for cal off
The array in pfcalI are dimensioned for the max number of sbc
(including alfa). Use Use pfcalI.pols[] to decide whether a particular
entry has data (0 --> no data).
Examples:
1. process the entire file. Tell corblauto to use (calOn-calOff)/calOff to
construct the mask of channels to use. Plot these fits.
file='/share/olcor/corfile...'
openr,lun,file,/get_lun
sl=getsl(lun)
ncals=pfcalonoff(lun,sl,pfcalI,/han,/blda,verb=-1)
2. The sl array that is passed in does not have to be the scanlist
from the entire file. You can create a subset and it will process
what it is passed (but the sl passed in must contain the calon,off
scans you want processed). For the example, only process the
source "EA1"..
ind=where(sl.srcname eq 'EA1',count)
if count gt 0 then begin
ncals=pfcalonoff(lun,sl[ind],pfcalI,/han,/blda,verb=-1)
endif
;
; You could have also made the selection on any other entry in sl[]
3. The scanOnAr is another way to process only a subset of the file. It
contains a list of calOn scans to process. In this case the sl argument
is ignored. To speed up the i/o you can still scan the file and
pass the scan list in as slar=slar keyword
Only process cal pair that starts with calOnSCan=517500001L
sl=getsl(lun) ; this is optional
scanOn=517500001l
ncal=pfcalonoff(lun,junk,pfcalI,scanOnAr=scanOn,/blda,verb=-1,slar=sl)
(See /pkg/rsi/local/libao/phil/Cor2/pfcalonoff.pro)
NAME:
pfcorimgonoff - make images of all on/off pairs in a file.
SYNTAX: sl=pfcorimgonoff(lun,sl=sl,ver=ver,clip=clip,col=col)
ARGS:
lun : int logical unit number that points to file.
KEYWORDS:
sl[] : {sl} scan list array. If the user has already created a
scan list of this file (getsl) then you can pass it in
here so the call will not do it again.
ver[2]: float vertical scale (min,max) for the line plot superimposed
on the image: (on/off -1)-median(on/off-1). The default
value is .005 Tsys.
clip[2]: float min/max clipping value for the image (on/off-1).;
The default value is .02 Tsys
col[2]: int columns of image to use for flattening in the time
direction.. 0 through nchannels-1.
RETURNS:sl[]: {sl} return the scan list made by this routine (or the
one passed in).
DESCRIPTION:
The routine will first scan the entire file looking for all of the on/off
pairs. It will then call corimgonoff() for each pair found. After
displaying the image the user will be prompted to continue or quit. It will
return the scan list made by the call. If you recall the routine with the
same file, then you can pass in this sl array to speed up the processing
(it takes a while to scan the file).
Before calling the routine be sure and call xloadct and setup the
lut for greyscale.
EXAMPLES:
xloadct .. then click on bw-linear for greyscale.
openr,lun,'/share/olcor/corfile.07jan02.a1511.2',/get_lun
.. show all of the images using the default settings. Scan the file
the first time.
sl=pfcorimgonoff(lun)
.. rerun the above using channels 900,950 to flatten the image
in the time direction. Rescale the line plot so the maximum is
.002 Tsys (the default was .005). Pass in the sl array we got from
the last call so we don't have to rescan the file.
sl=pfcorimgonoff(lun,sl=sl,ver=[-.002,.002],col=[900,950]
SEE ALSO:
corimgonoff
(See /pkg/rsi/local/libao/phil/Cor2/pfcorimgonoff.pro)
NAME:
pfposonoff - process all position switch onoffs in a file
SYNTAX: npairs=pfposonoff(filename,bf,tf,calsf,han=han,median=median,
scan=scan,dir=dir,sl=sl,maxrecs=maxrecs,
scltsys=scltsys,sclJy=sclJy,
skipscans=skipscans,
rcv=rcv,sav=sav,fsav=fsav,rsl=rsl,$
badIndAr=badIndAr,nbadindar=nbadindar)
ARGS :
filename: string .. file to process
KEYWORDS:
han : if set then hanning smooth the data
median : if set then use median rather mean when averaging a scan.
scan : long if set then start at this scan number
dir[] : string search for filename in these directories.
in this case filename should not have a path.
sl[] : {slist} Scan list to use when searching file. This
saves having to scan the file. Use this to speed
processing if you have already scanned the file (
see the Rsl keyword).
maxrecs : long The maximum number of recs to allow in a scan
the default is 300.
scljy : if set then scale to kelvins using the cals and then
use the gain curve to return the spectra in janskies.
The default is to return the spectra in kelvins.
scltsys : if set then scale to tsys rather than using the cals
to scale to kelvins (the default).
skipScans[] : long Any on position scan numbers found in the array
skipScans will not be processed.
rcv : int If supplied then only process reciver number rcv.
(1-327,2-430,3-610,5-lbw,6-lbn,7-sbw,8-sbh,9-cb,11=xb,
12-sbn)
sav : if keyword set then save the processed data
to ddmmyy.projid.N.sav . The keyword fsav
will let you switch this to a different save file.
fsav : string name for save file. Only used if /sav is set.
file will be forced to end in .sav
RETURNS:
npairs: int number of complete on/offs processed
-1 if illegal filename
-1 if scan option but scan not in file
bf[npairs] : {corget}.. return (on/off-1) average for each pair
scaled to: kelvins,janskies, or Tsys units.
tf[npairs] : {temps} .. array of temp structures (see corposonoff).
calsf[nbrds,npairs]: {cal} .. array of cal structures (see corposonoff).
Rsl[nscans]: {sl} .. return scan list that was used.
nBadIndAr : int .. number of scans not used because different type.
the indices into sl[] are returned in badindar.
badIndAr[nBadIndAr]: int .. hold indices into sl[] of on/off position
switch data that were not processed
because they had a different data structure.
The number of elements should be in nbadindar
DESCRIPTION:
pfposonoff will process all of the on/off pairs it finds in filename.
The keyword dir[] will search for this file in all of the directories in dir.
The routine will first scan the file and create a scan list array (see
getsl()). If you have already scanned the file, then using the sl=sl keyword
lets you pass in the scan list array and save the time to scan the file
(the scanlist array used is returned in the rsl=rls keyword).
If an incomplete on/off pair is found, it will skip forward and continue
searching. The processing can be limited to a subset of the entire file by :
1. scan= keyword. Start processing on this scan.
2. rcv=N keyword. Only process data from receiver number N.
The processed on/off-1 spectra will be returned in a array bf[npairs] of
{corget) structures. For this to work, all of the position switch data must
meet the following criteria:
1. All scans use the same number of boards.
2. For a single board (say board n)
a. All scans should have the same number of polarizations.
b. all scans should have the same number of lags.
(It is ok if different boards have different number of polarizations or
lags).
The keywords nbadIndAr and badIndAr will return the number and indices
into rsl[] of the patterns not included because the data type differed
from the first one found.
By default the units for the returned data is Kelvins. Setting
the keyword scljy will return the spectra in janskies (as long as there is
a gain curve for this receiver). Setting the keyword sclTsys will return
the spectra in units of Tsys.
The /sav keyword will store the processed data in an idl save file. The
default name is ddmmmyy.projid.n.sav where ddmmmyy.projid.n are taken
from the correlator datafile name. You can change the savefile name
with the fsav= keyword (it still must end in .sav). The variable names
stored in the file will be:
bf[npairs], tf[npairs], calsf[4,npairs] .. same as the returned data.
savefile :string The name of the save file used.
srcnamesf[npairs]:string An array holding the source names
The save file lets you recall this data into a later session of idl using
restore:
restore,'02may02.x101.1.sav'
You can also use corsavbysrc() to recall a set of these save files
and create arrays of processed data by src.
SEE ALSO: corposonoff, corsavbysrc
(See /pkg/rsi/local/libao/phil/Cor2/pfposonoff.pro)
NAME:
pfposonoffrfi - process all position switch onoffs in a file
SYNTAX: npairs=pfposonoffrfi(filename,bf,han=han,
scan=scan,dir=dir,sl=sl,maxrecs=maxrecs,
scltsys=scltsys,sclJy=sclJy,
skipscans=skipscans,
rcv=rcv,sav=sav,fsav=fsav,rsl=rsl,$
badIndAr=badIndAr,nbadindar=nbadindar,$
arSecPerChn=arSecPerChn,arRmsbyChn=arRmsbyChn)
ARGS :
filename: string .. file to process
KEYWORDS:
han : if set then hanning smooth the data
scan : long if set then start at this scan number
dir[] : string search for filename in these directories.
in this case filename should not have a path.
sl[] : {slist} Scan list to use when searching file. This
saves having to scan the file. Use this to speed
processing if you have already scanned the file (
see the Rsl keyword).
maxrecs : long The maximum number of recs to allow in a scan
the default is 300.
scljy : if set then scale to kelvins using the cals and then
use the gain curve to return the spectra in janskies.
The default is to return the spectra in kelvins.
scltsys : if set then scale to tsys rather than using the cals
to scale to kelvins (the default).
skipScans[] : long Any on position scan numbers found in the array
skipScans will not be processed.
rcv : int If supplied then only process reciver number rcv.
(1-327,2-430,3-610,5-lbw,6-lbn,7-sbw,8-sbh,9-cb,11=xb,
12-sbn)
sav : if keyword set then save the processed data
to ddmmyy.projid.N.sav . The keyword fsav
will let you switch this to a different save file.
fsav : string name for save file. Only used if /sav is set.
file will be forced to end in .sav
RFI EXCISION keywords:
smorfi: long When searching for rfi, smooth the data by
smorfi channels before searching. smorfi should
be an odd number. This is only used for searching
for the rfi, the returned data is not smoothed.
flatTsys: If set, then divide each spectra by
its median value before searching each freq channel
along the time direction. This will allow records with
total power fluctuations to be included.
adjKer[i,j]: float convol this array with the mask array for each
sbc. If the convolved value is not equal to the
sum of adjKer then also exclude this point. This
array lets you exclude points adjacent to bad
points (xdim is frequency, ydim is time). The
array dimensions should be odd.
frqChnNsig: float The number of sigmas to use when excluding a point.
linear fits are done by channel along the time axis.
Points beyond frqChnNsig will be excluded from the
averaging.
badNsig: float The rms/median is computed for each freq
channel (after excluding bad points). If badNsig is
provided, then any freq channels with rms/median
greater than 3. will be refit along the time direction
searching for bad points using badNsig rather than
frqChnNSig as the threshold for bad points.
An example would be frqChnSig=3, badNsig=2.
tpNsig: float The total power is computed vs time for the points
that pass the freqchannel test. A linear fit is
then done versus time. Any time points whose
residuals are greater than tpNSig will be ignored.
The default is 0 (this computation is not done)
RETURNS:
npairs: int number of complete on/offs processed
-1 if illegal filename
-1 if scan option but scan not in file
bf[npairs] : {corget}.. return (on/off-1) average for each pair
scaled to: kelvins,janskies, or Tsys units.
Rsl[nscans]: {sl} .. return scan list that was used.
nBadIndAr : int .. number of scans not used because different type.
the indices into sl[] are returned in badindar.
badIndAr[nBadIndAr]: int .. hold indices into sl[] of on/off position
switch data that were not processed
because they had a different data structure.
The number of elements should be in nbadindar
arSecPerChn[npairs]:{corget} Holds the number of secs avged for each freq channel
Use this as the weights for each freq bin
when combining different patterns.
arRmsbyChn[npairs]:{corget} rms/mean along each channel. 1 for each
pattern.
DESCRIPTION:
pfposonoffrfi will process all of the on/off pairs it finds in filename.
The keyword dir[] will search for this file in all of the directories in dir.
This routine tries to excise rfi before combining records. It calls
corposonoffrfi (rather than corposonoff).
The routine will first scan the file and create a scan list array (see
getsl()). If you have already scanned the file, then using the sl=sl keyword
lets you pass in the scan list array and save the time to scan the file
(the scanlist array used is returned in the rsl=rls keyword).
If an incomplete on/off pair is found, it will skip forward and continue
searching. The processing can be limited to a subset of the entire file by :
1. scan= keyword. Start processing on this scan.
2. rcv=N keyword. Only process data from receiver number N.
The processed on/off-1 spectra will be returned in a array bf[npairs] of
{corget) structures. For this to work, all of the position switch data must
meet the following criteria:
1. All scans use the same number of boards.
2. For a single board (say board n)
a. All scans should have the same number of polarizations.
b. all scans should have the same number of lags.
(It is ok if different boards have different number of polarizations or
lags).
The keywords nbadIndAr and badIndAr will return the number and indices
into rsl[] of the patterns not included because the data type differed
from the first one found.
By default the units for the returned data is Kelvins. Setting
the keyword scljy will return the spectra in janskies (as long as there is
a gain curve for this receiver). Setting the keyword sclTsys will return
the spectra in units of Tsys.
The /sav keyword will store the processed data in an idl save file. The
default name is ddmmmyy.projid.n.sav where ddmmmyy.projid.n are taken
from the correlator datafile name. You can change the savefile name
with the fsav= keyword (it still must end in .sav). The variable names
stored in the file will be:
bf[npairs], .. same as the returned data.
arSecPerChn[npairs]
arRmsByChn[npairs]
savefile :string The name of the save file used.
srcnamesf[npairs]:string An array holding the source names
The save file lets you recall this data into a later session of idl using
restore:
restore,'02may02.x101.1.sav'
You can also use corcmbsav() to recall a set of these save files
and create arrays of processed data by src.
SEE ALSO: corposonoffrfi, corcmbsav
(See /pkg/rsi/local/libao/phil/Cor2/pfposonoffrfi.pro)
NAME:
sl_mkarchive - create scan list for a set of files
SYNTAX: nscans=sl_mkarchive(listfile,slAr,slfilear,logfile=logfile)
ARGS:
listfile: string. filename holding names of files to scan
(one per line). semi-colon as the first char is a comment.
KEYWORDS:
logfile: string name of file to write progress messages. by default just
goes to stdout.
RETURNS:
slAr[nscans] {sl} One large scanlist array for all the files
slfileAr[m]: {slInd} One entry per file showing where each file starts/
stops in slAr.
nscans : long number of scans found
KEYWORDS:
DESCRIPTION:
Calls getsl() for every corfile filename in listfile. It creates one large
scanlist array for all of the scans. It also creates an index array telling
where each file starts and stops in the slar.
The slAr[] is an array of structures (1 per scan) containing:
scan : 0L, $; scannumber this entry
bytepos : 0L,$; byte pos start of this scan
fileindex : 0L, $; lets you point to a filename array
stat : 0B ,$; not used yet..
rcvnum : 0B ,$; receiver number 1-16
numfrq : 0B ,$; number of freq,cor boards used this scan
rectype : 0B ,$;1-calon,2-caloff,3-posOn,4-posOff
numrecs : 0L ,$; number of groups(records in scan)
freq : fltarr(4),$;topocentric freqMhz center each subband
julday : 0.D,$; julian day start of scan
srcname : ' ',$;source name (max 12 long)
procname : ' '};procedure name used.
The slFileAr is an array of structures (1 per corfile):
path : '' , $; pathname
file : '' , $; filename
i1 : 0L , $; first index into slAr array for this file
i2 : 0L } ; last index into slAr array for this file
EXAMPLE:
If the listfile contained:
; a comment
/share/olcor/calfile.19mar01.a1400.1
/share/olcor/calfile.20apr01.a1446.1
/share/olcor/calfile.20apr01.a1489.1
/share/olcor/calfile.20mar01.a1389.1
It would process all 4 files and return the data in slAr, and slfilear.
NOTE:
This routine is used to generate the save files that hold slar,slFileAr
by month. arch_gettbl uses these save files.
SEE ALSO:arch_gettbl,arch_getdata
(See /pkg/rsi/local/libao/phil/Cor2/sl_mkarchive.pro)
NAME:
sl_mkarchivecor - create cor scan list for archive
SYNTAX: nscans=sl_mkarchivecor(slAr,slFileAr,slcor,slfilecor,verbose=verbose,$
logfile=logfile
ARGS:
slAr[n]: {sl} sl array returned by arch_gettbl
slFileAr[m]: {slind} returned by arch_gettbl
RETURNS:
slCor[nscans] {corsl} Extended scanlist containing ra,dec and corinfo
slFileCor[m]: {slind} file array. same as slfilear unless some files
have disappeared between running slar,slcorar.
KEYWORDS:
verbose: If set then print each file as we process it.
logfile: string name of file to write progress messages. by default
messages only go to stdout.
DESCRIPTION:
sl_mkarchive creates a large scanlist for a set of files (usually a
months worth of data. These files are stored on disc and arch_gettbl
allows you to retrieve this data or a subset of it. The info in slAr
is a limited set that contains info common to all scans at ao (correlator
ri, radar, atm, etc..).
The sl_mkarchivecor will take the slAr and make a superset of this
structure. In addition to the {sl} info it will contain ra,dec info as well
as correlator info. You can then access this info with the arch_gettbl
command using the /cor keyword.
The slcor[] is an array of structures (1 per scan) containing:
scan : 0L, $; scannumber this entry
bytepos : 0L,$; byte pos start of this scan
fileindex : 0L, $; lets you point to a filename array
stat : 0B ,$; not used yet..
rcvnum : 0B ,$; receiver number 1-16
numfrq : 0B ,$; number of freq,cor boards used this scan
rectype : 0B ,$;1-calon,2-caloff,3-posOn,4-posOff
numrecs : 0L ,$; number of groups(records in scan)
freq : fltarr(4),$;topocentric freqMhz center each subband
julday : 0.D,$; julian day start of scan
srcname : ' ',$;source name (max 12 long)
procname : ' ',# ;procedure name used.
.. the extended cor info begins here...
projId : '' ,$; from the filename
patId : 0L ,$; groups scans beloging to a known pattern
secsPerrec : 0. ,$; seconds integration per record
channels : intarr(4),$; number output channels each sbc
bwNum : bytarr(4),$; bandwidth number 1=50Mhz,2=25Mhz...
lagConfig : bytarr(4),$; lag config each sbc
lag0 : fltarr(2,4),$; lag 0 power ratio (scan average)
blanking : 0B ,$; 0 or 1
azavg : 0. ,$; actual encoder azimuth average of scan
zaavg : 0. ,$; actual encoder za average of scan
raHrReq : 0.D,$; requested ra , start of scan J2000
decDReq : 0.D,$; requested dec, start of scan J2000.
Delta end-start real angle for requested position
raDelta : 0. ,$; delta ra last-first recs. Amins real angle
decDelta : 0. ,$; delta dec (last-frist)Arcminutes real angle
azErrAsec : 0. ,$; avg azErr Asecs great circle
zaErrAsec : 0. $; avg zaErr Asecs great circle
The routine sl_mkarchive is run monthly and the info is stored in
files by month. This routine can then be run to create the slcor data.
EXAMPLE:
The monthly processing would look like:
;
;1. create listfile from corflist.sc 010101 010131 excl.dat outside of idl.
;2. make the slAr, slFileAr and save it
nscans=sl_mkarchive(listfile,slAr,slfilear)
save,sl,slar,file='/share/megs/phil/x101/archive/sl/f010101_010131.sav'
;
;3. now make slcor ar.
arch_gettbl(010101,010131,slAr,slFilear)
nscans=sl_mkarchivecor(slar,slfilear,slcor,slfilecor)
slar=slcor ; so they have the same names
slfileAr=slfilecor ; so they have the same names
save,slar,slfilear,file=$
'/share/megs/phil/x101/archive/slcor/f010101_010131.sav'
SEE ALSO:sl_mkarchive,arch_gettbl,arch_getdata
(See /pkg/rsi/local/libao/phil/Cor2/sl_mkarchivecor.pro)
NAME:
wascheck - check if this call is for was data
SYNTAX: istat=wascheck(lun,file=file,b=b)
ARGS:
lun/desc: int or {wasdesc} ..
KEYWORDS:
b: {corget} if provided then this is a correlator buffer read in
with corget(). The was test will be on this buffer rather
than on the lun/desc. It tests for the existence of
b[0].(0).hf (fits header portion of the structure).
file: char if provided, then ignore lun and see if this file
ends in .fits.
RETURNS: 1 this is a was call (lun is a structure or file ends in .fits)
0 lun is an int, must be cor data
DESCRIPTION:
Return 1 if the lun is a was descriptor rather than a int. Or if
file ends in .fits
(See /pkg/rsi/local/libao/phil/Cor2/wascheck.pro)