NAME:
aomasexamples - Using the idl mas (mock spectrometer)
SYNTAX: none
At the Arecibo observatory (AO) the mock spectrometer run in
spectral line mode creates AO fits files (also known as cima fits V2.0). The
format is based on the sdfits format (with a few additions). Idl routines
have been written to access and process this data. All of these routines
begin with mas (Mock Arecibo Spectrometer).
Terminology
bm or box: there are 7 mock boxes in a spectrometer set.
-For alfa each of these is mapped to a beam.
-For single pixel observing these are mapped to different frequency bands.
band : each box has two separate 172 Mhz subbands.
-for alfa these two bands cover the 300 Mhz of alfa for a beam.
-for single pixel observing only band 1 is used.
group : there are two complete sets (or groups) of 7 boxes.
Box N of each group gets the same input IF. The different
groups can then vary the resolution, bandwidth, or integration time.
The data files:
The file naming convention is :
projId.yyyymmdd.bMsNgP.nnnnn.fits where:
projid : is the project id entered when the datataking gui (cima) is started.
yyyymmdd: is the date when the data file was openned (AST).
bMsNgP :bM,M=0..6 This is the box/bm number for this file b0..b6
:sN,N=0,1 The subband within each box. For single pixel observing
this is always s1. For alfa, this can be 0, or 1 for the lower or
upper 172 Mhz band.
:gP,P=0,1 This is the group number g0 or g1 for this file.
nnnnnn : 5 digit number. This gets incremented by 100 for each scan. If
a single scan requires more than 1 file (> 2gb) then the next
file is nnnnn +1.
.fits fits suffix.
; for more info on the fits file structure see:
; http://www.naic.edu/~phil/software/datataking/fits/fit_header.html
;
; Each scan starts a new datafile. If you do a position switch observation
;followed by a cal on, cal off you will end up with 4 files: position on,
;position off, calOn, caloff.
File locations:
The online datafiles are written to:
/share/pdataM/pdev/ where M= boxNumber+1 (it goes 1..7 sorry about that).
They may eventually be moved to /share/projid/xxx where projid is the
projid of the file (this does not always happen).
Starting idl:
To process the mas datasets:
idl
@phil (or whatever you need to add your base directory)
@masinit ..initialize for mas data. This calls @geninit and then
adds the ./mas directory path to idl's search path.
-----------------------------------------------------
NOTE: These routines use the goddard binary fits table routines to
access the file. At AO they are in the directory:
/pkg/rsi/idl/lib/locallib/astron/pro/fits_bintable. If you download
the AO idl routines, you will also need to get a copy of the
goddard routines and add the fits_bintable directory to the path.
;------------------------------------------------------------------------
SOME BASICS:
; open a file
;
file='/share/pdata1/pdev/x106_5.20100819.b0s1g0.00000.fits'
istat=masopen(file,desc)
;
; list contents of the file
;
maslist,desc
;
;
; read a row full of data (1 or more spectra):
; first rewind the file since maslist() read till the end of file
rew,desc
istat=masget(desc,b)
;
; look at the b data structure:
help,b
** Structure <88f81dc>, 8 tags, length=656472, data length=656467, refs=1:
H STRUCT -> MASFHDR Array[1]
NCHAN LONG 8192
NPOL LONG 2
NDUMP LONG 10
BLANKCORDONE INT 0
ST STRUCT -> PDEV_HDRDUMP Array[10]
ACCUM DOUBLE 0.0000000
D LONG Array[8192, 2, 10]
; Note that:
; nchan=8192 (8192 freq channels),
; npol=2 (2 polarizations)
; ndump=10 (there are 10 spectra in this row).
; d[8192,2,10] the data is dimensioned 8192 chan,by 2 pols, by 10 dumps)
;
;look at fits header for this row:
help,b.h,/st
IDL> help,b.h,/st
** Structure MASFHDR, 133 tags, length=888, data length=885:
TDIM1 STRING '(8192,1,1,2,10)'
TDIM2 STRING '(10,1,1,1,10)'
OBJECT STRING 'STOPPED'
CRVAL1 DOUBLE 1.2000000e+09
CDELT1 DOUBLE 21000.000
CRPIX1 DOUBLE 4097.0000
CRVAL2 DOUBLE 243.63913
CRVAL3 DOUBLE 15.457966
CRVAL4 DOUBLE -56.000000
CRVAL5 DOUBLE 79762.000
CDELT5 DOUBLE 0.10000000
AZIMUTH DOUBLE 285.42697
ELEVATIO DOUBLE 79.973986
; ... continues ...
;
; Plot the spectra
;
masplot,b
;
; this will overplot the 10 .1 second spectra in the row.
; red is polA, green is polB.
;
; look at the documentation for the idl routines:
explain,masdoc
----- Documentation for /pkg/rsi/local/libao/phil/doc/masdoc.pro -----
;NAME:
masdoc - routine list (single line)
gftget - input next galfacts timedome row from disc
gftgetfile - input an entire file
gftgetstat - input status info from file
gftopen - open galfacts decimation in time file.
masaccum - accumulate an array of buffers
masavg - read and average spectra
masavgmb - read and average spectra for multi beams
mascalonoff - compute cal scl factor from calon,caloff
; continues..
;
; look at the masplot routine:
;
explain,masplot
----- Documentation for /pkg/rsi/local/libao/phil/mas/masplot.pro -----
;NAME:
masplot - plot mas spectra
SYNTAX: masplot,b,freq=freq ,over=over,pollist=pollist,norm=normRange,off=off,$
retvel=retvel,restFreq=restFreq,velCrdsys=velCrdSys,$
mfreq=mfreq,colar=colar
ARGS:
b[n]: {b} mas structure from masget to plot
KEYWORDS:
freq[n]: float if provided then use this as the x axis
over: if set then overplot spectra
pollist: long pols to plot: 12 = polA,B, -1--> all pols available
norm[2]: float normalize the spectra. If two element array then is is the range
of frequency bins to use for the median (cnt from 0). If a single
element or /norm then use all the bins.
off: float If multiple spectra then increment each spectra by off.
chn : if set then plot vs channel number (count from 0)
smo : int smooth by this many channels
; continues..
;------------------------------------------------------------------------
; Summary of different routines:
;
; Open a file: masopen()
; close a file: masclose()
; Input a single row : masget()
; postion to a row in a file: maspos()
; input an entire file: masgetfile()
; plot one or more spectra: masplot()
; find a set of files: masfilelist()
; .. masonofffind() find position switch files
; .. masdpsfind() find double position switch files
;
; On off position switching:masposonoff()
; double position switching:masdpsp()
; Calibrate stokes data: masstokes()
;
; Accumulate spectra (after readin:):masaccum()
; Perform arithmetic on spectra: masmath()
; Compute rms/mean by channel for an array of spectra: masrms()
;
;------------------------------------------------------------------------
; automating file access:
; The datataking generates lots of files. You can automate the processing
; of these files used masfilelist(). This will select a subset of files
; on disc given the date, projid, beam, band, group, etc. It returns
; an array of structures (fnmI[] filenamem info structures).
; You can then pass elements of this array to masopen() and some other
; routines to automate the processing.
; A second routine:
; masfilesum() will take an fnmI[] array and read the headers for each of
; these files. You could then fine tune the selection process using the
; where() routine of idl:
; EG:
; find the a2489 files on 09oct09, all beams,band=1:
; n=masfilelist('',fnmI,proj='a2489',yymmdd=20091009,band=1,/appbm)
; n returns 448 files found
IDL> help,fnmI,/st
** Structure MASFNMPARS, 9 tags, length=64, data length=62:
DIR STRING '/share/pdata1/pdev/'
FNAME STRING 'a2489.20091009.b0s1g0.00000.fits'
PROJ STRING 'a2489'
DATE LONG 20091009
SRC STRING ''
BM INT 0
BAND INT 1
GRP INT 0
NUM LONG 0
;
; read the headers from all the bm =0 files
;
ii=where(fnmI.bm eq 0,cnt)
fnmI0=fnmI[ii]
nsum=masfilesum('',sumI,fnmI=fnmI0,/list)
;
; select the on source scans:
ii=where((sumI.h.obsmode eq 'ONOFF') and (sumI.h.scantype eq 'ON'),cnt)
;
; The above was just a demo for explaining things..
; All of the above can also be done using the routine masonofffind()
;
(See /pkg/rsi/local/libao/phil/mas/aomasexamples.pro)
NAME:
gftget - input next galfacts timedome row from disc
SYNTAX: istat=gftget(des,bon,boff,row=row,hdronly=hdronly)
ARGS:
desc:{descmas} from dftopen();
RETURNS:
b: structure holding the hdr and data input
istat: 1 ok
: 0 hiteof
:-1 i/o error..bad hdr, etc..
KEYWORDS:
row : if set then position to row before reading (count from 1)
if row=0 then ignore row keyword
hdronly : if set then just return the row header. no status or data.
DESCRIPTION:
Read the next row from a galfacts time domain decimation fits datafile pointed to by desc.
If keyword row is present, position to row before reading.
(See /pkg/rsi/local/libao/phil/mas/gftget.pro)
NAME:
gftgetfile - input an entire file
SYNTAX: istat=gftgetfile.pro(fname,bon,boff,nrows=nrows,hdronly=hdronly)
ARGS:
fname: string file to open
KEYWORDS:
hdronly: if set then just return the headers.
RETURNS:
istat: 1 got all the requested records
: 0 returned no rows
: -1 returned some but not all of the rows
bon[nrows]: {} array of structs holding the calon data
boff[nrows]: {} array of structs holding the caloff data
nrows : number of rows returned
(See /pkg/rsi/local/libao/phil/mas/gftgetfile.pro)
NAME:
gftgetstat - input status info from file
SYNTAX: istat=gftgetstat(desc,statOn,statOff,nrows=nrows,all=all)
ARGS:
desc: {} struct returned from masopen()
KEYWORDS:
row : long row to position to before reading count from 1.
nrows: long number of rows of stat data to read in .default=1 row
from current position.
all: if set then read in the stat info for the entire file.
RETURNS:
istat: n number of stat stuctures returned. One struct for each
spectra. if nrows were read, there will be nrows*spcPerRow
-1 error reading file
stat[n]:{stat} stat data from file.
(See /pkg/rsi/local/libao/phil/mas/gftgetstat.pro)
NAME:
gftopen - open galfacts decimation in time file.
SYNTAX: istat=gftopen(filename,desc,fnmI=fnmI,hdr=hdr)
ARGS:
filename: string filename to open (unless fnmI is specified)
fnmI : {] returned from masfilelist. if provided,
then ignore filename and use this as the
file to open.
KEYWORDS:
RETURNS:
istat: 0 ok
-1 could not open file.
desc : {} file descriptor to pass to the i/o routines.
hdr : [string] if present then return fits bintable header for the fille
(See /pkg/rsi/local/libao/phil/mas/gftopen.pro)
NAME:
masaccum - accumulate an array of buffers
SYNTAX: naccum=masaccum(bIn,baccum,avg=avg,scl=scl,new=new,mb=mb)
ARGS:
bIn[n]: {} stuctures to accumulate
KEYWORDS:
avg: if set then average when done
scl[n]: if provided then scale the bIn data by scl before
adding. This can be used to weight data by g/t.
scl can be a single number (then every spectra is
scaled by this
scl[n]: then every entry of bIn is scaled by its own
value.
if /mb is used then scl can be dim 1 or dim baccum[]
(which equals the first dim of Bin).
Note that polA,B of one bIn[i] always get the same scale
value.
mb: if set then bIn[nbm,nrecs] is a multi beam array. The 1st
dimension is the number of beams, the 2nd dimension are the
records to accumulate. In this case baccum[nbm] and
only the 2nd dimension of bIn is accumulated.
new: if set then allocate baccum. This is the first call
if not set then baccum passed in will be added to.
double: if set then make accum array double. needs to be used
when /new is called.
RETURNS:
naccum: > 0 number of accumulatins we accumulated in baccum (this does
not include scl).
: 0 none summed
baccum[nbm]: {} the accumulated spectra
The overflow status bits inf baccum.b.st will hold the max
overflow count for the records that were included in the accumulation.
DESCRIPTION:
Accumulate 1 or more records worth of data into baccum. If keyword
/new is set then allocate baccum before adding. The header values in
baccum will come from the first record added into baccum.
Each element of bin[i] will be added to baccum 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 masavg has been called on an accumlated bacccum.
It divides by bacccum.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 masplot is called, it will scale the data by 1./badd.b1.accum. before
plotting.
When calling masaccum with the new keyword, you can include the
/mb keyword. This will allocate baccum to be the same dimension as
the first dimension of binp. All future calls using baccum will add binp element
wise to baccum. This can be used when accumulating multiple maps.
Accumulated data must be of the same type (numlags, numbsbc, bw,etc..).
Example:
print,masget(lun,b)
masaccum,b,badd,/new
print,masget(lun,b)
coraccum,b,badd
masplot,badd
; Add n scans together element wise:
for i=0,n-1 do begin
print,masgetfile(lun,b)
masaccum,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 masgetfile).
print,masgetfile(lun,b,scan=scan)
masaccum,b,bsum,/new
masplot,bsum
;
(See /pkg/rsi/local/libao/phil/mas/masaccum.pro)
NAME:
masavg - read and average spectra
SYNTAX: istat=masavg(desc,navgspc,b,bIn=bIn,row=row,toavg=toavg,double=double)
ARGS:
desc: {} returned by masopen
navgspc: long number of averaged spectra to return
KEYWORDS:
bIn[]:{} if suppplied then take input data from bIn rather than
reading from disc (desc is ignored).
row: long row to position to before reading (cnt from 1)
toavg: long number of spectra to avg
double: if set then avg data as doubles. default is float.
RETURNS:
istat: 1 returned the requested number of averaged spectra
: 0 returned no average spectra
: -1 returned some but not all of the rows
b[navgspc]: {} array of structs holding the averaged data
DESCRIPTION
Read and average spectra. Each spectra will average toavg spectra
(if not supplied then average 1 spectra). Continue reading and averaging
spectra until navgspc averaged spectra have been done.
If blanking is enabled then scale the averaged spectra by the
number of ffts (so that all averaged spectra have the same mean value).
In this case the b.st.fftaccum field will hold the fft's accumulated
for the first spectra of the average (can't have routine the correct value
since this is a short rather than a float).
(See /pkg/rsi/local/libao/phil/mas/masavg.pro)
NAME:
masavgmb - read and average spectra for multi beams
SYNTAX:navgspcR=masavgmb(yymmdd,projid,filenum,b1,b2,$
toavg=toavg,navgspc=navgspc,row=row,grp=grp,band=band)
ARGS:
yymmdd: long date to process
projid: char string on first part of filename.
filenum: long the filenumber to process
KEYWORDS:
toavg: long number of spectra to avg. Default is the entire file
navgspc: long Number of avg spectra to return. Only used if toavg=
is specified. By default return all the averaged spectra
from the file. If requested number of averaged spectra
is greater than the number available, just return
the number we could avg.
row: long row to position to before reading (cnt from 1). Default
is start of file.
RETURNS:
navgspcR: long The number of average spc returned for each beam,band
: 0 did not find any files matching request
: -1 error accessing a file. no data returned
b1[navgspcR,nbeams]: {} array of structs holding the averaged data
for the lower 170 Mhz band in the first IF. This is usually
the higher 170 Mhz at RF since a high side lo flips the band.
b2[navgspcR,nbeams]: {} array of structs holding the averaged data
for the higher 170 Mhz band in the first IF.
(See /pkg/rsi/local/libao/phil/mas/masavgmb.pro)
NAME:
masazswinginp - input an azimuth swing.
SYNTAX: istat=masazswinginp(fnmI,fnum,azStart,azEnd,dat,azAr,freqAr,hwc=hwc)
ARGS:
fnmI[n]: {} hold info for files to search
fnum : long file number to process (one or more beams)
azStart : float deg. starting az
azEnd : float deg ending az
KEYWORDS:
hwc : if set then data has hardware winking cal. just return the cal Off data
Returns:
istat : -1 trouble
: >0 number of spectra returned
dat[nchan,nbms,npnts]:float spectra for swing
azAr[npnts] : azimuth for each spectra
freqAr[nchn] : for the spectra
DESCRIPTION:
input a mock azswing, return spectra and azpos for each spectra.
Also return the freqar.
The user inputs fnmI that holds the files to search for fnum.
Fnum is the file number to use. there can be more than 1 file, if
multiple beams were used.
Azstart,azend are the start, stop azimuth to use for the swing.
This should be in the order the data was taken,
On return, the spectra will be returned in increasing freq order
(as will freqAr). The center channel will be interpolated across (to get rid of dc).
The azAr and the order of the dat will be returned in increasing
azimuth order..
12meter notes: 29mar22
--> not needed since hdr_enc_az has the negative 51G
-reported azimuths are positive 0..359.000
-If you use a negative number for azstart or azend then
is will map the measured azimuth to -180 to 360.
eg. If azstart < 0 then it will map az[0] to 0 deg to -az (assumed CW)
if azEnd <0 then it will search for az 0 in the points and all folloing points will az-360.
kludge until a change 12m program to return the -180 to 360.
- added nudump
(See /pkg/rsi/local/libao/phil/mas/masazswinginp.pro)
NAME:
mascaldatafind - find patterns with cals and data
SYNTAX: n=mascaldatafind(projId,sumI,patI,yymmdd=yymmdd,appbm=appbm,
dirI=dirI,bm=bm,band=band,grp=grp,useSumI=useSumI)
ARGS:
projid: string project id to search for '.*' for all proj
[3]=2 cal off has to match cnt in cal on
KEYWORDS:
useSumI: int if true then user supplies sumI
selection keyords
yymmdd: long yyyymmdd limit to this date
appbm: apply bm number to each directory. This should normally
be set:
dirI[2]: string if data not in default /share/pdataN/pdev
see masfilelist for usage.
bm : int if supplied limit to this beam
band : int 0,1 limit to this band. Note if you have single pixel
data you probably should set band=1 or masfilelist may
not find the files.
grp : int limit to this group 0, or 1.
RETURNS:
n :int number proccessed patterns found. Each pointing pattern
will have 1 entry
multple bms.
patI[n]: info on processed patterns
useI[m]: summary info
patI STRUCTURE:
PATID LONG 34700117
CALONSCAN LONG 34700118
CALOFFSCAN LONG 34700119
DATAOBSMODE STRING 'FIXEDAZ'
DATASCAN LONG 34700117
DATAFILENUM1 LONG 200
NDATAFILES LONG 1
DATAROWS LONG 1200
DATADMPSROW LONG 1
USEGRP INT Array[2]
USEBANDBMGRP INT Array[2, 7, 2]
DESCRIPTION:
Find all of the patterns that have:
1. a data scan
2. a cal on ,off scans
3. and meet the keyword specs
The routine accepts the same parameters at masfilelist:
yymmdd=,dirI=,bm=,band=,grp=num=,appbm=appbm
The routine will pass back an array of patI structs holding info on
each pattern that meets the criteria. The info includes the
patid and and and numbers of the cals, data. You can then use these
values to search through the fnmI struct that is returned and then read the
files.
EXAMPLE:
1. Find the patterns with cals for 20101213:
yymmdd=20101213
projid='.*'
npat=mascaldatafind(projId,fnmI,patI,yymmdd=yymmdd,/appbm)
yymmdd=yymmd20100930,grp=0,band=1,/appbm)
2. recall mascaldatafind with user supplying fnmI
n=masfilelist('',fnmI, params)
npat=mascaldatafind(projId,fnmI,patI,/usefnmI)
(See /pkg/rsi/local/libao/phil/mas/mascaldatafind.pro)
NAME:
mascalonoff - compute cal scl factor from calon,caloff
SYNTAX: istat=mascalonoff(bon,boff,calI,edgeFract=edgeFract,$
usecalI=usecalI,mask=mask,cmpmask=cmpmask,okind=okind)
ARGS :
bon[n]:{}: mas struct hold cal on data.
If n > 1 then each record is assumed to have the same header
setup (nchan,freqrange,npol,etc)
the cntsToK conversion factor will be averaged over all
of these records.
boff[n]:{}: mas struct holding cal off data. boff[n] must have n = to
bon[n]
calI{} : struct holding cal value info and conversion factors.
Normally it is passed back to the caller. If useCalI
is set then the calValue info is passed in via calI.
The scale factors calI.cntsToK are always passed back
in calI.
KEYWORDS:
edgeFract[2]: float fract of the bandpass edge to ignore when averaging.
default: .06
1 entry: use both for low,high freq side
2 entery: [0] for lowFreq,[1] for high freq side
useCalI:int if true then use calI parameter for the cal value info.
If you are doing more than one call with the same
setup then this lets you skip the mascalval() call
after the first call.
mask[nchan]: int mask to use for computing cals. !=0 --> use this channel.
Should be the same order as the spectra eg
mask[0] corresponds the bon.d[0,xx]
cmpmask : int if true then compute the mask useing calon/caloff
It will also exclude echgeFract points from each
side (def to 6%). If mask= keyword is also present
then combine mask passed in with the bad channels we found
from cmpmask.
okInd : int if set then calI.indused is changed to and array of nchan.
It will then have a 1 for
channels used, 0 for channels skipped. This will allow you to
put calI into an array (if nchans is fixed).
If you use okInd then
you can not later pass okInd into this routine as calI (
since it needs the inides used to compute the calvalue.)
RETURNS:
istat: int 0 finished ok, -1 error
calI : {} returns calInfo and conversion factors calI.cntsToK[2]
** Structure , 9 tags, length=64, data length=64, refs=1:
CALVAL FLOAT Array[2] ; averaged cals pola,b
exposure float 0 ; integration time that cntstok corresponds to.
NUMPOL INT 2 ; number of pols 1 --> stokes I
CNTSTOK FLOAT Array[2] ; polA,b conversion cnts to Kelvins
NPNTS LONG 4096 ; if (npnts<0) then okInd was specified and
indused[nchan] has 1 for used,0 for not used
if npnts>= then indused is the indices used
for avg
FLIPPED BYTE 1 ; 1 --> spectra is flipped
EDGEFRACT FLOAT Array[2] ; edgeFract[0,1] fract to keep..spc order
USEMASK BYTE 1 ; 1 --> used mask,0--> used edgefract
INDUSED LONG Array[4096]; if okInd=0 (or not supplied) then:
indices into spc for channels used.
same order as the spectra,
if okInd eq 1 then the array will be of length
nchan and will have a 1 for
ind used, 0 for index excluded, and npnts will
be < 0.
DESCRIPTION:
Compute the scale factors to convert from spectrometer counts to
kelvins using calon,offs. The scale factors are returned in calI.cntsToK[2]
for polA,polb
Examples:
istat=masgetfile('',bcalOn,filename=fnameOn)
istat=masgetfile('',bcalOff,filename=fnameOff)
istat=mascalonoff(bcalon,bcaloff,calI)
scale bcalOff to kelvins
TsysA=mean(bcalOff.d[*,0,*]*calI.cntsToK[0])
TsysB=mean(bcalOff.d[*,1,*]*calI.cntsToK[1])
;
; If you want to convert individual channels to degK, you must remove
; the IF band pass:
;
bpc=masmath(bcaloff,/avg)
bpc[*,0]=bpc[*,0]/mean(bpc[calI.freqInd[0]:calI.freqInd[1],0])
bcalOffK[*,0,*]= bcalOffK[*,0,*]/(
NOTES:
Current restrictions:
1. assumes calon,caloff have same integration time.
2. fixed to work with stokes data
3. need to check that polAdding works..
(See /pkg/rsi/local/libao/phil/mas/mascalonoff.pro)
NAME:
mascalonoffspc - compute cal scl factor from calon,caloff
SYNTAX: istat=mascalonoffspc(bon,boff,calI,edgeFract=edgeFract,brmsar=brmsar,$
usecalI=usecalI,mask=mask,cmpmask=cmpmask)
ARGS :
bon[n]:{}: mas struct hold cal on data.
If n > 1 then each record is assumed to have the same header
setup (nchan,freqrange,npol,etc)
the cntsToK conversion factor will be averaged over all
of these records.
boff[n]:{}: mas struct holding cal off data. boff[n] must have n = to
bon[n]
calI{} : struct holding cal value info and conversion factors.
Normally it is passed back to the caller. If useCalI
is set then the calValue info is passed in via calI.
The scale factors calI.cntsToK are always passed back
in calI.
KEYWORDS:
edgeFract[2]: float fract of the bandpass edge to ignore when averaging.
default: .06
1 entry: use both for low,high freq side
2 entery: [0] for lowFreq,[1] for high freq side
brmsar[2] : {} holds the calon,caloff rms spectra computed in masgethwcal.pro
if cmpmask=2 then you must include this variable as input.
useCalI:int if true then use calI parameter for the cal value info.
If you are doing more than one call with the same
setup then this lets you skip the mascalval() call
after the first call.
mask[nchan]: int mask to use for computing cals. !=0 --> use this channel.
Should be the same order as the spectra eg
mask[0] corresponds the bon.d[0,xx]
cmpmask : int =1 then compute the mask using calon/caloff
It will also exclude echgeFract points from each
side (def to 6%). If mask= keyword is also present
then combine mask passed in with the bad channels we found
from cmpmask.
=2 compute mask usig =1 and also fit to the rms/chan from the
for avged spectra excluding any points > 3sigma from a linear fit to the
rmsbychan vs freq
RETURNS:
istat: int 0 finished ok, -1 error
calI : {} returns calInfo and conversion factors calI.cntsToK[2]
** Structure , 9 tags, length=64, data length=64, refs=1:
CALVAL FLOAT Array[2] ; averaged cals pola,b
exposure float 0 ; integration time that cntstok corresponds to.
NUMPOL INT 2 ; number of pols 1 --> stokes I
CNTSTOK FLOAT Array[2] ; polA,b conversion cnts to Kelvins
NPNTS LONG 4096 ; if (npnts<0) then okInd was specified and
indused[nchan] has 1 for used,0 for not used
if npnts>= then indused is the indices used
for avg
FLIPPED BYTE 1 ; 1 --> spectra is flipped
EDGEFRACT FLOAT Array[2] ; edgeFract[0,1] fract to keep..spc order
USEMASK BYTE 1 ; 1 --> used mask,0--> used edgefract
INDUSED LONG Array[4096]; if okInd=0 (or not supplied) then:
indices into spc for channels used.
same order as the spectra,
if okInd eq 1 then the array will be of length
nchan and will have a 1 for
ind used, 0 for index excluded, and npnts will
be < 0.
DESCRIPTION:
Compute the scale factors to convert from spectrometer counts to
kelvins using calon,offs. The scale factors are returned in calI.cntsToK[2]
for polA,polb
Examples:
istat=masgetfile('',bcalOn,filename=fnameOn)
istat=masgetfile('',bcalOff,filename=fnameOff)
istat=mascalonoffspc(bcalon,bcaloff,calI)
scale bcalOff to kelvins
TsysA=mean(bcalOff.d[*,0,*]*calI.cntsToK[0])
TsysB=mean(bcalOff.d[*,1,*]*calI.cntsToK[1])
;
; If you want to convert individual channels to degK, you must remove
; the IF band pass:
;
bpc=masmath(bcaloff,/avg)
bpc[*,0]=bpc[*,0]/mean(bpc[calI.freqInd[0]:calI.freqInd[1],0])
bcalOffK[*,0,*]= bcalOffK[*,0,*]/(
NOTES:
Current restrictions:
1. assumes calon,caloff have same integration time.
2. fixed to work with stokes data
3. need to check that polAdding works..
(See /pkg/rsi/local/libao/phil/mas/mascalonoffspc.pro)
NAME:
mascalscl - Scale from spectrometer cnts to K using cals
SYNTAX: istat=mascalscl(bspc,calI,code,bK,bpc=bpc)
ARGS :
bspc[n]:{}: mas struct holding spectra to convert to kelvins.
calI{} : cal info structure returned from mascalonoff().
code: int how to do the band pass correction:
-1 - no bandpass correction. Use the avg cntsToK
computed by mascalonoff()
0 - use average of bspc[n] for bandpass correction.
1 - use bpc= keyword for bandpass correction.
bpc has not been normalized
2 - use bpc= keyword for bandpass correction.
bspc has already been divided by bpc.
bpc has not been normalized.
Use this for position switching:
bspc=bposOn/bposOff - 1
bpc =bposOff
3 - return only total power. No bandpass correction is
needed. bK will be an array of floats [npol,n]
KEYWORDS:
bpc{}: mas struct holding band pass correction spectra to use.
see code
RETURNS:
istat: 0 ok, -1 error
bK[n]: code=3: float[npol,n] will have averages over each bandpass
else : bk[n] masstructs scaled to kelvins.
If bpc
DESCRIPTION:
Scale spectra from spectrometer counts to Kelvins. Before calling
this routine you must call mascalonoff() to get the calI structure
loaded.
The code variable tells how to process the data.
code=3 returns only total power averaged over each spectra. This
does not need a bandpass correction.
The other codes returne mas structures with spectra scaled to kelvins.
these codes need a band pass correction. It can be bspc itself or
it can be passed in bpc= keyword.
code=2 is for position switching data. In this case, bspc is
bposOn/bposOff (the -1 is optional). You need to pass in bpc=bposOff
for the band pass correction.
code=-1 will scale the spectra to the average cntsToK value computed
by mascalonoff. The bandpass shape will remain in the returned
spectrum. Use this to scale the calOff to kelvins.
Examples:
;0. get the cal ons, cal offs:
istat=masgetfile('',bcalOn,filename=fnameOn)
istat=masgetfile('',bcalOff,filename=fnameOff)
istat=mascalonoff(bcalon,bcaloff,calI)
;1 scale the cal offs to Kelvins
istat=mascalscl(bcalOff,calI,0,bcaloffK)
warning: here you're dividing calOff,by CalOff so you
don't have any frequency dependence..:w
;2 position switching scale to kelvins:
istat=masgetfile('',bposOn,filename=fnamePosOn)
istat=masgetfile('',bposOff,filename=fnamePosOff)
bavgOn=masmath(bposOn,/avg)
bavgOff=masmath(bposOff,/avg)
bavgonoff=masmath(bposOn,bposOff,/div)
istat=mascalscl(bavgonoff,calI,2,bpc=bavgOff)
NOTES:
Current restrictions:
1. fixed to work with stokes spectra.. not yet tp
2. need to check that polAdding works..
(See /pkg/rsi/local/libao/phil/mas/mascalscl.pro)
NAME:
mascalsclspc - Scale from spectrometer cnts to K using cal spectra
SYNTAX: istat=mascalsclspc(bspc,calI,code,bK,avgK=avgK)
ARGS :
bspc[n]:{}: mas struct holding spectra to convert to kelvins.
calI{} : cal info structure returned from mascalonoffspc().
KEYWORDS:
RETURNS:
istat: : 0 ok, -1 error
bK[n]: : spectra returned in degrees K
avgK[ndump*nrecs,npol]: average spectral values in degK
It will use the mask in calI to exclude any bad channels.
DESCRIPTION:
Scale spectra from spectrometer counts to Kelvins. Before calling
this routine you must call mascalonoffspc() to get the calI structure
loaded.
There are two sets of cal routines:
mascalxxx() and mascalxxxspc().
The first set (mascalxx()) computes averages over bandpasses to generate
the conversion factors to go from spetrometer counts to degK. They are normally
used when you want to make sure that the calvalue computations do not add to
the noise of the final spectral results.
;
The second set (mascalxxxspc()) does not average over freq. It retains the
freq dependence of the cal values. This gives a better idea of tsys freq dependence
(especially for wider band mesaurements). The drawback is that the cal computation has the
potential to increase Tsys (if the cal integration is much shorter than the
normal data acquistion). These routines can be used with a hardware winking cal where
1/2 the time the cal is on, and 1/2 the time the cal is off. They used the cal deflection
for the band pass correction. If there is any frequency dependence in the signal before
the cal insertion (eg.. standing waves..) they will not be cancelled.
Examples:
;0. get the winking cal data
istat=masgethwcal(fnmI[ifile],calIhw,bonar,boffar,bavgar=bavgar,brmsar=brmsar,$
/avgrow,/dctinp
; get the cal values for calib files
istat=mascalvalspc(bavgar[0].h,calI)
; compute cal scale factors using avg spectra
; exclude 10% edges for sb, 6% for xb
; cmp mask using calration and the rms/chan from brmsar
cmpmask=2
istat=mascalonoffspc(bavgar[0],bavgar[1],calI,edge=edge,cmpmask=cmpmask,brmsar=brmsar)
; now scale calOff avg spectra to degK .. (or you could do the 40ms caloff spectra
istat=mascalsclspc(bavgar[1],calI,avgK=avgK)
NOTES:
(See /pkg/rsi/local/libao/phil/mas/mascalsclspc.pro)
NAME:
mascalval - return the pol A/B cal values/info for a mas struct
SYNTAX: stat=mascalval(hdr,calI,date=date,edgeFract=edgeFract,
mask=mask,okind=okind)
ARGS:
hdr: {b.hdr} header for spectra input via masget
KEYWORDS:
date[2]: intarray [year,dayNum] if provided, then compute the calvalues
that were valid at this epoch. Default is to use date in hdr
edgeFract[2]:float fraction of the band to discard at the edges.
edgeFract= single .. fraction to discard on both sides.
edgeFract=[2] .. edge[0] -fraction on lowFreqedge
edge[1] -fraction on hiFreqedge.
Default value=.06
Note that hi/low is in frequency (even if the
band was flipped).
mask[nchan]: int if supplied then mask saying which frequencies
to use when doing the average.
1=use, 0= don't use this freq channel.
nchan must match channels in spectra.
The freq order of freqMask is the same as the spectra.
If the spectra are flipped (hdr.cdelt1<0) then
freqMask[0] is the highest freq,instead of the lowest.
The same mask is used for polA and polB
okInd : int is 0 then
calI.npnts >=0 is good points to use for total power
calI.indused[npnts] has indices of freq chan use has
!= 0 calI.npnts <0 is good points to use for total power
calI.indUsed[nchan] with 1=used bin, 0=don't use bin
RETURNS:
calI:{} structure containing:
calI.calVal[2] :polA, b cal value averaged over the band
calI.numPol 1 :if only stokes I data returned. The
cal values will be averaged and returned in
both calval[0] and calVal[1].
calI.exposure: zeroed. filled in by mascalonoff
calI.cntsToK[2].. zeroed. filled in by mascalonoff
calI.npnts 4000: >0 points used for averages
calI.npnts 4000: <0 -points used for averages
and indused[nchan] has 0,1 for bins to use.
calI.freqRange[2]: freq Mhz used to compute the average
calI.flipped : 1 if freqRange[0]>freqRange[1]
calI.edgeFract[2]: fraction to drop at each edge when
computing spectra. Ignored if usemask
freq order same as spectra
calI.usemask : if 1 then mask was supplied to generate
indused[]
calI.indUsed[m] : calI.npnts> 0 m=npnts indices into freq array
calI.indUsed[m] : calI.npnts< 0 m=nchan with1=use,0=don't use
stat: int .. -1 error, 0 got the values ok.
DESCRIPTION:
Return the cal values in degrees K for the requested spectra. The
calvalues will be averaged over the frequency band excluding
edgefract fraction on each edge. If mask= is included, then
only the channels with mask[ichan]=1 will be used in the average.
By default the date in the header is used to determine
which cal values should be used. You can override this with the date
keyword.
The data is returned in the calI structure.
The following is returned with calI.freqInd[0] < calI.freqInd[1] so
you can say mean(b.d[calI.freqInd[0]:calI.freqInd[1]). This will
make freqRange[0]>freqRange[1] if flipped is set. Returned values
with this order are:
calI.freqInd[]
calI.freqRange[] .. if flipped set then freqRange[0] > freqRange[1]
calI.edgeFract[] .. These will now correspond to freqInd
It will be the reverse of the input if flipped is 1.
calI.flipped 1 if freqRange[0]>freqRange[1]
calI.indUsed[nchanUsed] index into freq for channels used.
EXAMPLE:
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 average is performed over the available range.
SEE ALSO:
mascalscl, gen/calget gen/calval.pro, gen/calinpdata.pro
(See /pkg/rsi/local/libao/phil/mas/mascalval.pro)
NAME:
mascalvalspc - return the pol A/B cal values/info for a mas struct over the freq
SYNTAX: stat=mascalvalspc(hdr,calI,date=date)
ARGS:
hdr: {b.hdr} header for spectra input via masget
KEYWORDS:
date[2]: intarray [year,dayNum] if provided, then compute the calvalues
that were valid at this epoch. Default is to use date in hdr
RETURNS:
calI:{} structure containing:
calI.nchan number of freq chan
calI.freq[nchan] freq (mhz) for the data
calI.npol number of pols. if 4 pols, we just return a,b
calI.flippped if true then the hdr channels were flippped in freq.
(the cal values will be in the same order).
calI.calVal[nchan,npol]:polA, b cal value for the freq range of spectra (in degK)
stat: int .. -1 error, 0 got the values ok.
DESCRIPTION:
Return the cal values in degrees K for the requested spectra. The
measured calvalues will be interpolated to the freq range of the specra.
If the requested freq range is outside the cal values, then the cal values
at the edge will be repeated.
By default the date in the header is used to determine
which cal values should be used. You can override this with the date
keyword.
This routine can be used in place of mascalval() if you want the cal values
across the freq band. This routine is normally called by mascalonoffspc().
NOTE:
mascalonoffspc,mascalsclspc, gen/calget1
(See /pkg/rsi/local/libao/phil/mas/mascalvalspc.pro)
NAME:
masclose - close a mas file for i/o
SYNTAX: masclose,desc,all=all
ARGS:
desc: {masdescr} - descriptor to close (returned by masopen)
KEYWORDS:
all: if set then close all open descriptors.
DESCRIPTION:
Files opened with masopen() need to be closed with masclose() so that
the resources are freed up.
EXAMPLE:
filename='/share/pdata/pdev/phil/071106//testfits.20071107.b0s0.00000.fits'
istat=masopen(filename,desc)
.. process the data in the file
masclose,desc .. this closes the file when done with the processing.
(See /pkg/rsi/local/libao/phil/mas/masclose.pro)
NAME:
mascmppwrrms - compute total power exluding rfi via rms.
SYNTAX: istat=mascmppwrrms(bar,tpI,fixblank=fixblank,norm=norm,startmask=startmask,$
maskar=maskar,brms=brms,$
)
ARGS:
bar[N]:{} array of data from masgetm
KEYWORDS:
fixblank: if set the correct for any blanked spectra so they
all have the standard integration time.
norm : if set then normalize tp to median value
startMask[nchan] : if supplied then this is the mask to begin with.
we then add to this. 0-> skip,1-> used
These points will always be excluded
RETURNS:
istat: int 0 ok, -1 error
tpi:{} structure hold total power info
maskAr[nchan,npol]: float mask used
brms : return rms by chan if supplied
(See /pkg/rsi/local/libao/phil/mas/mascmppwrrms.pro)
NAME:
mascrossfit2d - 2d gauss fit to az,el strips done on 12meter.
SYNTAX: nfit=mascrossfit2d(cIAr,fitI,verb=verb,remMod=remMod,pol=pol,$
zfitar=zfitar,$
wait=wait,dbg=dbg,nosidelobe=nosidelobe,rot=rot,rembase=rembase)
ARGS:
cIAr[ncross]:{} array of structs hold strip info for 2 strips in each cross
see mascrossinp()
KEYWORDS:
verb : if set then plot results each fit
remMod : if set then remove current model before fitting.(not yet implemented)
pol : int 1->polA,2=polB.. default add pols
wait : if set and verb set, then wait for keystroke before moving to
the next fit.
nosidelobe: int if set then don't exclude sidelobes from the fit. use this for weak sources.
def is to loop fitting twice, where the 2nd loop excludes the sidelobe regions.
rot : int if set then include rotation of beam ellipse in the fit. by default don't rotate.
the 12m beam is not elliptical, rotating just causes problems.
(The 305m had an elliptical beam).
Warning: if /rot used then returned results (widths,pnterr are along
the rotated major,minor axis of the ellipse (not the az,el directions).
rembase: int if set then remove a linear baseline from the az and el strips before the
2d gaussfit is done. The fit coef[0] will no longer contain tsys (you'll need
to look at the baseline fit info.
RETURNS:
nfit : int number masmon2dfits,
fitI[nfit]: {} array of structs holding fitinfo for each cross.
zfitAr :fltarr(pntsStrip,2,nbms,nfit) .. holds evaluated fit. If /rembase used then the
baselines will have been removed from zfitAr
DESCRIPTION:
Perform a 2d gaussfit to the az,el strip pairs for each cross taken on the 12meter.
The cross is done with constant rate great circle offsets relative to the
source we are tracking.
Strip 1: az strip. constant great circle rate in az, while el offset (relative to sourcer)
doesn't change
Strip 2: el strip. constant great circle rate in el, while az offset relative to source
does not change.
Note: the fit offsets are relative to the source. Since the source is moving in az,el
the encoder el will change during the az strip.
-----------------------------------
Some of the keyword options:
pol=1,2 (default add pols)
Normally the routine averages polA,B prior to fitting.
If the pol= keyword is used then it ;will fit either polA(1) or polB(2).
nosidelobe: (def:loop twice excluding sidelobes locations)
By default the routine iterates the fit 2 times
1stloop: 2d gaussfit everything
2ndLoop: from first fit, find the beam center and width. Then compute where the sidelobes
are located. On the 2nd loop exclude these sidelobes from the fitting.
This is important for strong sources casA,cygA, crab.
For weak sources the sidelobes are not detectable. you should set /nosidelobe.
This will give you a better fit by including the area of the sidelobes.
rot: (default no rotation).
if /rot is set then the fit will include rotating the x,y axis to align with the
major,minor axis of the elliptical beam (this was needed for the old 305m elliptical beam).
The bmwidths, pnterrs will then be relative to the ellipse axis rather than the az,el
directions.
The 12meter beam is circular, so the rotation angle is not well defined. Including /rot
will cause the bmwidths and pnterrors to point in random directions (unless the beam
ends up becoming elliptical?..)
---------------------------------
The fitI structure contains:
help,fiti1,/st
SRCNM STRING '3C286' ; srcname
STRIPLEND FLOAT 1.00000 ; strip length in deg
PNTSTRIP INT 600 ; samples in strip
PATID LONG 221500426 ; pattern id for cross
TMST LONG64 Array[2] ; starttm (secs 1970) az,el strips
OFFARD DOUBLE Array[600, 2] ; strip offset from source (deg) for az,el used for fitting
ELOFF DOUBLE Array[600, 2] ; elevation offset from middle of cross. used for elSlope in fit.
AZREQCEND DOUBLE 299.98128 ; central az for cross
ELREQCEND DOUBLE 60.868295 ; central el for cross
NCOEF INT 8 ; number of coefin fit
NGOODFITS LONG 7 ; number of good fits(in the Nbeams)
REMOVEBASELINE INT 1 ; not 0 if we removed a linear baseline first from strips
FITIBM STRUCT -> Array[7]; fitinfofor each beam used
---------------------------------
info from each frequency band fit (up to 7 bands (also called bms.. alfa holdover:)
help,fiti.fitibm,/st
Z DOUBLE Array[600, 2] ; z values used for az,el fits
; if baseline removes,then it is nolonger in z
CFRMHZ DOUBLE 8219.0000 ; center freq MHz for this band
COEF DOUBLE Array[8] ; coef from fit (see below)
BLFIT STRUCT -> Array[1]; blfit baseline fitinfo (see below(
SIGCOEF DOUBLE Array[8] ; sigma for each coef from the fit
SIGFIT DOUBLE 0.032022892 ; sigma from the fit
CHISQ DOUBLE 0.0010306230 ; chisq from the fit (but we don't give correct errors)
TROUBLE INT 0 ; if not 0 then fit didn't converge
---------------------------------
if /rembase is used then a linear baseline will be removed from the az,el strips before the
the 2d gaussfit (so the 2d fit coef[0] will be close to 0). The linear fits info
is stored in .blfit.az,el . You can extract the Tsys from thes structs.
help,fiti.fitibm.blfit,/st
AZ STRUCT -> Array[1] az baseline fit info (if used)
EL STRUCT -> Array[1] el baseline fit info (if used)
help,fiti1.fitibm.blfit.az,/st
COEF DOUBLE Array[2] bl=coef[0] + offarD[*,x]*coef[1]
AVGVAL DOUBLE 99.322503 avg value of fit
and similar for ..blfit.el..
You can reconstruct the baselinf from:j
azbaseline=coef[0] +fitI[i].offarD[*,0]*coef[1]
elbaseline=coef[0] +fitI[i].offarD[*,1]*ceof[1]
The average baseline fit values (Tsys) are stored int
blfit.az.avgVal
blfit.ez.avgVal
To get tsys you should probably average these to values.
---------------------------------
The 2dgaus fit is:
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2 - yp^2/sigyp^2] + a7*el
xp,yp have been rotated to align the x,y axis with the major,
minor axis of the possibly ellipitcal beam.
note: rotation is only done if /rot is supplied. The 12m data shouldn't be rotated
(the beam is not elliptical)
where:
eg:
xm=(x-x0)
ym=(y-y0)
if /rot was supplied then x,y axis are rotated to the major, minor axis of the
elliptical bm
if /rot not supplied then the rotation is not done (th=0._
xp = xm*cos(th) + ym*sin(th)
yp =-xm*sin(th) + ym*cos(th)
The fit parameters used in the fit equation are:
a0 - constant
a1 : amplitude
a2 : xoffset ,az coord direction offset is az units
a3 : yoffset ,el coord direction offset in el units
a4 : sigx^2 ,in prime coordinate system. See returned units
a5 : sigy^2 ,in prime coordinate system: see returned units
a6 : theta ,rotate unprimed to primed aligned along ellipsoid of beam,Deg
a7 : elslope ,The el slope in amplitude units per el unit
RETURNED UNITS:
fitI[i].coef[0] = offset (or tsys) in input units. if /rembase is used then this will
be close to 0 (baselines removed before 2dfit).Look at blfit info for Tsys.
fitI[i].coef[1] = Amplitude of source deflection from fit.
fitI[i].coef[2] = az error in az units. Note that since the strips
are done az greatcircle strips, the az error in a great
circle error.
fitI[i].coef[3] = el error in el units.
fitI[i].coef[4] = az fwhm from fit (we convert to fwhm from sigmas when
returning the fit info.
fitI[i].coef[5] = el fwhm from fit
fitI[i].coef[6] = theta in deg. rotation angle to get az,el axis to
align with gaussian axis. will be 0 unless /rot is supplied.
fitI[i].coef[7] = slope of baseline in el direction
dz/del
You should check fiti[i].trouble. if it is not zero, you should ignore the fit.
Note:
the errors are computed at:
(x-x0) (y -y0)
: assume the strip is -10 to +10. No error would be 0
: assume the peaks is at +3 in x (after the center) . then x0 =3
: (x-x0) would put the peak back at zero offset.
: so add these offsets to the current model to have the center of the strip
be at the peak.
:-> Need to check this..
30aug21: added patid to struct
09sep22:
1.used weights to exclude sidelobes
- fit twice
- 1. first equal weights
- 2.find center,bmwd use to set weights for sidelobe to 0.
2. when doing robfit_poly for baseline, tsys, i was using mask to
not include main beam and sidelobes. but fit routine only uses mask
on first iteration, after that i just keeps points within n sigma.
- change to use 1st 1/6, last 1/6 of strip
1 bm/ 4bms sidelobe,mainbm,sidelobe/ 1bm
10mar23:added rot= keyword. only rotate if /rot is set. New default is no rotation
14mar23:--> WARNINGS. the returned struct fitI has changed. It now contains
fitI.removebaseline and
fiti.fitibm.blcoef(2,2) ({a0,a1},{az,el}) for linear baseline fits using fit
(See /pkg/rsi/local/libao/phil/mas/mascrossfit2d.pro)
NAME:
mascrossinp - input cross data, compute totpwr
SYNTAX: istat=mascrossinp(cFI,cI,baz=baz,bel=bel,brms=brms,edge=edge,$
inpmask=inpmask,$
freqExcl=freqExcl,maskUsed=maskUsed)
ARGS:
cFI : {} cross file info (see masfindcross)
KEYWORDS:
edge[2,nbms] : float fraction of edge to exclude on each side when computing total power.
def .06 (6%)
inpmask[nchan]:intarr if supplied then use this as the mask (don't bother computing one)
it will still include edge and any freqexcl[] values. Use this if you precompute
a mask using many input scans.(see masmakemask)
freqExcl[2,m]:float MHz. if supplied then don't include these freq in total power computation
You can specify up to m freq sections
freqExcl[0,..] is min Freq,freqExcl[1,..] max freq for this section
For each beam, it will exclude all sections within the current band.
RETURNS:
istat : 1 ok
-1 error
cI : {} cross info returned here.
maskUsed: intarr(nchan,nbms) If supplied then return the mask of good channels
we computed. 1-> used, 0-> not used
The same mask is used for polA and polB
-->NOTE:
DESCRIPTION:
The routine will input the spectral data files for a cross pattern and then compute
the total power. It will process a single cross on each call.
The user inputs:
- a cFI (CrossFileInfo structure) for the cross to process (see masfindcross() routine).
- edge= keyword
-can be used to exclude a fraction of the spectra on each edge of the spectra.
- default value is .06
- The fraction to exclude can be 0 to .99 .
- The edge dimension can be:
edge=.06 .. exclude 6% of spectra on each band edge for all bands
edge=[.1,.06] exclude 10% on the left edge, 6% on right edge for all bands
edge=[2,bands] specify the left,right edge for all bands
- freqExcl= freqExcl[2,m] freq range sections to exclude in total power computation. Tested
separately for each beam.
The routine will then:
For each bm:
1. input the azimuth and el strips
2. interpote across the DC freq channel
3. using the 1st 1/3 and last 1/3 of the azimuth strip, compute rms/mean by channel
- we don't use the el strip since tsys can change vs el.
4. compute a mask for freq channels to use in the total power computation.
- intialize the mask using the edge= and freqExcl= keywords.
- do a robust linear fit to the rms/mean by channel
- exclude any points > 3 sigma from the fit.
- use the same mask for polA,B and az,el strips
5 if cals are present compute the counts to Kelvins scale factor using the same mask as the
strips
6. compute the total power for the az, el strips for polA and polB using the above mask.
and convert to kelvins if cals are present
Return the cI structure
If the following keywords are present then return:
let n be the number of rows in the fits file for a strip
baz=baz .. return the n mas structures holding the individual az strip spectra
bel=bel .. return the n mas struxtures holding the individual el strip spectra
brms=brms. return the 1 mas brms structure computed from the az strip
maskUsed= return maskUsed[nchan,nbms] int array. If a channel=1 it was included, 0
means it was not included in the total power measurement.
EXAMPLE:
; get the file list
proj='x101'
yyyymmdd=20220406l
; find the files
ncross=masfindcross(yyyymmdd,proj,cFI,/verb)
; loop processing the crosses
; you could take a subet by source or rcvr , or just do everything
; eg ii=where(cFI.srcName='3C144',ncross)
; you also might want to keep track of the freq channels excluded
for icross=0,ncross-1 do begin
istat=mascrossinp(cFI[icross],cI,maskused=maskUsed)
if icross eq 0 then begin
cIAr=replicate(cI,ncross)
nbms=cI.nbms
nchan=cI.bmI[0].h.chan_in ; freq chan used
maskAr=fltarr(nchan,nbms)
endif
cIAr[icross]=cI
maskAr+=maskUsed
endfor
maskAr/=ncross ; to get a fraction 0..1
;
; plot the the total power az,el strips for a cross
;
plot,cIar[0].bmI[0].tpAz[*,0] ; polA
oplot,cIar[0].bmI[0].tpAz[*,1],col=colph[2] ; polB
----------
The cI structure:
the structure shown bel1ow used 7 freq bands (bms) for each cross:
help,cI,/st
.. info common to all beams
RCV STRING 'xb' : receiver used xb or sb
SRCNM STRING '3C461' : source name
AZLEND DOUBLE 1.2000000 ; azimuth strip length deg
ELLEND DOUBLE 1.2000000 : elevation strip length deg
CROSSCEND DOUBLE Array[2] ; az,el center for cross
STRIPCEND DOUBLE Array[2, 2] ; [i,j] i= az,el strip center, j=azstrip,elstrip
AZPOSD DOUBLE Array[1200, 2] ; azposition [deg] for each sample az,el strips
ELPOSD DOUBLE Array[1200, 2] ; elposition [deg] for each sample az,el strips
AZOFFD DOUBLE Array[1200] ; az offset azimuth strip .great circle
ELOFFD DOUBLE Array[1200] ; el offset el strip .great circle
SMPSTRIP LONG 1200 ; number of samples in each strip
SMPTM DOUBLE 0.10000000 ; sample time (secs)
MAXTRKERRD DOUBLE Array[2] ; maxEncoder tracking error [deg] az,el strip
NBMS INT 7 ; number of beams this cross.
units int 1,2,3 ; 1 - spectrometer counts, 2=calunits, 3=degK
.. beam specific info
BMI STRUCT Array[7] ; beam specific info (see bel1ow)
The beam specific info will depend on whether cals were present
with no cals:
help,cI.bmI,/st
BM INT 0 : beam number 0..6
CFR DOUBLE 8221.0000 ; rf center frequency MHz
HDRAZ STRUCT -> MASFHDR Array[1]: first header for azimuth strip
HDREL STRUCT -> MASFHDR Array[1]: first header for el strip
FPNTS FLOAT Array[2] : fraction of freq chan used for tp: polA,B
FLIST STRING Array[2] : filenames for az,el strips
TPAZ DOUBLE Array[1200, 2] : computed total power azStrip:polA,B
TPEL DOUBLE Array[1200, 2] : computed total power elStrip:polA,B
with cals present
help,cI.bmI,/st
BM INT 0 : beam number 0..6
CFR DOUBLE 8221.0000 ; rf center frequency MHz
HDRAZ STRUCT -> MASFHDR Array[1]: first header for azimuth strip
HDREL STRUCT -> MASFHDR Array[1]: first header for el strip
FPNTS FLOAT Array[2] : fraction of freq chan used for tp: polA,B
FLIST STRING Array[2] : filenames for az,el strips
FLISTCAL STRING Array[2] : filenames for calon,off
CALI STRUCT -> Array[1] ; cal info
TPAZ DOUBLE Array[1200, 2] : computed total power azStrip:polA,B
TPEL DOUBLE Array[1200, 2] : computed total power elStrip:polA,B
and the calI struct will have:
IDL> help,ci.bmi.cali,/st
CALVAL FLOAT Array[2] : cal value in degK (pola,b). if =1 then not yet measured
NUMPOL LONG 2 : number of pols
EXPOSURE FLOAT 0.100000 ; sample time for spectra
CNTSTOK FLOAT Array[2] ; convert spectrometer counts to degK
NPNTS LONG 900 ; nfreq points used
FLIPPED BYTE 0 ; was the spectra flipped in freq
EDGEFRACT FLOAT Array[2] ; edge fract excluded each side. same for polA,B
USEMASK BYTE 1 ; 0- just used mask, 1 - used mask
MASK BYTE Array[1024] ; resultant mask 1=chanused, 0 not used. includes mask and edgefract
Notes on Positions
The az,el positions come from the fits header h.azimuth, h.elevatio. These values have been interpolated to crval5 of the hdr.
These positions do not have the model included (so back computing should give ra,dec without needing the model offsets.
Finally the reported values have been interpolated to the center of each spectral/totalpower sample.
(See /pkg/rsi/local/libao/phil/mas/mascrossinp.pro)
NAME:
masdpsfind - find scans for dps patterns
SYNTAX: n=masdpsfind(projId,
patI,npat,upatIdAr,nsrc,usrcAr,$
yymmdd=yymmdd,appbm=appbm,dirI=dirI,bm=bm,band=band,grp=grp)
ARGS:
projid: string project id to search for
KEYWORDS:
yymmdd[m]:lonarr yyyymmdd limit to any of the dates in this array
srcNm[l] :strarr limit to source names in this array. Must match
object name in fits header
appbm: apply bm number to each directory. This should normally
be set:
dirI[2]: string if data not in default /share/pdataN/pdev
see masfilelist for usage.
bm : int if supplied limit to this beam
band : int 0,1 limit to this band. Note if you have single pixel
data you probably should set band=1 or masfilelist may
not find the files.
grp : int limit to this group 0, or 1.
RETURNS:
n :int number proccessed patterns. Each pointing pattern
can generate multiple processed patterns if you use
multple bms.
patI[n]: info on each pattern
npat : int number of pointing patterns executed.
nsrc : int number of unique sources
srcAr[nsrc]: string list of sources used.
patI STRUCTURE:
- there will be one entry for each pattern and bm used
If you have 4 bms then a single patId will have 4 entries
in patI[]
patI[i].srcNm - source name
patI[i].nscans - # of scans in pattern ..
patI[i].bm - bm number 0..6
patI[i].grp - grp 0 or 1
patI[i].flist[maxScans] - list of filenames:
same orde as scanTypeAr
patI[i].sumI[maxScans] - summary info for each scan
DESCRIPTION:
Find all of the dps patterns for the given project and constraints.
The routine accepts the same parameters at masfilelist:
yymmdd=,dirI=,bm=,band=,grp=num=,appbm=appbm
(except that yymmdd= can be and array of dates)
The routine will pass back an array of patI structs holding info on
each pattern.
There is a distinction between a pointing pattern and a processed
dataset. You can have multiple processed datasets for each pointing
pattern if you used multiple beams or groups.
EXAMPLE:
1. Find the dps info for project a2516 on 20100930 group 0.
n=masdpsfind('a2516',dpsI,npat,upatIdAr,nsrc,usrcAr,$
yymmdd=20100930,grp=0,band=1,/appbm)
dpsI[n] - holds the info
npat - number of pointing patterns found
upatIdAr[npat] - unique pattern id for each pointing pattern.
This can include multiple dpsI[] entries if multiple
bms taken on each pointing pattern.
nsrc - Number of unique on source names found
usrcAr[nsrc] - list of source names.
2. now process all of the patterns found using group 0 bm 1.
ii=where((dpsI.grp eq 0 ) and (dpsI.bm eq 0),cnt)
for i=0,cnt-1 do begin
j=ii[i] ; index in dpsI for this dataset
istat=masdpsp(dpsI[j].flist,b)
if i eq 0 then b0ar=replicate(b[0],cnt); generate large array
b0ar[i]=b
endfor
this assumes that all of the bm=0 data sets have the same
number of channels (since we are trying to store them all in an array).
3. plot out the first set of results
masplot,b0ar[0]
(See /pkg/rsi/local/libao/phil/mas/masdpsfind.pro)
NAME:
masdpsp - double position switching
SYNTAX: istat=masdpsp,flist,b,fdir=fdir,flux=flux,fnmI=fnmI,
srcOn=srcOn,srcOff=srcOff,$
bpon=bpon,bpoff=bpoff,src_onoff=src_onoff
ARGS:
flist[4]: string srcOn,srcOff,bpOn,bpOff filenames
KEYWORDS:
fdir : string If supplied then prepend this to each entry in flist.
You could also just include the directory in flist.
flux : float bandpass calibrator flux
fnmI[4]: {} srcOn,srcOff,bpOn,bpOff from masfilelist. If supplied use
this rather than flist[].
RETURNS :
istat : 1 ok, -1 error reading a file
b : {masget} (srcOn-srcOff)/(bpOn-bpOff) * flux averaged.
srcOn[n] :{masget} the srcOn records
srcOff[n]:{masget} the srcOff records
bpOn[m]:{masget} the bandpass on records
bpOff[m]:{masget} the bandpass off records
src_onoff:{masget} srcOn/srcOff averaged
DESCRIPTION:
Perform double position switching processing:
b=avg(srcOn-srcOff)/avg(bpOn-bpOff) * fluxBpCalibrator
The routine will also return the individual records for srcOn,srcOff,
bpon,bpoff if you supply the keywords as well as the averaged
src_onoff = avg(srcOn/srcOff)
SEE ALSO
masdpsfind()
(See /pkg/rsi/local/libao/phil/mas/masdpsp.pro)
NAME:
masdpsproc - process multiple dps patterns
SYNTAX:npat=masdpsproc(ProjId,nband,bar,patI,srcNm=srcNm,cfr=cfr,dateAr=datear,$
epscfr=epscfr,exclpatid=exclpatid,wait=wait,$
barp=barp)
ARGS:
projId: string project id to match.
KEYWORDS:
srcNm[m]:strarr limit to sources in srcNm array. The names
should match those in the fits.object header.case matters.
cfr :float only include patterns whose average freq is cfr(in Mhz).
The routine averages all of cfr's of a pattern.
also see epscfr=
dateAr[l]:lonarr only include patterns whose file dates match a
date in datear[]. The format is yyyymmdd.
epscfr :float Include patterns with average cfr within epsCfr of
the cfr keyword. Units are MHz. The default value
is 10. Mhz.
exclPatId[n]: int exclude any pattern id's in this array.
wait : int if true then wait for a keyboard entry between
each plot.
RETURN VALUES:
npat:int number of patterns we found
nband:int number of frequency bands in each pattern
bar[nband,npat]: {} array of mas struct holding dps processed patterns
patI[nband,npat]: {} array of structs holding info on each pattern.
barP[nband,npat]: {} array of mas struct holding dps processed patterns
after polarizations have been added.
TERMINOLOGY:
spectra: a single integrated spectrum. Typically 1 second.
scan : a set of spectra with the same receiver setup and telescope motion,
taken for a requested number of spectra..
eg. the on source, or off source are separate scans. The header holds
a h.scan_id to identify each scan
pattern: A number of scans that are taken to form a particular pattern. In
dps there are at least 4 scans in a pattern: calibrator on,off and
source on,off. The header contains hdr.pattern_id
band : a single scan,pattern can have up to 14 frequency bands (corresponding
to the 14 mock boxes.
DESCRIPTION:
Find and process a number of dps (double position switch) patterns.
This routine uses masdpsfind and masdpsp.
The user specifies the project id to search for. The caller can further limit
the patterns used by:
srcNm[m] : this is an array of source names. For a pattern to be included
its hdr.object must match one of these names. Case in important.
dateAr[l]: a list of dates that each pattern must match. Format is
yyyymmdd. The files includes must have one of these dates.
cfr,epscfr: Mhz. Average the cfr of all of the bands in a pattern. To be
included, it must be withing epsCfr(mhz) of the supplied cfr.
exclPatId[n]: exclude any pattern id that is in this array. You can use
this to skip bad patterns.
After finding all of the patterns that meet the requirements, the routine
will process one pattern at a time. Each processed pattern will be plotted to
the screen. If the wait keyword is set, then the user will have to hit a key
to continue after each plot.
The routine returns the processed pattern data in bar[nbands,npat]. If
barp=barp keyword is present, then then barp[nbands,npat] variable will contain
the polarization averaged spectra.
---------------------------------------------------------------------
EXAMPLE:
process the a2772 5-6 ghz data taken 07jun13 -> 12jun13.
projId ='a2765'
scrNm ='ARP220'
cfr =5500.
yymmddar=[20130607,$
20130610,$
20130611,$
20130612]
npat=masdpsproc,projId,nband,bar,patI,srcNm=srcNm,cfr=cfr,dateAr=dateAr,$
barp=barp
The returned values are:
npat=20
nband=14
bar[14,20]
barp[14,20]
patI[14,20]
The patI struct contains
IDL> help,pati
PATI STRUCT = -> Array[14, 20]
IDL> help,pati,/st
** Structure <8e1d68>, 7 tags, length=7664, data length=7568, refs=2:
SRCNM STRING 'ARP220'
BM INT 0
GRP INT 0
NSCANS INT 4
PATID LONG 315800057
FILELIST STRING Array[6]
SUMI STRUCT -> Array[6]
The SumI holds summary info for a scan. It includes the fits header
IDL> help,pati.sumi,/st
** Structure <91e808>, 13 tags, length=1256, data length=1241, refs=2:
OK INT 1
FNAME STRING '/share/pdata1/pdev/a2765.20130607.b0s1g0.00200.fits'
FSIZE LONG64 20643840
NROWS LONG 300
DUMPROW LONG 1
DATE LONG 20130607
BM INT 0
BAND INT 1
GRP INT 0
NUM LONG 200
H STRUCT -> MASFHDR Array[1] <--- standard fits header
HMAIN STRUCT -> PDEV_HDRPDEV Array[1]
HSP1 STRUCT -> PDEV_HDRSP1 Array[1]
SEE ALSO:
masdpsfind(), masdpsp()
(See /pkg/rsi/local/libao/phil/mas/masdpsproc.pro)
NAME:
masfilelist - get list of mas files
SYNTAX: nfiles=masfilelist(dirBase,fnmiar,projId=projId,year=year,mon=mon,$
day=day,yymmdd=yymmdd,bm=bm,band=band,grp=grp,num=num,$
appBm=appBm,pdev=pdev,psrfits=psrfits,dirI=dirI,verb=verb)
ARGS:
dirBase: char base directory to start search
if '' then use default /share/pdata
KEYWORDS:
projId : string beginning of filename.if null then any
You can also use a reg exp as eg:
proj='x101_7_[12]'
year : long year. < 0==> any year
mon : long mon . < 0--> any month
day : long day . < 0--> any day
yymmdd : lond 080615.. date for data. If supplied, ignor year,mon,day
format is yymmdd or yyyymmdd. If you want to search
and entire month yymm00 or yyyymm00 will also work.
bm : int beam to use 0..6. if not present or < 0 use all beams
band : int band to get, 0 or 1 .. not present or <0 --> all
grp : int grp to get 0 or 1. .. not present or < 0 --> check
both g0,g1. grp 0 files on psrv1-7, grp1 on 8-14
num : long sequence number. If not present or < 0 get all
appBm: if set then append (bm+1) to the end of dirBase when
defining the initial base directories
pdev : if set then look for .pdev files rather than .fits
psrfits: if set then assume filename contains srcname
dirI : strarr(2): directory info. If the files are in a non standard place
then pass in location in dirI
dirI[0]=prefix for dir: /share/pdata
Num 1..n goes here
dirI[1]=postfix for dir: /pdev
RETURNS:
ntot : long number of .fits files found
fnmIAr[ntot]:{} an array of structures holding info on the
mas files found
DESCRIPTION:
Search for mas fits files or psrfits files (with the /psrfits keyword).
The search starts in /share/pdataN by default. You specify the base directory and the
filename characteristics to search for (via the keywords).
Example:
The spectrometer files are nfs mounted to
/share/pdataN/pdev/
This would require dirBase="/share" which would make the search to
broad. The appBm keyword can narrow this down. If set, it will append
all of the beam numbers +1 to to the end of dirBase before searching..
so dirBase=/share/pdata
and the searches would go thru
/share/pdata1, /share/pdata2, .. etc..
The returned structure contains:
nmAr='/share/pdata/mas/agc110443/x107.20070123.agc110443.b0a.00000.mas'
istat=masfilelist(nmAr,fnmI)
help,fnmI,/st
* Structure masfnmInfo 8 tags, length=60, data length=60:
DIR STRING '/share/pdata1/pdev/'
FNAME STRING 'philtest.20080219.b0s0g0.00000.fits'
PROJ STRING 'philtest'
DATE LONG 20080219
src string '' if psrfits file
bm int 0
band int 0
grp int 0
num long 0
(See /pkg/rsi/local/libao/phil/mas/masfilelist.pro)
NAME:
masfilesum - get summary info for files
SYNTAX: nsum=masfilesum(flist,sumI,desc=desc,fnmI=fnmI,list=list,patid=patid)
ARGS:
flist[n]:strarr filenames for summary. Includes directory.
KEYWORDS:
desc :{} descriptor returned from masopen. If supplied then ignore flist
and use desc as the file to look at.
fnmI[n]:{} structure returned by masfilelist. If this keyword is
present then ignore flist and use fnmI for the files
to use.
list : if set then output oneline summary to stdout of file
patid : if set then print pattern id rathern than scan id with list
RETURNS:
nsum: long number of summaries we return. Don't count files that had trouble being read
sumI[n]: {} summary info for mas files.
DESCRIPTION:
Return summary info for a list of files. You can specify the files with
the argument flist or the keyword fnmI.
The returned summary info contiains:
help,sumI,/st
OK INT 1
FNAME STRING '/share/pdata1/pdev/x108.20080530.b0s0g0.00000.fits'
FSIZE LONG64 84389676 .. bytes in file
NROWS LONG 4 .. rows in fits file
DUMPROW LONG 319 .. dumps in first row
DATE LONG 20080530 .. yyyymmdd
BM INT 0 .. 0..6 for alfa
BAND INT 0 .. 0,1 low,hi band (at if)
GRP INT 0 .. grp 0,1
NUM LONG 0 .. file number .nnnnnn.fits
H STRUCT 1st row of fits file (without data,stat)
HMAIN STRUCT pdev main header
HSP1 STRUCT pdev sp1 header
The fits file row header contains all the info for the first row
(except for the spectral and status data).
EXAMPLE:
1. get 1 file
fname='/share/pdata1/pdev/x108.20080530.b0s0g0.00000.fits'
nsum=masfilesum(fname,sumI)
2. Get all the x108 files under /share/pdata1 (online files) for 30may08
n=masfilelist(junk,fnmI,projid="x108",year=2008,mon=5,day=30,/appbm)
nsum=masfilesum(junk,sumI,fnmI=fnmI)
742 files returned..
NOTE:
nsum will always equal the number of elements in flist or fnmI. You need
to check sumI.ok eq 1 to see if the data is ok. This lets you keep track
of bad files (say with 0 rows).
(See /pkg/rsi/local/libao/phil/mas/masfilesum.pro)
NAME:
masfindcross - find cross files
SYNTAX: ncross=masfindcross(yyyymmdd,proj,cfI,rcv=rcv,verb=verb,dirI=dirI)
ARGS:
yyyymmdd: long date to search for
proj : string project name for files
KEYWORDS:
rcv : string limit to recevier. "sb" or "xb". default: all
verb : if set then print progress as you search
dirI[2] : strarr in case files moved .. (see masfilelist )
RETURNS :
ncross : the number of crosses we found
cfI[ncross]: {} array of structs holding Cross File Info.
DESCRIPTION:
This routine is written for the 12meter data.
What is a cross pattern:
A cross pattern will consist of:
- an optional calOn, caloff scan at the beginning.
- an azimuth scan
- an elevation scan.
So there can be 2 or 4 scans in a pattern.
There can als be multiple frequency bands (bms) in a scan (1 to 7).
Each scan/bm is stored in separate file so a pattern can have up to 4*7=28 files
associated with it.
All scans in the pattern will have the same pattern_id (equal to the
scannumber of the first scan in the pattern). This is used to group scans into
the pattern.
This routine will scan all files for a day and extract the files the
cross files belonging to the specified projectd id. It will the group all
of the files for a particular pattern into a cfI (crossFileInfo) structure. An
array of these structures will be returned. This aray can then be passed to
the cross processing routines.
Some tests are done on the scans found for each pattern:
- it must have an az and el scans.
- the number of beams for the az,el scan must be the same.
- it can optionally have calon,caloff scans.
- the number of beams in the calon,caloff scans must be the same
- If the scans in a pattern don't meet these requirements then it
is skipped. An example might be the azimuth strip completing, but the
source sets before the elevation strip.
--> Warning..this only finds files for 1 day. if the pattern straddles
midnight then it will be skipped.
An example of the returned struct array is shown below. These crosses
were made at xband and had 7 freqbands (bms) for each cross.
IDL> help,cFI,/st
PATID LONG 125300002 .. unique pattern id for pattern
SRCNM STRING '3C461' .. source name
NBMS INT 7 .. number of beams in this pattern
BMLIST INT Array[7] .. bm numbers 0..6
CFRAR DOUBLE Array[7] .. center freq (MHz) each band
FLISTCROSS STRING Array[2, 7] .. filenames: [az,el strips] for the 7 beams
SUMICROSS STRUCT -> Array[1] .. sumI struct for azstrip 1st bm
HASCAL INT 0 .. 1 if pattern has cal..if not ignore below
FLISTCAL STRING Array[2, 7] .. filenames [calOn,calOff] for the 7 beams
SUMICAL STRUCT -> Array[1] .. sumI struct for calOn of 1st beam
The sumI struct hold more detailed info that may be handy. Some of the info
is duplicated above. See masfilesum for a more complete description.
help,cFI.sumIcross,/st
OK INT 1
FNAME STRING '/share/pdata1/pdev/x101.20210910.b0s1g0.0010'...
FSIZE LONG64 38856960
NROWS LONG 120 --> number of rows in fits file
DUMPROW LONG 10 --> number of spectra/row
DATE LONG 20210910
BM INT 0
H STRUCT -> MASFHDR Array[1] --> fitsheader for the file
HMAIN STRUCT -> PDEV_HDRPDEV Array[1]
HSP1 STRUCT -> PDEV_HDRSP1 Array[1]
Example using it:
proj='x101'
yymmdd=20220406l
ncross=masfindcross(yymmdd,proj,cFI)
; process the crosses
mascrossproc(cFI,cI) ; to be written:)
(See /pkg/rsi/local/libao/phil/mas/masfindcross.pro)
NAME:
masfindmap - find the files that belong to an radecmap
SYNTAX: nstrips=masfindmap(proj,yyyymmdd,patternId,mapfI,$
noappbm=noappbm,dirI=dirI,bm=bm,band=band,grp=grp)
ARGS:
proj: string project id to search for. files will start with this. eg 'x101'
yyyymmdd[]:lonarr date(ast) to search for map. It is normally single element.
if your map spanned local midnite then you can enter an array of 2 dates.
pattern_id: long pattern id for the map.
KEYWORDS:
noappbm: if the data for separate beams has all been moved to a single directory then this
keyword should be set. Normally you will ignore it.
dirI[2] : strarr if the data has been move to non-standard directories then you can specify the
location here. see masfilelist for more details.
bm : int limit the returned data to a single bm/mock. numbers are 0 ..6
band : int limit the returned data to a single mock band. Single pixel observations
only take data in band 1. so setting band=1 will make things run a little faster.
grp : int limit data returned to grp 0 or 1. Multiple groups will be taken when doing
piggyback observing. Setting grp=0 will cause it to run a little quicker.Currently
the 12meter only uses grp=0.
RETURNS : nstrips
>0 number of strips in map that were found
=0 no data was found
-1 error trying to determine the cal,data layout..
DESCRIPTION:
Search for the files that belong to radecmap taken with the 12meter. The user specifies
the date yyyymmdd, project name, and pattern id. The routine will search all mock files on disc
that match these specs. The information is returned in a mapfI structure. The mapfI struct
can then be passed to masgetmap to load the actual map data into memory.
note: this routine returns info on the files,headers that belong to the map.It does not load
the data (see masgetmap() for that).
You can limit the file info that is returned by set the bm= keyword to a single bm.
mock files look like: x101.20220217.b0s1g0.nnnnn. b0 is the first bm or freq band recorded. bms go
0 to 6 depending on which mocks you used to record the data.
The mapfI struct contains:
IDL> help,mapfi,/st
** Structure <1c4b658>, 7 tags, length=1099240, data length=1085957, refs=1:
PATTERN_ID LONG 204800064 ; pattern id input by caller
SRC STRING '3C461' ; source name found in the file headers
NSTRIPS LONG 25 ; number of strips found in the map
NBMS LONG 7 ; number of bms recorded
CALPERSTRIP INT 2 ; number of cals per strip.can be 0,1,or 2
STRIPCODE INT 6 ; how was the map recorded;
1 - no cals fired
2 - winking cal embeded in the data
3 - cals only fired at start and end of the map
4 - cals fired only at start of each strip.
5 - cals fired only at end of each strip.
6 - cals fired at start and end of each strip
STRIPI STRUCT -> Array[25]; info for each strip found
The stripI structure contains info for each strip of the map. the entries will depend on
the number of cals fired during each strip:
the example below has cals fired at the start end end of each strip
IDL> help,mapfI.stripI,/st
** Structure <2586f38>, 4 tags, length=43968, data length=43437, refs=2:
STRIPNUM INT 0 strip number 0 .. nstrips-1
SUMICAL1 STRUCT -> Array[2, 7] cals before the strip
SUMIDAT STRUCT -> Array[7] the files hold the strip data
SUMICAL2 STRUCT -> Array[2, 7] the cals after the strip
for the cals, the first index is the calOn,caloff
the 2nd index is the beams
Note: when the cal is fired, 2 files are generated for each beam. a calON and cal off.
If a different number of cals is used, the the structures will change a bit.
stripcode=3 (cals at start and end of map only);
mapfI has 2 extra fiels containing the file info for the cals at start,end of the map:
calstart STRUCT -> Array[2, 7] . sumI for cals before map
calEnd STRUCT -> Array[2, 7] . sumI cals after map
for the stripI struct.
- if only one cal/strip then you have:
STRIPNUM INT 0 strip number 0 .. nstrips-1
SUMICAL STRUCT -> Array[2, 7] cals for stripo
SUMIDAT STRUCT -> Array[7] the files hold the strip data
- if no cals are fired then you'd just have:
STRIPNUM INT 0 strip number 0 .. nstrips-1
SUMIDAT STRUCT -> Array[7] the files hold the strip data
-----------------------------------------------
One question is where to find the pattern id.
the pattern id format is: ydddnnnnn.
y is the last digit of the year
ddd is the daynumber in the year
nnnnn is a number that increments by 100 each scan.
looking in the logfile it is the first scan after the line printing the MAPRADEC START:
2022-02-17_14:54:54 MAPRADEC
2022-02-17_14:54:54 src: ra,dec:232327.90 584842.00 rate:119.67 strip1:1 numStrips:25 code:0x30
2022-02-17_14:54:57 Mock scan starting filenumber:06100
2022-02-17_14:55:00 Observation begins Thu Feb 17 14:55:00 2022. scan:204800064
the pattern_id is 204800064
You could also use the masdmp program for particular day, proj,
masdmp --date=20220217 --proj=x101 --bm=0 --patid
It will dump all of the scans for the day. setting --bm=0 limits it to just 1 bm
(all bms have the same papattern id). scrolling through the output you'd find your map:
Proj Fnum BSG SOURCE PATID RA DEC C nSpc POL Chn PatNm TopFrq Bw RCV
x101 6000 b0s1g0 3C461 204800061 231557.6 574841 J 3 2 2048 RADECMAP 8219.0 172.0 xb
x101 6100 b0s1g0 3C461 204800064 231557.3 574842 J 6 2 2048 CAL 8219.0 172.0 xb
x101 6200 b0s1g0 3C461 204800064 231557.3 574841 J 6 2 2048 CAL 8219.0 172.0 xb
we started on 204800064 (the previous scans was testing)
-----------------------------------------------
EXAMPLE:
Suppose you made a map of 3C461 (casA) firing the cals at the start and end of each strip. The
system had been configured to over 1ghz bandwidth using 7 mocks.
proj='x101'
yyyymmdd=20220217
patid=204800064
nstrips=maspatfind(yyyymmdd,proj,patid,mapfI,band=1,grp=0)
you could then call masgetmap(mapfI,mapI) to load the actual data (once i write masgetmap :).
The sumI struct hold file info (see masfilesum).
help,sumI,/st
OK INT 1
FNAME STRING '/share/pdata1/pdev/x108.20080530.b0s0g0.00000.fits'
FSIZE LONG64 84389676 .. bytes in file
NROWS LONG 4 .. rows in fits file
DUMPROW LONG 319 .. dumps in first row
DATE LONG 20080530 .. yyyymmdd
BM INT 0 .. 0..6 for alfa
BAND INT 0 .. 0,1 low,hi band (at if)
GRP INT 0 .. grp 0,1
NUM LONG 0 .. file number .nnnnnn.fits
H STRUCT 1st row of fits file (without data,stat)
HMAIN STRUCT pdev main header
HSP1 STRUCT pdev sp1 header
(See /pkg/rsi/local/libao/phil/mas/masfindmap.pro)
NAME:
masfreq - return freq (or vel) array for a spectra
SYNTAX: retData=masfreq(hdr,retvel=retvel,restFreq=restFreq,velCrdSys=velCrdSys)
ARGS:
hdr : {} fits header b.h returned from masget();
retvel : if set then return the velocity array in the velocity coordinate system
restFreq:float if retvel is set, use this for the rest freq rather than the
value in the header. Units are MHz.
velCrdSys: string:user can specify a different velocity coordinate system from
that in the header.
Values are: 'T': topocentric,'G':geocentric,'B':barycenter,
'L;: lsr
RETURNS:
retDat[]: frequency array in Mhz for the points in the spectra
or the velocity array in km/sec
DESCRIPTION:
For returning velocity:
1. compute topocentric frequencies:freqTop
2. compute antenna velocity in velocity coord sys projected in ra,dec direction.
3. compute freq in velCrdSys :freqVcs=freqTopo*(1-antVelProj/c)
4. vel=c*(restFrq/freqVcs -1)
- user can specify a different restFreq from what is in the header
using the keyword restFreq
- currenly returns vel optical
- also assumes velocity << c.
- Uses the topocentric freq and the projected antenna velocity
for the reqested ra,dec. This comes from the jpl ephm every sec from pointing.
It does not use the velocity in the header
history:
25oct11: switched to used dindgen instead of findgen.
was always double anyway (from the constants).
(See /pkg/rsi/local/libao/phil/mas/masfreq.pro)
NAME:
masfreqdesc - return freq array for a spectra given the descriptor
SYNTAX: freqAr=masfreqdesc(desc,skycfr=skycfr,lo2offset=lo2offset,
lenfft=lenfft)
ARGS:
desc: {} returned by pdevopen
KEYWORDS:
skycfr: double sky cfr of band
lo2Offset: double offset lo2 from center of band. default: bw/2.
lenfft : long if timedomain data, then specify the len of the fft
DESCRIPTION:
This is an old routine prior to the fits header having the correct
frequency info. Now you should use masfreq(hdr).
RETURNS:
freqAr[]: frequency array in Mhz for the points in the spectra
(See /pkg/rsi/local/libao/phil/mas/masfreqdesc.pro)
NAME:
masgainget - return the gain given an mas header
SYNTAX: stat=masgainget(hdr,gainval,date=date,az=az,za=za,onlyza=onlyza)
ARGS:
hdr[n]: {hdr} header for data
KEYWORDS:
date[2]: intarray [year,dayNum] if provided, then compute the gain value
at this epoch (ast).
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 dataset.
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 mas record and then get the gain value
print,masgget(desc,b)
istat=masgainget(b.h,gain)
.. gain now has the gain value in K/Jy
input an entire file and compute the gain for all
records of dataset
print,masgetfile(filenam,bar)
istat=masgainget(bar.h,gain)
gain is now a array = to number of rows in file
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 date 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().
SEE ALSO:
gen/gainget gen/gaininpdata.pro
(See /pkg/rsi/local/libao/phil/mas/masgainget.pro)
NAME:
masget - input next mas rowfrom disc
SYNTAX: istat=masget(des,b,row=row,hdronly=hdronly,blankcor=blankcor,
float=float,double=double)
ARGS:
desc:{descmas} from masopen();
KEYWORDS:
row : if set then position to row before reading (count from 1)
if row=0 then ignore row keyword
hdronly : if set then just return the row header. no status or data.
blankcor : if set then correct for any a/d blanking. Blanked
spectra (which have fewer fft accums) will be scaled
to the unblanked value. This implies /float.
float : if set then force returned data to be floats
float : if set then force returned data to be doubles
RETURNS:
b: structure holding the hdr and data input
istat: 1 ok
: 0 hiteof
:-1 i/o error..bad hdr, etc..
DESCRIPTION:
Read the next row from a mas fits datafile pointed to by desc.
If keyword row is present, position to row before reading.
(See /pkg/rsi/local/libao/phil/mas/masget.pro)
NAME:
masgetcol - get a col from the fits file
SYNTAX: n=masgetcol(desc,colNm,data,rows=rows)
ARGS:
desc: {} struct returned from masopen()
colNm: string name of column. Same as help,b.h,/st structure tags
(except for DATE-OBS, MJD-OBS) see below..
KEYWORDS:
rows[2]: long rows to return. If not supplied then return all rows of
file. If single element then just return that row.
If two elements then first,last row to return.
Row numbering starts at 1.
RETURNS:
istat: n number of rows returned, -1 if error.
data[n]: date returned from file for this column.
DESCRIPTION:
Return individual columns from the fits files. The cols correspond
to the elements in the header struture returned by masget(desc.b) : b.h
The colNm passed in should match one of the strucuture tag names of
b.h. You can list them via: help,b.h,/st. The one exception is the
DATEXXOBS AND MJDXXOBS tagnames (see below).
By default all of the rows of the file are returned. You can limit
this by using the rows= keyword to limit the number of rows returned.
Row numbers starts at 1.
This routine is much faster than cycling through the file with
masget() since the spectral/stat data is not read in.
NOTES:
1. It can not be used to return the data stored in the heap
. stat, or data. see masgetstat or masget for that.
2. The two fits file column names: DATE-OBS and MJD-OBS are mapped to
idl structure names .DATEXXOBS, .MJDXXOBS since - is an illegal
structure tag name. You need to input the actual fits column names
for these two: DATE-OBS and MJD-OBS
3. use upper case.
4. If rows are specified they start at 1.
5. see desc.totrows for the number of rows in the file if you need it.
6. To find the column names, use masget(desc,b) then help,b.h,/st
(See /pkg/rsi/local/libao/phil/mas/masgetcol.pro)
NAME:
masgetfile - input an entire file
SYNTAX: istat=masgetfile(desc,b,avg=avg,median=median,ravg=ravg,$
float=float,double=double,filename=filename,fnmI=fnmI,$
blankcor=blankcor,remdc=remdc,$
; returned...
tp=tp,descRet=descRet,azavg=azavg,zaavg=zaavg,hdr=hdr)
ARGS:
desc: {} returned by masopen.. file to read (unless filename present)
KEYWORDS:
avg: if keyword set then return the averaged data
It will average over nrows and ndumps per row.
the returned data structure will store the averaged
spectra as a float. use /double if you want a double
see also : median
median: if set then return the median value for the file. This
reads the entire file into memory so you need room to
hold it.
ravg: if set then average each row
float: the averaged data should be returned as float.
if any averaging is done then the default is float.
double: The data should be returned as doubles.
filename: string if present then ignore desc. openfile,read, then close
the descriptor.
fnmI: {} filname Info struct returned from masfilelist. If present
then open, read, then close this file.
blankcor: if set then do the correction for blanking. The
blanked spectra will be corrected to have the same number of
spectral accumulations as the unblanked spectra.
If /avg,/ravg,or /median then blankcor is also enabled.
remdc : if set then interpolate across the dc channel for each spectra
RETURNS:
istat: 1 got all the requested records
: 0 returned no rows
: -1 returned some but not all of the rows
: -2 could not open the fie
b[n]: {} array of structs holding the data
tp[n*m,2]: float array holding the total power for each spectra/sbc
m depends on the number of spectra per row
descRet{}: if supplied and filename or fnmi is present, then the file descriptor will
be left open and returned in descRet
azavg : float return the average azimuth
zaavg : float return the average za.
hdr[] : string if supplied then return header from fileopen
(See /pkg/rsi/local/libao/phil/mas/masgetfile.pro)
NAME:
masgetfiles - input a set of files
SYNTAX: istat=masgetfiles(fnmIAr,bar,avg=avg,median=median,
float=float,double=double,fnmAr=fnmAr,blankcor=blankcor,$
tpAr=tpAr,remdc=remdc)
ARGS:
fnmIAr[nfiles]: {} from masfileslist. files to input
KEYWORDS:
avg: if keyword set then return the average of each file rather than the individual recs.
It will average over nrows and ndumps per row.
the returned data structure will store the averaged
spectra as a float. use /double if you want a double
see also : median
median: if set then return the median value for the file. This
reads the entire file into memory so you need room to
hold it.
float: the averaged data should be returned as float.
if any averaging is done then the default is float.
double: The data should be returned as doubles.
fnmAr[nfiles]: string if present then ignore fnmIAr and use fnmAr to find the files
it should hold the complete filenames.:w
blankcor: if set then do the correction for blanking. The
blanked spectra will be corrected to have the same number of
spectral accumulations as the unblanked spectra.
If /avg,/ravg,or /median then blankcor is also enabled.
remdc : If set then interpolate across the dc channel in the spectra
RETURNS:
istat: 1 got all the requested files
: 0 found no files
: -1 did not correctly read one file
bAr[n,nfiles]: {} array of structs holding the data
tp[n*m,2,nfiles]: float array holding the total power for each spectra/sbc
m depends on the number of spectra per row
descRet{}: if supplied and filename or fnmi is present, then the file descriptor will
be left open and returned in descRet
DESCRPTION:
Input a set of files. This is normally used to input files from different
frequency bands that were taken simultaneously. The files must have the same
format and the same number of rows.
You would normally used fnmIAr from masfilelist() to specify the files.
If just have a list of filenames (with directories) then use the fnmAr= keyword
(See /pkg/rsi/local/libao/phil/mas/masgetfiles.pro)
NAME:
masgethwcal - input 25hz hardware cal data
SYNTAX: istat=masgethwcal(fnmI,calI,bonar,boffar,fnmAr=fnmAr,avgrow=avgrow,verb=verb,$
incrfreq=incrfreq,$
bavgAr=bavgAr,brmsAr=brmsAr,avgmed=avgmed,$
sumI=sumI,dcintp=dcintp,double=double,$
row=row,nrows=nrows,hdrar=hdrar)
ARGS:
fnmI[n] :{} files to input. returned from masfilelist
KEYWORDS:
fnmar[n] : int array of filenames to use. ignore fnmI
avgrow : int if set then return bonar,boffar averaged to
a row's worth of data (usually 1 second). If not
set then data is averaged to 1 cal cycle (on 20ms,off20ms)
freqorder: int if set then return bonar[*,n] in freq order
default is order in fnmI.
verb : if present then print 1 line file summary
incrfreq: if set then make all data returned in the
arrays increasing freq. For those bands with decreasing
freq, the data is flipped and the header is adjusted
accordingly.
dcintp : if set then interpolate across the dc channel
double : if set then return data as doubles. def:floats
avgmed : if set then use median when computing row and file averages
row : long start at row=row (count from 1) (def: start of file=1)
nrows : long number of rows to process. def: all rows starting from row
RETURNS:
istat : int >=0 nn= number of files with hdw winking cal returned
< 0 error:
-1 = can't open file
-2 = file has different number of rows
calI[nn]: {} holds info about each file input
bonAr[m,nn]: {} mas struct of cal on,m number of cal cycles
boffar[m,nn]:{} mas struct of cal off
bAvgAr[nn,2]:{} if present then also return averaged data for each file
[*,0]=calon, [*,1]=caloff
brmsAr[nn,2]:{} if present then compute rms by chan for each file
[*,0=calOn,1=calOff . If rowavg is provided, the rms
is computed after the row avg.
sumI[nn]: {} summary info each file
hdrAr[nn]:{} if provided then return ascii hdrar from masopen (had digrms)
DESCRIPTION:
input hardware winking cal (25Hz) data from files. Any files that are not
taken with the hardware winking cal (hdr.calsrc=2) are skipped.
User can specify the files using flist[n] or with the fnmI= keyword (see
masfilelist).
The program uses masgetwcal() to input the data. The processing is:
1. separate the calons, cal offs for each row.
average the ons and the offs separately for row.
When doing the average, discard the first and last spectra of
each cal transition. The 25hz cal has 20millisecs on, 20millisecs off.
if you sampled at 1 millisec then
avg 2-19 on, 2-19 off. This is to not include the transitions.
Any adc blanking will be corrected in the averaging.
2. Load the bonAr, boffAr first index with the row averaged values.
The average will normally be 1sec. or 25 cycles.
3. Each file must have the same number of rows...
4. fill in the calI[i] struct for info about this file.
For this to work, the spectra need have integrations that are
divisible into 20 milliseconds. Since we throw out the 1st and last
you need at least 4 samples (5 ms each) 2ms per spectra is probably better.
This routine assumes that the integration of each spectra is a
multiple of 1 millisecond (it rounds the exposure to the nearest millisec).
(See /pkg/rsi/local/libao/phil/mas/masgethwcal.pro)
NAME:
masgetm - read multiple rows from a mas fits file
SYNTAX: istat=masgetm(desc,nrows,b,row=row,avg=avg,ravg=ravg,tp=tp,$
blankcor=blankcor,float=float,double=double, $
azavg=azavg,zaavg=zaavg)
ARGS:
desc: {} returned by masopen
nrows: long rows to read
KEYWORDS:
row: long row to position to before reading (cnt from 1)
if row=0 then ignore row keyword
avg: if keyword set then return the averaged data
It will average over nrows and ndumps per row.
the returned data structure will store the averaged
spectra as a float. use /double to get a double returned
ravg: if set then average within a row
hdronly: if set then only headers (row) will be returned.no data,stat.
blankcor: if set then rescale any a/d blanking spectra to have the same
number of accumulations as the unblanked data. This
implies /float. /avg or /ravg will also turn on blankcor
float : if set then return spectra as floats instead of
default type (int,long).
RETURNS:
istat: 1 got all the requested records
: 0 returned no rows
: -1 returned some but not all of the rows
b[n]: {} array of structs holding the data
tp[n*m,2]: float array holding the total power for each spectra/sbc
m depends on the number of spectra per row
azavg : float If /avg then the average azimuth
zaavg : float If /avg then the average zenith angle
DESCRIPTION:
Read in multiple rows from the requested data file. Return the
data in b[m].
If /blankcor is set then rescale any a/d blanked spectra to have the
same number of accumulations as the unblanked data. This will return
the data as floats.
If /avg is set, then average all of the rows into a
single spectra. If /ravg is set then average the dumps in each row
returning a single spectra for each row. Both of these averages
will do the blankcor before averaging.
If tp= is specified, also return the total power for each spectra.
The /float keyword will force the returned data to be floating point
(instead of the native datatype read from the file).
Any type of averaging (or /blankcor) will automatically return
float data.
NOTES:
- The last row of the file could have fewer dumps than the rest of the file.
If last row is read and it has different dumps than the previous rows
in the request, then it will be ignored.
(See /pkg/rsi/local/libao/phil/mas/masgetm.pro)
NAME:
masgetmap - input data from a radecmap
SYNTAX: nstrips=masgetmap(mapFI,mapDI,edgefract=edgefract,median=median,$
calIAr=calIAr,cntsToKar=cntsToKar)
ARGS:
mapFI : {} structure returned by masfindmap(). contains file info for map
(since we may not have the correct cal values yet.
KEYWORDS:
edgefr[2]:fltarr fraction of spectra to ignore on each edge. default:[.06,.06]
if 1 number is supplied then it will be used for both edge.
median : if set then use the median rather than the mean for:
averaging cal records
averaging spectra to get total power.
remdc : if set then interpolate across dc for each spectrag
noscl : if set the don't tpar to kelvins
RETURNS:
istat : int num strips. < 0 --> error
mapDI : {} struct holding the map data that was input
calIAr[nsets,nstrips,nbms] {} return cal info.
if no cals then nothing is returned
if cals start/end map then nstrips=1
if cals in strip then all 3 dimensions returned.
cntsToKAr[nstrips,npols,nbms]: fltarr scale factors used to convert to kelvins.
if cal only at start, end map than nstrips=1
dividing tpAr by these values will return to spectrometer counts
(if you have questions about the cal scaling)
DESCRIPTION:
Input an radecmap taken with the 12meter. The user passes in the mapFI struct
(returned by masmapFI. The routine inputs the map and returns the info in the
mapDI structure.
the mapdI struct contains:
IDL> help,mapdI,/st
** Structure <12d5d08>, 3 tags, length=234465008, data length=234289829, refs=1:
NSTRIPS LONG 25
BAR STRUCT -> Array[77, 25, 7]; nsmp,nstrips,nbms
TPAR FLOAT Array[77, 25, 2, 7] ; nsmp,nstrips,npol,nbms
raPos dblarr array[77,25] ; raLocation deg
decPos float array[77,25] ; decPosition deg
------
Notes:
terminology
nsmp - samples per strip
nsets- number of cals per strip. 1 or 2
nstrips- number of ra strips in map
npol - number of polarizations (usually 2).. we don't add pols
nbms - number of separate freq bands recorded
The ordering of the above arrays is:
[nsmp or nsets,nstrips,npols,nbms]
if an array is missing one of these, just drop it from the order
--> need to fix for winking cal and multiple dumps/record
(See /pkg/rsi/local/libao/phil/mas/masgetmap.pro)
NAME:
masgetmapwc - input data from a winking cal radecmap
SYNTAX: nstrips=masgetmapwc(mapFI,mapDI,edgefract=edgefract,median=median,$
degGainVar=degGainVar,verb=verb)
ARGS:
mapFI : {} structure returned by masfindmap(). contains file info for map
KEYWORDS:
edgefr[2]:fltarr fraction of spectra to ignore on each edge. default:[.06,.06]
if 1 number is supplied then it will be used for both edge.
median : if set then use the median rather than the mean for:
averaging cal records
averaging spectra to get total power.
remdc : if set then interpolate across dc for each spectrag
noscl : if set the don't scale tpar to kelvins
degGainVar:int order of polynomial when fitting caldif for gain variation in time.
def:3
verb : if set then output progress of routine.
RETURNS:
istat : int num strips. < 0 --> error
mapDI : {} struct holding the map data that was input
DESCRIPTION:
Input an 12meter radecmap taken with the hardware winking cal. The user passes in
the mapFI struct (returned by masmapFI). The routine inputs the map and returns the
info in the mapDI structure.
the mapdI struct contains:
IDL> help,mapdI,/st
** Structure <12d5d08>, 3 tags, length=234465008, data length=234289829, refs=1:
NSTRIPS LONG 25
BAR STRUCT -> Array[77, 25, 7]; nsmp,nstrips,nbms
TPAR FLOAT Array[77, 25, 2, 7] ; nsmp,nstrips,npol,nbms
raPos dblarr array[77,25] ; raLocation deg
decPos float array[77,25] ; decPosition deg
------
Notes
------
1. Assume there is a gain variation coming after the cal injection that varies much slower than the
cal cycle (40ms).
2. assume that Tcal (from the diode) is stable with time (it's temperature stabilized).
Let
calOn: (Tsky + Trcv + Tcal)*Giflo(stable)*Gvar(t)
calOff: (Tsky + Trcv )*Giflo(stable)*Gvar(t) (= Tsys)
calDif=calOn-caloff=Tcal*Giflo(stable)*Gvar(t)
3.If we divide caloff (tsys) by CalDif at each 40ms cal cycle we can cancel Gvar(t) since we assume it
varies much more slowly than 40ms.
- But this will blow up the noise since the noise/mean value for calDif is a factor of 2 to 3
bigger than noise on calOff (since the noise is the same and Tcal but you put the ripples of the cal into the off, on banpasses .. but they are much less than
the iflobanpass..maybe later use some edge spectra caloff on the maps.
2. using spectra bon, boff.
bavgOn=avg(bon)
bavgOff=avg(bavgoff)
bavgCal=bavgon - bavgoff
bavgCalN=normalize bavgCal to unity (since we want spec cnts still)
flatten all of the spectra with bavgcalN
bonF=bon/bavgcalN
boffF=boff/bavgcalN
3. find channels with rfi;
- compute rms/mean by channel for a strip
- do a linear fit and exclude channels >3sigma
also exclude edgefract channels (4% on each each side)
4. compute TpOnAr, tpOffAr using good channels from bonF, boffF
5. compute tpCal(t) =tpONar(t) - tpOffar(t)
6. nth order poly fit to tpCal(t)
7. divide tponar,tpoffar by fit
8. Get calInK(f) for channels used (avg)
9. multiply TponAr,tpOffAr by calin kel to switch to degK
Notes:
terminology
nsmp - samples per strip
nsets- number of cals per strip. 1 or 2
nstrips- number of ra strips in map
npol - number of polarizations (usually 2).. we don't add pols
nbms - number of separate freq bands recorded
The ordering of the above arrays is:
[nsmp or nsets,nstrips,npols,nbms]
if an array is missing one of these, just drop it from the order
--> need to fix for winking cal and multiple dumps/record
(See /pkg/rsi/local/libao/phil/mas/masgetmapwc.pro)
NAME:
masgetscanhwc - input scan data using hardware winking cal
SYNTAX: nbms=masgetscanhwc(fnmI,scanI,edgefract=edgefract,median=median,$
degGainVar=degGainVar,verb=verb,row=row,nrows=nrows,filelist=filelist)
ARGS:
fnmI[n] : {} from masfilelist. files in scan (if multiplebeams)
currently assumes 1 file/scan/bm
KEYWORDS:
edgefr[2]:fltarr fraction of spectra to ignore on each edge. default:[.06,.06]
if 1 number is supplied then it will be used for both edge.
median : if set then use the median rather than the mean for:
averaging cal records
averaging spectra to get total power.
remdc : if set then interpolate across dc for each spectrag
noscl : if set the don't scale tpar to kelvins
degGainVar:int order of polynomial when fitting caldif for gain variation in time.
def:3
row : long . if supplied then limit scan to maxrow recs
nrows : long . if supplied then input nrows rows. def: entire file (starting
at row)
filelist[n]: strarr if supplied then ignore fnmI,use filelist
verb : if set then output progress of routine.
RETURNS:
istat : int num of beams (freq bands of data returned) < 0 --> error
scanI : {} struct holding the scan data that was input
DESCRIPTION:
Input an 12meter scan taken with the harware winking cal. The user passes in
an fnmI[nmbs] array with the file info for the files in theh scan. It currently works
for multiple bms/scan but only 1 file/bm.
This routines works similar to masgetmapwc() except that it inputs a single scan rather than
and entire map.
the scanI struct contains:
IDL> help,scani,/st
** Structure <14bbc08>, 7 tags, length=3143128, data length=3140282, refs=1:
NBMS LONG 7
NRECSFILE LONG 380
NCALCYCLESFILE LONG 9500
HDRAR STRUCT -> MASFHDR Array[380] hdrAr for bm0 of each rec
XH FLOAT Array[380] ; index 0..380-1
XHI FLOAT Array[9500] ; use to interpolate hdrAr data to calcycle sample times
BMI STRUCT -> Array[7] ; info from each beam.
IDL> help,scani.bmi,/st
** Structure <1522548>, 11 tags, length=391256, data length=391230, refs=2:
TPCALDIF FLOAT Array[9500, 2] totPwr calOn - calOff raw units
TPCALDIFFIT FLOAT Array[9500, 2] fit to tpcaldif
TPARONK FLOAT Array[9500, 2] totPwrcalon in deg K
TPAROFFK FLOAT Array[9500, 2] totPwrcaloff in degk
BONAVG STRUCT -> Array[1] ; avg calon spc degK (no bpc)
BOFFAVG STRUCT -> Array[1] ; avf caloff spc degK (no bpc)
BRMSAR STRUCT -> Array[2] rms/mean for calon,caloff spectra
MASKAR INT Array[512] 1-> channel used,0 not used in totpwr computation
EDGEFRACT FLOAT Array[2] 0..1. fraction of edge to exclude in totPwr computation
NCHANUSED LONG 465 # channels used for totpwr computation
BWUSED FLOAT 156.240 bandwidth used MHz.
CNTSTOK FLOAT Array[9500, 2] convert counts to degK for each calcycle
CALKUSED FLOAT Array[2] cal values used in degK
Examples:
* Suppose you use the hardware winking cal for an on, off observation
You could then call this routine for the Position on scan and the position off scan separately.
* scanI.xh and scanI.xhI can be used to interpolate header information sampled once a second
to the time resolution of the cal cycles
eg:
azPosI=interpol(scanI.hdarar.azimuth,scanI.xh,scanI.xhI)
Notes
------
how the gain correction works: see masgetmapwc() description.
Some steps in the processing:
Output spectra
scanI.bmI[nbms].bonavg .. is the average of each calon cycle spectra after converting to K
with the correction factor for each calcycle.
scanI.bmI[nbms].boffavg .. is the average of each caloff cycle spectra after converting to K
with the correction factor for each calcycle
(See /pkg/rsi/local/libao/phil/mas/masgetscanhwc.pro)
NAME:
masgetstat - input status info from file
SYNTAX: istat=masgetstat(desc,stat,nrows=nrows,all=all)
ARGS:
desc: {} struct returned from masopen()
KEYWORDS:
row : long row to position to before reading count from 1.
nrows: long number of rows of stat data to read in .default=1 row
from current position.
all: if set then read in the stat info for the entire file.
RETURNS:
istat: n number of stat stuctures returned. One struct for each
spectra. if nrows were read, there will be nrows*spcPerRow
-1 error reading file
stat[n]:{stat} stat data from file.
(See /pkg/rsi/local/libao/phil/mas/masgetstat.pro)
NAME:
masgetwcal - Inp Row, compute CalOn,calOff for winking cal
SYNTAX: istat=masgetwcal(desc,bon,boff,toavg=toavg,bin=bin,row=row,$
nrows=nrows,lastonly=lastonly,double=double)
ARGS:
desc : {} from masopen()
KEYWORDS:
row: long position to row before reading..`
count from 1. row=0 --> start current position
toavg: int number of calcycles to average. This must divide into the
number of calcycles in a row. If it is greater than the
number of calcycles in a row, then all of the calcyles in
a row will be averaged (eg.. toavg=999) will probably
sum all of them.
nrows: int number of rows to return. def=1
If more rows are requested than in file, then the only
the number of rows left is returned. No error is reported.
lastonly: int if set then only drop the last sample of each group of
calons, or cal offs. This works if the cal phase was set to
50% (where the transition occurred).
double : if set return spectra as doubles. default is float
RETURNS:
istat: n number of time samples in structure
-2 winking cal not enabled..
-3 record does not have at least 1 cal cycle
-4 spectra in row not a multiple of calCyleLen
or requested toavg does not divide into
cals cycles in row.
-5 calon,caloff out of sync.. not 10,10...
bin[n] : {} if supplied then also return the input data.
DESCRIPTION:
Normally winking cal data is taken at a high data rate with
multiple winking cal cycles within 1 row of data. This routine will
input a rows worth of winking cal data and compute Tsys + Calon,Tsys +calOff
When computing the cal, ignore the first and last spectra of each
calon or caloff.
If there are 10 specta with calOff followed by 10 spectra with calOff
then sum spectra 2-9 of each before computing the calon,caloff
NOTES:
1. The output spectra are summed (not averaged) over ncalon-2 and ncaloff-2
spectra.
2. If blanking has been enabled, then the spectra with fewer than normal
counts will be scaled up normalSum/blankedSum. The values
3. For now, the routine wants an integral number of calCycles in
a row...
4. The .st status portion of the structure will have the adcOverflow,pfbOverflow
and the fftaccum adjusted for the number of integrations (summed).
5. the header entry exposure will be increased to reflect the
accumulation time in the returned spectra. eg
suppose:
- 10 millisecond dumps
- 50 millisecond cal on, cal off
- average 10 cal cycles
- the exposure time goes from .01 secs to:
10ms*3calon*10avg=300 millisecs
5. The accumulation of the calOn, caloff cycles reduces the number of
spectra in the output data structures. eg. suppose you have
2048 channel, stokes data with 500 spectral samples per row, and
20 spectral time samples in the cal cycle (10 caloff , 10 calon). Then
the dimensions of the input, output data are:
bin.d[2048,4,500] .. input
bsum.d[2048,4,25] .. 500/20
bdif.d[2048,4,25] .. 500/(20)
6. If blanking is enabled and bin=bin returns the input structure
note!! the blanked spectra have already been scaled to the average values
(it is not the raw input buffer).
(See /pkg/rsi/local/libao/phil/mas/masgetwcal.pro)
NAME:
masimgdisp - display a set of correlator records as an image
SYNTAX: img=masimgdisp(b,nsigclip=nsigClip,clip=clip,pollist=pollist,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,desc=desc,han=han,maxrecs=maxrecs,mytitle=mytitle,$
hlind=hlind,hlval=hlval,hldash=hldash,hlvlines=hlvlines,$
useind=useind,ln=ln,chn=chn,freqincrease=freqincrease)
ARGS:
b[nrecs]: {corget} correlator data to make image of
RETURNS:
img[nchns,nrecs]: float image of last sbc displayed (before clipping).
KEYWORDS:
nsigClip:float if suppled then clip data to [-sig,sig]*nsigClip
This is done after any averaging for -zx,-zy.
If this is supplied, then clip= is ignored.
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].
pollist: long A single number whose decimal digits specify the pols
to display. eg 1,12,1234 .. pols,3,4 are stokes u,v
The default is pola,polb.
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 masimgdisp 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:{masget} 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
desc : {} if provided,then routine will input data into
b[] before making image. desc comes from masopen().
It starts reading from the current row. The number
of spectra in image defaults to 600(for polA and 600
for polB). You can change this with the numspc keyword
numspc: int number of spectra to display from b (or to read from the
file. The default if b is provided is the entire array.
If reading from the file then we will use about 600 spc.
han: if set and scan keyword set, then hanning smooth the
data on input.
mytitle:string user supplied tittle instead of scan,srcname,az,za
az,za at top of the plot.
yr[2]: float specify min,max label values for y axis. Default is
spectral count
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
freqincrease: If set then make the output image increasing frequency. This is
done just before output so all other keywords, values should use
the native order of the data eg. bpc should match the input
data. The flipping is only done on the displayed image.
The returned image is in the default freqorder
extra_=e this allows you to input keywords that will be
passed to the plotting routine. eg title=title..
DESCRIPTION:
masimgdisp will display a set of mas spectra as an image.
By default it will make a separate image for each polA and polB.
The polList keyword lets you specify polA or polB (and someday maybe
stokes U,V.
You can input the data outside of this routine (eg with masgetm) or
use the desc,numspc keywords to have masimgdisp input the data directly.
By default, the image is bandpass normalized by the average of all the
records (sbc/avg(pol) - 1). If the median keyword is used then
avg(sbc) is replaced by median(pol). The bpc keyword allows you to input
a separate mas record to use as the normalization.
The col keyword lets you also flatten the image in the time
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 time 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
pol displays. The multiple pol displays are always scaled to fit
inside the window (700 by 870 pixels).
The desc and numspc keywords can be used to input the data directly from
a file. In this case the desc comes from the masopen() routine and numspc
is the number of spectra to put in the image.
After displaying the image, use xloadct to manipulate the color table.
The routine returns the last image displayed (before clipping).
EXAMPLES:
input data and make an image
istat=masgetm(desc,30,b,/han0 .. input 30 rows with hanning smoothing
.. assume there are 10 spec/row
1. display the image of pola and polB
img=masimgdisp(b)
2. display only polB, scale y by 2, and x by -2
img=masimgdisp(b,pol=0,zx=-2,zy=2)
3. display pola,polB, clip to .015 Tsys , median filter the
bandpass correction:
img=masimgdisp(b,/median,clip=[.015])
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=masimgdisp(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/mas/masimgdisp.pro)
NAME:
masinterpfreqbands - interpolate multiple freq bands to 1 spectra
SYNTAX: masinterpfreqbands,bar,df,dI,f0=f0,f1=f1
ARGS:
bar[n]: {} array of mas structs holding multiple frequencies
df : float delta freq (Mhz) for final interpolated spectrum
ARGS:
f0 : MHz if supplied then use this as the lowest interpolated chan. Def is
the minimum in the initial data.
f1 : MHz if supplied then use this as the maximum interpolated chan. Def is
the maximumin the initial data.
RETURNS:
Di : {} struct holding freq and interpolated data (see below for definintion).
TERMINOLOGY:
bms: when using alfa, the 7 mock boxes were called bms.This terminology
has carried over into single pixel observing even though they all point
at the same location.
bands: With alfa observing, bands refered to the upper and lower 150 Mhz freq bands
taken with one mock box. For single pixel observing the lower 150 Mhz band
is not used and the term band is often used in place of bms to refer to
masinterpfreqbands,reform(bar[i,*]),df,dI
different frequency bands.. (sorry for the confusion:)
DESCRIPTION:
Single pixel observing normally has up to 7 different frequency bands (or bms) that
are typically overlapped. This routine will interpolate these bands to a single spectra.
The user supplies the mas (bar) struct to combine and well as the requested frequency
spacing for the final spectra.
Most times the bands are overlapped with the overlapping region having the bandpass filter
falloff.
Prior to 28feb23 the excloverlap keyword specified how much of this overlap region should be
excluded from the edge of each band. this caused problems when the bandpasses were
still present. The interpolation could then jump back and forth between adj values from
the different bandpasses.
The current version just splits the overlap region in 2, taking the bottom part from the
lower freq spectra, and the upper half from the higher freq spectra.
The routine will now work with a record with multiple dumps. If you have multiple records to interp,
you must call this routine multiple times
If there is no overlap between a set of bands then that portion of the spectrum
will be a linear interpolation between the endpoints.
The returned structure dI will have:
dI.npnts ; points in a spectrum
dI.npol ; 2nd dimension of data
dI.freq[npnts]: dbl frequency array
dI.d[npnts,npol] ; data
or dI.d[npnts,npol,ndump] if multiple dumps per record
-->Note:
This routine just does interpolation between points (idl interpol) it does not do averaging.
If you specify a much larger freqspacing than the input spectra then the returned signal will not
be averaged over the wider freq step.
(See /pkg/rsi/local/libao/phil/mas/masinterpfreqbands.pro)
NAME:
masisp12m - return 1 if p12m data
SYNTAX: istat=masisp12m(hdr)
ARGS:
hdr: {} mas.h header
RETURNS:
istat: 0 305m data, 1 =p12m data
use date_obs from header. > 2021_01 is p12m
(See /pkg/rsi/local/libao/phil/mas/masisp12m.pro)
NAME:
maslist - one line summary listing of file(s)
SYNTAX maslist,desc,fnmI=fnmI,flist=flist,sumI=sumI,patid=patid
ARGS:
desc : {} descriptor from masopen() for file to list. Ignored if fnmI or flist
keywords present
KEYWORDS:
fnmI[n]:{} structure returned by masfilelist. If this keyword is
present then use this for the files to list
flist[]:strarr array of filenames to list.
patid : if set then print pattern id rather than scan id
RETURNS:
sumI[n]: {} summary info for mas files.
DESCRIPTION:
Make a 1 line summary listing for file. You can specify the files with:
1. fnmI=fnmI[] array fileInfo structures returned from masfilelist()
2. flist=flist[] strarr of filenames to list
3. desc a descriprtor from the call to masopen().
If more than one of the above is specified, then the order is
1 then 2 then 3.
The 1 line listing contains:
SOURCE SCAN RA DEC C nSpc POL Chn PatNm TopFrq Bw RCV
Cur-pos 931000143 180506.0 85109 J 3000 2 8192 RUN 1450.0 172.0 alfa
SEE ALSO
masfilesum()
(See /pkg/rsi/local/libao/phil/mas/maslist.pro)
NAME:
masmakemask - make a mask for good channels using a dataset.
SYNTAX: istat=masmakemask(fnmI,maskI,verb=verb,nsig=nsig,edge=edge,flist=flist)
ARGS:
fnmI[n]: {} array holding files info to process (see masfilelist())
KEYWORDS:
nsig : float number of sigmas that separate good from bad channels
edge[2]: fltarr fraction of each edge of bandpass to ignore def:[.06,.06]
verb : if set then plot rms/mean by channel for each file.
flist[n]: strarr filenames to use (including path).
if supplied then ignore fnmI[] and use flist
RETURNS:
istat : int 1 - ok, -1 error
maskI : {} struct holding the mask info (see below).
DESCRIPTION:
Given a set of datafiles, this routine will compute a mask of good freq channels to use
(when computing total power, or fitting a baseline). The processing steps are:
(assume there are N files)
for each file in fnmI[N]
- input the file
- compute rms/mean by freq channel using the spectra in the file.
channels with rfi will have a larger rms/mean .
- compute the avg spectra for the file
After all N files are input:
- compute the average rms/mean for the nfiles. call it:rmsAvg
- compute the peak value of rms/mean for each file (take the peak from each channel):call it: rmsPk
Using these 2 averages, find the channels with a large rms/mean:
-Do a robust linear fit to rmsAvg(nchan) and rmsPk(nchan).
Any fit residuals greater then nsig will be flagged as bad.
maskI returns the following info
IDL> help,maskI,/st
NCHAN LONG 1024 ; number of freq channels
NFILES LONG 31 ; number of files used
NSIGUSED FLOAT 2.00000 ; value of nsig used for clipping
MASKAVG STRUCT -> Array[1] ;holds maskInfo for the avgerage mask (see below)
MASKPK STRUCT -> Array[1] ;holds maskinfo for the peakHold mask (see below)
BRMSAVG STRUCT -> Array[1] ;average rms/mean spectra (display with masplot,brmsavg
BRMSPK STRUCT -> Array[1] ;peakHold rms/mean spectra
BAVG STRUCT -> Array[1] ;average spectra
The maskavg, maskpk struct holds:
help,maski.maskavg,/st
SIGA FLOAT 8.52951e-05 ;polA sigma from robust fit to rms/mean avg (or rms/mean peak)
SIGB FLOAT 8.04953e-05 ;polB sigma from robust fit to rms/mean avg (or pk)
NGOODA LONG 382 ;# of good chan in polA
NGOODB LONG 502 ;# of good chan in polB
NGOOD LONG 327 ;# of good chan in polA and B
MASK INT Array[1024] ;mask for channels 0=bad, 1=goood,
The routine makes a mask using the average of the rms/mean by files.
It makes a second mask using the max (in each chan) of the peakHold rms/mean spectra
The pkMask will most likely exclude more points than the average mask.
The /verb keyword will plots:
Each rms/mean will be plotted as it is processed.
- you should set the ver scale to whatever the limits will be (from the radiometer equation):
eg: 172MhzBw/1024chan rms/mean = .008 so ver,0,.02 will work
- after plotting all the files it will plot the average rms/means
top frame: rms/mean avg
- white + :polA all points
- red + :polB all points
- green + :polA good freq channels
- yellow + :polB good freq channels
If you see green or yellow points that stick up above the rms/mean avg then you may
want to recall the routine with a smaller value for nsig (def=2). But be careful, if you
make nsig too small, then you may run out of good points.
Some requirements:
- each file should have the same number of channels.
- you probably don't want to mix different freq bands.
EXAMPLE:
Suppose we took a number of sband crosses on 3C295 on 11nov22
- there is a calOn,off at the start of each cross.
- use the calOff's to compute the rfi mask.
proj='x101'
yymmdd=20221111l
nfiles=masfilelist('',fnmI,proj=proj,yymmdd=yymmdd,bm=0,band=1,grp=0,/appbm)
;.. get summary info (it includes the headers) to select the cal Offs
istat=masfilesum('',sumI,fnmI=fnmI)
ii=where((sumI.h.obsmode eq 'CAL') and (sumI.h.scantype eq 'OFF'),nfiles)
; make a subset of only the cal offs (note sband only has 1 mock band)
fnmICaloff=fnmI[ii]
; make the mask, setup hor ,ver for /verb plotting.
; exclude 8%. 10% from the 2 edges (bw of rf < mock bw)
hor
ver,0,.02
edge=[.08,.1]
istat=masmakemask(fnmICalOff,maskI,/verb,edge=edge0i
; use the peak hold spectra
mask=maskI.maskavg.mask
; call stripchart routine to monitor with mask
masmonstripch,mask=mask,proj=x101
(See /pkg/rsi/local/libao/phil/mas/masmakemask.pro)
NAME:
masmath - perform math on correlator data
SYNTAX: bout=masmath(b1,b2,sub=sub,add=add,div=div,mul=mul,$
avg=avg,median=median,polAvg=polAvg,$
sadd=sadd,ssub=ssub,sdiv=sdiv,smul=smul,$
norm=norm,mask=mask,double=double)
ARGS:
b1[n] : {masget} first math argument
b2[n or 1]: {masget} second math argument
KEYWORDS:
2 Vector
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
1 vector
avg: if set then bout= avg(b1)
median: if set then bout= median(b1)
polavg: if set then avg the pols this can be called
avg, or med.
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 if:
:b1,b2, and add,sub,mult,div
normalize (mean) each elm of b2 before op.
:b1, /norm and no add,sub,mul,div
Normalize (mean) each elm of b1
:b1,b2 and /norm and no add,sub,mul,div
Error
:scaler op then /norm does nothing
mask :{masmask} (not implemented yet)
a mask to use with norm keyword (see masmask).
double: if set then make data array double rather than float.
RETURNS:
bout[n]: {masget} result of arithmetic operation.
with d not float rather than long
DESCRIPTION:
Perform arithmetic on mas data structures. The operation performed
depends on the keyword supplied.
VECTOR OPERATIONS:
-op=/add,/sub,/mul,/div
bout=b1 op b2 on a spectra by spectra operation.
If /norm is also set then each spectra of
b2 is normalized (mean) before the operation.
SCALER OPERATIONS:
-op= /ssub,sadd,smul,sdiv
: bout=b1* op
: Scaler can be a single value. This number will be used
for all operations,
: Scaler can be fltarr(npol) equal to the number of pols
Then each polarizion will used its on value:
eg: bout.d[*,ipol]=b1.d[*,ipol] op scalVal[ipol]
MISC OPERATIONS:
-op=/norm with no add,sub,mul,div specified
:just normalize each spectra of b1 (mean):
eg: bout[i].d[*,ipol]=b1[i].d[*,ipol]/mean(b1[i].d[*,ipol])
-op= /avg,/median
: compute the average or median by pol
eg: bavg.d[*,ipol]=mean(b1.d[*,ipol],dim=2) avg over spectra
EXAMPLES:
let bon[300] be position switch on data
boff[300] be position swith off data
1. avg on,off
bonavg=masmath(bon,/avg) or masmath(bon,/median)
boffavg=masmath(boff,/avg) or masmath(boff,/median)
2. on/off -1
bonoff=masmath(bonavg,boffavg,/div)
bonoff=masmath(bonoff,ssub=1.)
3. Normalize each 1 sec on by the 1 sec off
bonoffM=masmath(bon,boff,/div,/norm)
4 Normalize each bon rec to unity
bonN=masmath(bon,/norm)
5. use normalized average off to remove bandpass for each
rec of bon
bonFlat=masmath(bon,boffavg,/div,/norm)
6 remove a smoothed bandpass from each spectra of bon
;; compute bsmo somehow.. then
bonSmo=masmath(bon,bsmo,/sub)
(See /pkg/rsi/local/libao/phil/mas/masmath.pro)
NAME:
masmkfnm - create filename for fnmI struct
SYNTAX: fnmINew=masmkfnm(fnmI)
ARGS:
fnmI: struct fnmI struct from masfilelist
RETURNS:
fnmIN: {} structure holding the new filename
DESCRIPTION:
masmkfnm will take the parsed info from fnmI and create a new
filename and load it into a new fnmI struct. This allows you to
create a new fnmI struct by just changing one of the parsed parameters.
This routine will then fill in the .fname entry to reflect the new values.
The format is:
dir/proj.date.obs.brdId.seqnum.mas
where obs is optional
* Structure PDEVFNMPARS, 8 tags, length=60, data length=60:
DIR STRING '/share/pdata1/pdev/'
FNAME STRING 'if2noise.20081107.b0s0g0.00000.fits'
PROJ STRING 'if2noise'
DATE LONG 20081107
SRC STRING ''
BM INT 0
BAND INT 0
GRP INT 0
NUM LONG 0
(See /pkg/rsi/local/libao/phil/mas/masmkfnm.pro)
NAME:
masmkstruct - make a mas structure from template
SYNTAX: bnew=masmkstruct(btempl,float=float,double=double,ndump=ndump,$
npol=npol,nelm=nelm)
ARGS:
btempl: {masget} template struct to use
KEYWORDS:
keywords can override settings of templ
float: make data float
double: make data double
ndump : int make ndump spectra per entry
npol : int make npol pols
nelm : int make nelm entries in array
RETURNS:
bnew[nelm]:{masget} struct with accum entry
DESCRIPTION:
Create masget struct using btempl as a template. Override
things in btempl with the keywords.
(See /pkg/rsi/local/libao/phil/mas/masmkstruct.pro)
NAME:
masmon - monitor mas files.
SYNTAX: masmon,bsg=bsg,num=num,projid=projid,date=date,noavg=noavg
norm=norm,pollist=pollist,desc=desc,dirI=dirI,
noappbm=noappbm,flip=flip
and masplotkeyords,usepixwin=usepixwin
ARGS: none
KEYWORDS:
bsg: string beam,band,group to monitor. Format is bBsSgG where:
B =0..6 is the beam
S =0..1 is the band (0=1450 band, 1=1300 band)
G =0..1 is the group.
default:'b0s0g0'
num: int the file number to start with. (The default is the most recent)
projid: string project id (first part of filename) to monitor. The default
is the most recent one written.
date : long The date (in the filename) to start monitoring. The
default is the current date. If you enter "*" it will allow
any date. the format is yyyymmdd or "*"
noavg: If set then don't average multiple spectra within one row.
By default each row is averaged to one spectra.
desc : {} You can return from masmon with the current file still
open (using the x command from the menu). The file descriptor
will be returned here (this is actually a structure used by the
masxx routines). It is then your responsibility to close the
file via masclose,desc or masclose,/all.
dirI : strarr(2): directory info. If the files are in a non standard place
then pass in location in dirI
dirI[0]=prefix for dir: /share/pdata
Num 1..n goes here
dirI[1]=postfix for dir: /pdev
noappbm: if set then no bm number between pre,pos dir name
eg.. /share/pdata/pdev
flip : if set then flip the freq array of the plot.. use if
hdr info is wrong
PLOTTING
norm: if set then normalize spectra before plotting
pollist: int pols to plot: 1,2,12,123,1234. -1 --> all pols available
off: float if multiple spectra/plot (/noavg and multiple spectra per
row) then separate spectra by this amount vertically.
masplot keywords
--> any masplot keyword can also be used.. it will be passed on to masplot. eg
retvel
velcrdsys
restfreq
-----------------------------------------------------------------------------
DESCRIPTION:
Monitor online and offline mas files plotting each row in the file. The
files to monitor can be selected by:
1. if no arguments are supplied, use the most recent files from
beam0,band0,group0.
2. Specify the projid,date using the startup keywords.
3. Use the menu inside masmon to select a different set of files to monitor.
- Enter the menu once masmon is started by hitting any key.
-----------------------------------------------------------------------------
EXAMPLES:
masmon - start monitoring most recent file for b0s0g0
masmon,bsg=b0s1g0 .. start with beam0, band1 , group0
masmon,projid='a2130' .. use a2130 files
masmon,projid='a2130',/noavg .. do not average the 2 millisecond spectra
masmon,pollist=12 .. only only plot pola and B, ignore, U,V.
-----------------------------------------------------------------------------
USING THE INTERNAL MENU:
Enter the internal memory by hitting any key.
The menu is:
Cmd CurVal function
--> PLOTTING:
p 1234 pols to plot (1..4)
h low high set horizontal plot range ..no args --> autoscale
v low high set vertical plot range ..no args --> autoscale
delay secs secs to delay between plots..
--> POSITIONING:
bsg b0s0g0 beam,band,grp to use (change all at once)
bm 0 beam to plot
band 0 band to display (0 = 1450, 1=1300)
d 20081011 date to use. First number of data displayed
num 00801 number.. file number to display
r 1 position to row. maxRows: 99
l (all) list files this date (all --> all dates)
--> MISC:
hdr display last header read
s 0 step mode. 1-on,0-off
z stop in procedure (for debugging)
q to quit
x quit leaving the current file open. return description in desc keyword
otherKeys continue plotting
Enter the cmd and a possible value. The program will then continue.
The menu options are:
p 1234 .. pol to plot. By default all pols in the file are plotted. The
.. The pols are numbered: 1-polA,2-polB,3-U,4=V. Any combination
.. of the pols can be displayed (entered as a one number eg 12..)
Example:
p 12 .. just plot pols 1(a) and 2 (b)
h,v .. These will change the horizontal,vertical scale of the plots
Example:
h 1400 1450 display frequencies 1400 thru 1450 (if they are in the band)
delay .. secs to delay between plots. Use if the plotting is going
too fast.
Example:
delay .1 Delay .1 seconds between plots
bsg b0s0g0 .. Select the beam,band,group to display.
.. By default beam0,band0,grp0 is displayed.
Example:
bsg b0s1g0 .. will display the band centered at 1300 Mhz.
d 20080826 .. Move to first file on this date using the current
project id.
n 100 .. Move to this filenumber for the current date and project id.
If the number does not exist, you remain on the current file.
l (all) .. list files for this projid, data. If all is included then
list files for this projid for all dates.
r 1 .. position to a row in the current file (count from 1)
hdr .. print last header read
s 0/1 .. step mode will stop and display the menu after each spectra is
plotted. s 1 turns on step mode. s 0 turns it off.
z .. stop in masmon routine (for debugging)
q .. exit masmon.
(See /pkg/rsi/local/libao/phil/mas/masmon.pro)
NAME:
masmon2dfits - plot 2d gauss fits from mascrossfit2d
SYNTAX: masmon2dfits,fitiAr,delay=delaysecs,freqtouse=freqtouse,pwr=pwr,$
remtsys=remtsys
ARGS:
fitiar[n] : {} array of fit structs from mascrossfit2()
freqToUse: Mhz freq to use. def=all bands
pwr : if set then data in in power units. def: degK
(used for yaxis labels)
remTsys : if set then remove the constant Tsys term from the data
before plotting.
DESCRIPTION:
The mascrossfit2d() routine returns an array of structures with the fitting
info. This routine will plot the input data and resulting fit. Each plot will
have the data from 1 freq band and fit. There will be 3 frames on a plot:
- upper frame: az strip data vs offset with the 2d fit for az overplotted.
- middle frame: el strip data vs offset with the 2d fit for el overplotted.
- bottom frame: the az and el data and fits are plotted vs the el offset from
the center of the cross. The fitting routine includes a
linear term in elevation when fitting.
There are a number of parameters that will control what and how things
are plotted. After entering the routine, hitting any key will bring up a
menu that lets you change these parameters. You can then enter a command
and value.
CMD ARG Description
i index move to index (0 to nfits-1= 209)
src all srcToPlot,all or
3C274 3C286 3C348 3C405 CasA
f all Freq to plot or all. frqList:
8223 8363 8505 8647 8789 8931 9073
delay secs delay between plots. def:.1 secs
r 0,1 remove tsys:0,1
s 0,1 single step (0,1)
v min,max .. set vertical scale, blank is autoscale
d debug. stop in routine
q quit monfits
any otherChar continue plotting
currently stepping: 0
i : lets you skip around in the fitiar[] array
src : you can limit the plots to a single source. Entering "all" will
i plot all sources
f : limit plots to this freq. Entering "all" will return to all frequencies.
delay: Put a delay of delay secs between each plot (fractional secs ok).
r : 1-> remove tsys from data and fits before plotting. 0.. don't remove
Since it uses the fit value,if the fit failed, the plot may not appear.
s : 1-> single step the plots. It will redisplay the menu after each plot and wait for
you to enter a key.
0-> go back to continuous mode.
v : vmin,vmax set the vertical scale for the plots. Entering a blank
will return to autoscaling.
d : go to debug mode. will stop within the routine.
enter .continue to continue on.
q : quit this routine.
The routine uses masplot2dfit() to plot a single fit.
Example:
; if you've saved the fits somewhere
restore,'fitiar.sav' ,/verb
% RESTORE: Portable (XDR) SAVE/RESTORE file.
% RESTORE: Save file written by phil@megs3, Fri Oct 7 08:55:20 2022.
% RESTORE: IDL version 7.1 (linux, x86_64).
% RESTORE: Restored variable: FITIAR.
masmon2dfits,fitiar
.. plottings starts
.. you can resize the window with the cursor to get a better view
.. hit space bar to bring up the menu:
CMD ARG Description
i index move to index (0 to nfits-1= 209)
src all srcToPlot,all or
3C274 3C286 3C348 3C405 CasA
f all Freq to plot or all. frqList:
8223 8363 8505 8647 8789 8931 9073
delay secs (0-> no delay)
r 0,1 remove tsys:0,1
s 0,1 single step (0,1)
v min,max .. set vertical scale, blank is autoscale
d debug. stop in routine
q quit monfits
any otherChar continue plotting
currently stepping: 0
v -.2 2 .. change vertical scale -.2 to 2 K
(See /pkg/rsi/local/libao/phil/mas/masmon2dfits.pro)
NAME:
masmon2dplot1 - plot 1 fit from 2d fits
SYNTAX: masmon2dplot1,fiti,ibm,pwr=pwr,lab=lab,remtsys=remtsys
ARGS:
fitI: {} structure returned from mascrossfit2d()
ibm : int index (0..6) for bm to plot in fitI.fitibm[ibm].
KEYWORDS:
pwr : if set then data units are spectrometer pwr. Def: degK
lab : string add to title of plot
var[2]: float set vertical scale. def: autoscale
remTsys: if set then remove constant Tsys value before plotting
DESCRIPTION
This routine is used by masmon2dfits()
(See /pkg/rsi/local/libao/phil/mas/masmon2dplot1.pro)
NAME:
masmonimg - monitor mas files via dynamic spectra.
SYNTAX: masmonimg,bsg=bsg,num=num,projid=projid,date=date,dirI=dirI,$
nspc=nspc,pollist=pollist,desc=desc,zx=zx,zy=zy,noappbm=noappbm
ARGS: none
KEYWORDS:
bsg: string beam,band,group to monitor. Format is bBsSgG where:
B =0..6 is the beam
S =0..1 is the band (0=1450 band, 1=1300 band)
G =0..1 is the group.
default:'b0s0g0'
num: int the file number to start with. (The default is the most recent)
projid: string project id (first part of filename) to monitor. The default
is the most recent one written.
date : long The date (in the filename) to start monitoring. The
default is the current date. Entering "*" will allow any date.
The format is yyyymmdd or "*"
dirI : strarr(2): directory info. If the files are in a non standard place
then pass in location in dirI
dirI[0]=prefix for dir: /share/pdata
Num 1..n goes here
dirI[1]=postfix for dir: /pdev
noappbm: if set then no bm number in directory path.
use for /share/pdata/pdev etc..
IMAGE:
nspc: int the number of spectra for each image. The default is
500 (or the max number in the file). It will try and
round the number to a multiple of spectra per row.
pollist: int pols for image. 1,2,12, -1 --> all pols available
if stokes, it does not plot u,v
zx : int if single pol then zoom (=n) or contract (-n) x axis by this
amount.
zy : int if single pol then zoom (=n) or contract (-n) y axis by this
amount.
RETURNS:
desc : {} You can return from masmon with the current file still
open (using the x command from the menu). The file descriptor
will be returned here (this is actually a structure used by the
masxx routines). It is then your responsibility to close the
file via masclose,desc or masclose,/all.
-----------------------------------------------------------------------------
DESCRIPTION:
Monitor online and offline mas files by making a dynamic spectra of the data. The
files to monitor can be selected by:
1. if no arguments are supplied, use the most recent files from
beam0,band0,group0.
2. Specify the projid,date using the startup keywords.
3. Use the menu inside masmon to select a different set of files to monitor.
- Enter the menu once masmon is started by hitting any key.
-----------------------------------------------------------------------------
EXAMPLES:
masmonimg - start monitoring most recent file for b0s0g0
masmonimg,bsg=b0s1g0 .. start with beam0, band1 , group0
masmonimg,projid='a2130' .. use a2130 files
masmonimg,pollist=1 .. images of only polA
-----------------------------------------------------------------------------
USING THE INTERNAL MENU:
Enter the internal memory by hitting any key.
The menu is:
Cmd CurVal function
--> IMAGES:
p 12 pols for images (1,2,12)
h low high set horizontal range for images (freq low,hihg), blank --> all
nspc 500 number of spectra per image. blank--> auto
c 3 clip to +/- nsig.
delay secs secs to delay between images..
--> FILE POSITIONING:
bsg b0s0g0 beam,band,grp to use (change all at once)
bm 0 beam to plot
band 0 band to display (0 = 1450, 1=1300)
d 20081011 date to use. First number of data displayed
num 00801 number.. file number to display
r 1 position to row. maxRows: 99
l (all) list files this date (all --> all dates)
--> MISC:
hdr display last header read
s 0 step mode. 1-on,0-off
z stop in procedure (for debugging)
q to quit
x quit leaving the current file open. return description in desc keyword
otherKeys continue plotting
Enter the cmd and a possible value. The program will then continue.
The menu options are:
p 12 .. pols for image. By default both pols are used (but in stokes mode
u,v are note used.
.. The pols are numbered: 1-polA,2-polB. Any combination
.. of the pols can be displayed (entered as a one number eg 12..)
Example:
p 12 .. use pols a and polB.
h .. change the frequency range of the image. Units are Mhz. The default
is the entire frequency range.
Example:
h 1400 1450 display frequencies 1400 thru 1450 (if they are in the band)
delay .. secs to delay between plots. Use if the plotting is going
too fast.
Example:
delay .1 Delay .1 seconds between plots
bsg b0s0g0 .. Select the beam,band,group to display.
.. By default beam0,band0,grp0 is displayed.
Example:
bsg b0s1g0 .. will display the band centered at 1300 Mhz.
d 20080826 .. Move to first file on this date using the current
project id.
n 100 .. Move to this filenumber for the current date and project id.
If the number does not exist, you remain on the current file.
l (all) .. list files for this projid, data. If all is included then
list files for this projid for all dates.
r 1 .. position to a row in the current file (count from 1)
hdr .. print last header read
s 0/1 .. step mode will stop and display the menu after each spectra is
plotted. s 1 turns on step mode. s 0 turns it off.
z .. stop in masmon routine (for debugging)
q .. exit masmon.
(See /pkg/rsi/local/libao/phil/mas/masmonimg.pro)
NAME:
masmonstripch - stripchart recording of total power
masmonstripch,bsg=bsg,num=num,projid=projid,date=date,maxpnt=maxpnt,$
rec=rec,v1=v1,$
rdstep=rdstep,pltDelay=pltDelay,sub=sub,win=win,polList=polList,$
median=median,maskar=maskar,colar=colar,noappbm=noappbm,tming=tming,dirI=dirI
-> Warning.. i kludged this documentation from masmon and an older version of masstripch
I need to verify that i've updated the documentation correctly
ARGS:
KEYWORDS:
bsg: string beam,band,group to monitor. Format is bBsSgG where:
B =0..6 is the beam
S =0..1 is the band (0=1450 band, 1=1300 band)
G =0..1 is the group.
default:'b0s0g0'
num: int the file number to start with. (The default is the most recent)
projid: string project id (first part of filename) to monitor. The default
is the most recent one written.
date : long The date (in the filename) to start monitoring. The
default is the current date. Entering "*" will allow any date.
The format is yyyymmdd or "*".
polList:int pols to use. 1,12,-1 (all pols a,b def)
maxpnt: long max number of points to display at once. default: 3600
v1[2] : float min,max value for top plot
sub : if set then subtract running mean
v1[2] : float min,max value for top plot
sub : if set then subtract running mean
pltDelay: float number of seconds to wait after each plot. Use this
to slowly scan a file that already exists.
win : long window number to plot in. Default is 0
median : if set, use the median to compute the power
maskar[nchn]:int if provided then use to mask off chan
with rfi. use chan where maskAr[] > 0
colar[2] :int color indices to use for plotting polA,B
def: 1-white,2-red
1-black,2-red,3-green,4=blue,5=yellow,6=
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 2 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).
v1[2] : float min,max value for top plot
sub : if set then subtract running mean
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.
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:
Cmd CurVal function
--> PLOTTING:
p 1234 pols to plot (1..4)
h low high set horizontal plot range ..no args --> autoscale
v low high set vertical plot range ..no args --> autoscale
delay secs secs to delay between plots..
--> POSITIONING:
bsg b0s0g0 beam,band,grp to use (change all at once)
bm 0 beam to plot
band 0 band to display (0 = 1450, 1=1300)
d 20081011 date to use. First number of data displayed
num 00801 number.. file number to display
r 1 position to row. maxRows: 99
l (all) list files this date (all --> all dates)
--> MISC:
hdr display last header read
s 0 step mode. 1-on,0-off
z stop in procedure (for debugging)
q to quit
x quit leaving the current file open. return description in desc keyword
otherKeys continue plotting
EXAMPLES:
1. monitor the online datataking:
masmonstripch,bsg='b0s1g0',proj=proj
2. set fewer points per plot for better resolution. Set the top vertical
scale to .8,2.
masmonstripch,desc,maxpnt=600,v1=[.8,2]
3. If you want to position to a different filenumber:
hit any character
num newFileNumbe r
4. If you want to monitor a file offline, with no wait between
updates, and 500 points plotted:
masmonstripch,bsg=bsg,proj=proj,maxpnt=500,v1=[.5,3]
5. Do the same thing but wait 1 second between plots and read 200 points at
a time:
masmonstripch,bsg=bsg,proj=proj,lun,maxpnt=200,v1=[.5,3],pltDelay=1
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/mas/masmonstripch.pro)
NAME:
masmostrecentfile - find most recent mas file
SYNTAX:stat=masmostrecentfile(projId,ldate,bsg,lnum,baseDir,newFI,$
curFI=curFI,flist=flist,oldest=oldest,dirI=dirI,noappbm=noappbm,$
src=src,psrfits=psrfits)
ARGS:
projid: string project to search for. "*" for any project. Project
if the first part of the filename (projId.yyyymmdd.xxx)
ldate : string date to search for. Format is: yyyymmdd. Use "*" for any
date.
bsg: string BeamSubbandGroup to search for. eg:
b0s0g0 for beam0, subband0, group 0.
lnum: string file number to search for. Should be 5 digits eg: 00100.
Use "*" for any file number.
KEYWORDS:
curFI:{} file Info structure. If supplied then file the next file after
this file.
oldest:int if set then return the oldest rather than the newest file.
Use this to find the first file of a day.
dirI:strarr[2] If supplied then this specifies the prefix and suffix for the
directories to search. The prefix is the directory name prior to
to directory number and the suffix is the directory path
following the number. The default is:
dirI[0]='/share/pdata', dirI[1]="/pdev/"
noappbm: if set then dirI[0],dirI[1] do not have a beam number
between them. use for data in /share/pdata/pdev
src: '' if supplied then require this source name (for psrfits)
psrfits: if set then search for psrfits files that have a srcname
RETURNS:
stat: 1: found file, -1 no file.
basdir: string base dir for files we searched for: eg
/share/pdata1/pdev/
newFI: {} File info structure for file that was found.
flist:strarr list of files found from the ls command.
DESCRIPTION:
Find the most recent file that fits the specifications provided by the
user. This routine is normally used for monitoring where you want to
find the next file automatically. The routine uses ls filename where
filename may contain some wildcards.
Examples for using it.
1. If monitoring a days worth of files:
a. set projid,ldate,bsg, lnum=*, and set oldest
on return set curFI=newFI
b. When ready to move to the next file, pass in curFI=curFI to find the
next file.
2. for online montoring of most recent file do the same as 1. do not set oldest.
(See /pkg/rsi/local/libao/phil/mas/masmostrecentfile.pro)
NAME:
masonofffind - find onoff position switch patterns
SYNTAX: n=masonofffind(projId,
patI,npat,upatIdAr,nsrc,usrcAr,$
yymmdd=yymmdd,appbm=appbm,dirI=dirI,bm=bm,band=band,grp=grp)
ARGS:
projid: string project id to search for
KEYWORDS:
yymmdd: long yyyymmdd limit to this date
appbm: apply bm number to each directory. This should normally
be set:
dirI[2]: string if data not in default /share/pdataN/pdev
see masfilelist for usage.
bm : int if supplied limit to this beam
band : int 0,1 limit to this band. Note if you have single pixel
data you probably should set band=1 or masfilelist may
not find the files.
grp : int limit to this group 0, or 1.
RETURNS:
n :int number proccessed patterns. Each pointing pattern
can generate multiple processed patterns if you use
multple bms.
patI[n]: info on processed patterns
npat : int number of pointing patterns executed.
nsrc : int number of unique sources
srcAr[nsrc]: string list of sources used.
patI STRUCTURE:
- there will be one entry for each pattern and bm used
If you have 4 bms then a single patId will have 4 entries
in patI[]
patI[i].srcNm - source name
patI[i].nscans - # of scans in pattern ..
patI[i].bm - bm number 0..6
patI[i].grp - grp 0 or 1
patI[i].flist[maxScans] - list of filenames:
same orde as scanTypeAr
patI[i].sumI[maxScans] - summary info for each scan
DESCRIPTION:
Find all of the onoff patterns the given project and constraints.
The routine accepts the same parameters at masfilelist:
yymmdd=,dirI=,bm=,band=,grp=num=,appbm=appbm
The routine will pass back an array of patI structs holding info on
each pattern.
There is a distinction between a pointing pattern and a processed
dataset. You can have multiple processed datasets for each pointing
pattern if you used multiple beams or groups.
EXAMPLE:
1. Find the onoff info for project a2516 on 20100930 group 0.
n=masonofffind('a2516',patI,npat,upatIdAr,nsrc,usrcAr,$
yymmdd=20100930,grp=0,band=1,/appbm)
patI[n] - holds the info
npat - number of pointing patterns found
upatIdAr[npat] - unique pattern id for each pointing pattern.
This can include multiple patI[] entries if multiple
bms taken on each pointing pattern.
nsrc - Number of unique on source names found
usrcAr[nsrc] - list of source names.
2. now process all of the patterns found using group 0 bm 1.
ii=where((patI.grp eq 0 ) and (patI.bm eq 0),cnt)
for i=0,cnt-1 do begin
j=ii[i] ; index in patI for this dataset
istat=masposonoff(patI[j].filelist,b,...params for masposonoff)
if i eq 0 then b0ar=replicate(b[0],cnt); generate large array
b0ar[i]=b
endfor
this assumes that all of the bm=0 data sets have the same
number of channels (since we are trying to store them all in an array).
3. plot out the first set of results
masplot,b0ar[0]
(See /pkg/rsi/local/libao/phil/mas/masonofffind.pro)
NAME:
masopen - open mas file for reading
SYNTAX: istat=masopen(filename,desc,fnmI=fnmI,hdr=hdr)
ARGS:
filename: string filename to open (unless fnmI is specified)
fnmI : {] returned from masfilelist. if provided,
then ignore filename and use this as the
file to open.
KEYWORDS:
RETURNS:
istat: 0 ok
-1 could not open file.
desc : {} file descriptor to pass to the i/o routines.
hdr : [string] if present then return fits bintable header for the fille
(See /pkg/rsi/local/libao/phil/mas/masopen.pro)
NAME:
masparsfnm - parse a mas file name
SYNTAX: istat=masparsfnm(filename,fnmI)
ARGS:
filename: string mas filename to parse
RETURNS:
istat: int 1 if valid mas name, 0 if not a valid mas name
fnmI : {} structure holding the parsed information.
The format is:
dir/proj.date.obs.brdId.seqnum.mas
where obs is optional
EXAMPLE:
filename='/share/pdata/mas/agc110443/x107.20070123.agc110443.b0a.00000.mas'
istat=masparsfnm(filename,fnmI)
IDL> help,fnmI,/st
* Structure PDEVFNMPARS, 8 tags, length=60, data length=60:
DIR STRING '/share/pdata1/pdev/'
FNAME STRING 'if2noise.20081107.b0s0g0.00000.fits'
PROJ STRING 'if2noise'
DATE LONG 20081107
BM INT 0
BAND INT 0
GRP INT 0
NUM LONG 0
(See /pkg/rsi/local/libao/phil/mas/masparsfnm.pro)
NAME:
maspatfind - find data taking patterns
SYNTAX: n=maspatfind(projId,obsModeAr,scanTypeAr,minScanPat,recCntMatch,
patI,npat,upatIdAr,nsrc,usrcAr,$
yymmdd=yymmdd,appbm=appbm,dirI=dirI,bm=bm,band=band,grp=grp)
ARGS:
projid: string project id to search for
obsModeAr[n]: strarr obsmode for each step in pattern
scanTypeAr[n] : strarr scan type name for each scan in pattern
minScanPat : int minimum number scans required for pattern.
obsmodear,scantypear should have the optional
scans at the end. eg for on,off you could optionally
have a cal.
recCntMatch[n]: int index in obsModeAr that has to have the same
number of records.
eg: suppose onoff position swithch with cals.. then
recntMatch=[-1,0,-1,2] index from 0, -1 --> don't care
[0]=-1 pos on don't care any number is on .. on
[1]=0 pos off has to match cnt in on.
[2]=-1 cal on don't care
[3]=2 cal off has to match cnt in cal on
KEYWORDS:
yymmdd: long yyyymmdd limit to this date
appbm: apply bm number to each directory. This should normally
be set. If yymdd=0 then any date will match.
dirI[2]: string if data not in default /share/pdataN/pdev
see masfilelist for usage.
bm : int if supplied limit to this beam
band : int 0,1 limit to this band. Note if you have single pixel
data you probably should set band=1 or masfilelist may
not find the files.
grp : int limit to this group 0, or 1.
RETURNS:
n :int number proccessed patterns. Each pointing pattern
can generate multiple processed patterns if you use
multple bms.
patI[n]: info on processed patterns
npat : int number of pointing patterns executed.
nsrc : int number of unique sources
srcAr[nsrc]: string list of sources used.
patI STRUCTURE:
- there will be one entry for each pattern and bm used
If you have 4 bms then a single patId will have 4 entries
in patI[]
patI[i].srcNm - source name
patI[i].nscans - # of scans in pattern ..
patI[i].bm - bm number 0..6
patI[i].grp - grp 0 or 1
patI[i].flist[maxScans] - list of filenames:
same orde as scanTypeAr
patI[i].sumI[maxScans] - summary info for each scan
DESCRIPTION:
Find all of the patterns the given project and constraints.
The routine accepts the same parameters at masfilelist:
yymmdd=,dirI=,bm=,band=,grp=num=,appbm=appbm
The routine will pass back an array of patI structs holding info on
each pattern.
There is a distinction between a pointing pattern and a processed
dataset. You can have multiple processed datasets for each pointing
pattern if you used multiple beams or groups.
EXAMPLE:
1. Find the dps info for project a2516 on 20100930 group 0.
n=masdpsfind('a2516',patI,npat,upatIdAr,nsrc,usrcAr,$
yymmdd=20100930,grp=0,band=1,/appbm)
patI[n] - holds the info
npat - number of pointing patterns found
upatIdAr[npat] - unique pattern id for each pointing pattern.
This can include multiple patI[] entries if multiple
bms taken on each pointing pattern.
nsrc - Number of unique on source names found
usrcAr[nsrc] - list of source names.
2. now process all of the patterns found using group 0 bm 1.
ii=where((patI.grp eq 0 ) and (patI.bm eq 0),cnt)
for i=0,cnt-1 do begin
j=ii[i] ; index in patI for this dataset
istat=masdpsp(patI[j].flist,b)
if i eq 0 then b0ar=replicate(b[0],cnt); generate large array
b0ar[i]=b
endfor
this assumes that all of the bm=0 data sets have the same
number of channels (since we are trying to store them all in an array).
3. plot out the first set of results
masplot,b0ar[0]
BUGS:
1. this routine assumes a scan (for 1 bm) always fits in 1 file.
(See /pkg/rsi/local/libao/phil/mas/maspatfind.pro)
NAME:
maspatidday - get the pattern id's for a day and obsmode
SYNTAX: npatid=maspatidday(yymmdd,proj,obsmode,patIday,src=src,bmtotest=bmtotest,$
patid=patid)
ARGS:
yymmdd long date to test
proj string project id to test
obsmode string the pattern must include this obsmode type:see list below
(it can contain additional obsmodes eg ONOFF with CAL
The exception is if patid= keyword is used.
KEYWORDS:
src: string if supplied then the pattern must also be on this source.
def: any source
bmToTest: int which bm (mock box to check) 0..6
patid : long if present then ignore obsmode.
return info for patid.. If src= is supplied then this patid
must have src as its object..
RETURNS:
npatid int number of unique patterns we found
patIdDay: {} struct holding the info
Search for all patterns for date yymmdd and project=proj.
To be included the pattern must include the obsmode=obsmode (see below)
If the keyword src is included then the pattern must also have src for the
source observered.
The routine only scans a single mock bm (by default bm0). the keyword bmtotest=
will use that beam instead (0..6).
If the patid= keyword is present, then the routine ignores obsmode
and searches for the single patid=patid (in bmToTest). If src= is supplie then
this pattern id must have src as its object.
The info is returned in the structure patIday
help,patIday,/st
NPAT LONG 2 number of patterns found
PROJ STRING 'sun' project id used
SRC STRING 'sun' source requested
PATI STRUCT Array[2] holds the info on each pattern
the patI struct holds:
help,patIday.patI,/st
PATID LONG 211700002 patid for this pattern
FNUM0 INT 0 first filenumber for pattern
FNUM1 INT 2800 last filenumber for pattern
HDR0 STRUCT MASFHDR Array[1] fits header first scan, bm=bmToTest
HDR1 STRUCT MASFHDR Array[1] fits header last scan
obsmode values are:
ONOFF - on off position switching
RADECMAP - map driving in ra, stepping in dec
ON - just tracking a source
CAL - calon,off will probably be part of another pattern. eg, onoff
CROSS - az,el cross on a source
Note:
this routine will not work with pulsar observations since the psrfits files
have a different header
(See /pkg/rsi/local/libao/phil/mas/maspatidday.pro)
NAME:
masplot - plot mas spectra
SYNTAX: masplot,b,freq=freq ,over=over,pollist=pollist,norm=normRange,median=median,$
off=off,chn=chn,smo=smo,retvel=retvel,restFreq=restFreq,velCrdsys=velCrdSys,$
mfreq=mfreq,colar=colar,edgefr=edgefr,flip=flip
ARGS:
b[n]: {b} mas structure from masget to plot
KEYWORDS:
freq[n]: float if provided then use this as the x axis
over: if set then overplot spectra
pollist: long pols to plot: 12 = polA,B, -1--> all pols available
norm[2]: float normalize the spectra. If two element array then is is the range
of frequency bins to use for the mean (cnt from 0). If a single
element or /norm then use all the bins. If keyword median set then
use mean rather than mean
median: if set then use median rather than mean when normalizing
off: float If multiple spectra then increment each spectra by off.
chn : if set then plot vs channel number (count from 0)
smo : int smooth by this many channels
** to plot versus velocity
retvel: if set (/retvel) then plot versus velocity (see masfreq())
restFreq: float Mhz. if /retvel then use restFreq as the restFreq when
computing the velocity. The default is to use the hdr restFreq.
velCrdSys:string If supplied then use this as the velocity coord system. The
default is to use what is in the header. The values are:
"T":topcentric,"G":geocentric,"B":Barycenter,"L",lsr
mfreq : multiple frequencies.if set then the b[n] array has different
freqs (single pixel).
It we replot a new x axis each time. User needs to
setup the hor value ahead of time.
Note: if you also select /retvel, you'd better supply
restFreq.
colar[4]: int color indices for pol1,2,3,4.1-black,2-red,3-green,4-blue
5-yell w..
edgefr : float if provided, the fraction of band on each edge to not plot
flip : if set then flip the freq array. use when hdr info is wrong.
DESCRIPTION:
Plot the various spectr in b. If b.accum > 0 the scale by one over this
when plotting (but don't change the data).
(See /pkg/rsi/local/libao/phil/mas/masplot.pro)
NAME:
masplotmb - plot mas multi beam data
SYNTAX : masplotmb,b,bmlist=bmlist,vel=vel,rest=rest,off=ploff,pol=pol,
over=over,nolab=nolab,extitle=extitle,newtitle=newtitle,
col=col,sym=sym,chn=chn,$
xtitle=xtitle,ytitle=ytitle,flag=flag,lnsflag=lnsflag
ARGS:
b[nspc,nbms]: {masget} data to plot
KEYWORDS:
bmlist: int This is the beams to plot numbering beams 0..6.
To plot beams 0,2,4,6 use 2360 (note that beam 0 can't come
first). By default all beams are plotted.
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)
norm : if set then normalize the spectra for a given pol,beam before
plotting..
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
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
DESCRIPTION:
masplotmb will plot the mas spectra. The data should come from
masavgmb(). If one dimension, then the dimension is over beams. If
two dimensions then b[nspc,nbeams].
You can plot:
- any combination of sbc using pol=pol and bmlist
- by topocentric frequency (defaul), rest frequency, or by velocity
- make a strip plot if b is an array by using off=ploff.
EXAMPLES:
masavgmb(yymmdd,projid,filenum,b1,b2,...)
masplotmb,b1 plot all sbc
masplotmb,b,brd=0 plot brd0 (beam 0)
masplotmb,b,brd=1 plot brd1 (beam 1)
masplotmb,b,brd=2,pol=2 plot brd2 (beam 2),polB only
masplotmb,b,brd=3,pol=1 plot brd3 (beam 3) pola only
masplotmb,b,brd=30,/vel plot brd0,3 (beam0,3) by velocity
masplotmb,b[1],brd=0,/over,col=[5,5]
NOTE:
use hor,ver to set the plotting scale.
oveplotting only works if you display 1 brd at a time.
SEE ALSO:
masavgmb. hor,ver in the generic idl routines.
(See /pkg/rsi/local/libao/phil/mas/masplotmb.pro)
NAME:
maspos - position to row in fits file
SYNTAX: istat=masposrow(desc,row)
ARGS:
desc:{} returned from masopen().
row : long row to position to. Count rows from 1
RETURNS:
istat: 0 ok
(See /pkg/rsi/local/libao/phil/mas/maspos.pro)
NAME:
masposonoff - process a position switch scan
SYNTAX: istat=masposonoff(flist,b,bonrecs,boffrecs,sclcal=sclcal,
sclJy=sclJy,median=median,double=double,smooff=smooff,$
sclmask=sclmask,calI=calI,edgefract=edgefract,nocheck=nocheck,remdc=remdc)
ARGS :
flist[4,n]: string filenames for on,off,calon,caloff
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.
han : if set then hanning smooth the data. not yet supported..
median : if set then median filter rather than average
the poson,posoff,calon,caloff records.
double : if set then return doubles. default is float.
smooff : if supplied then smooth the off by this number of channels
before dividing.
-1,0,1 no smoothing
> 1 smooth by averaging
< -1 median fitler by this many channels.
dividing.
sclmask:int used when computing scale factor that average frequency.
These averages include:
avg(tcal[freq]),avg(calOn[freq]-caloff[freq]) and
1./avg(posOff[freq])
0 - (default)
linear fit to calon/caloff then all points 3 sigma from
the fit are flagged as bad. These points will
not be included in the averages.
1 - Use sclmask=0 (calon/caloff) plus check the posOn,posOff
data.
:compute rms/mean along each channel for the posOn
and posoff datasets then avgerage the two rms's.
:robust linear fit to rms throwing out freq channels
of more than 3 sigma.
2 - exclude 6% of channels on each side of the bandpass.
All other channels are used for the averages.
The result of sclMask is returned in calI.indused. These are the
indices into the freq array that were used for the averages.
edgefrac[2]:float fraction of edge on each side to exclude when
computing cal scaling.
nocheck : if set then don't bother to check the type of files.
remdc : if set then interpolate across the dc channel in the middle of the spectra
RETURNS:
istat int 1: completed ok.
0: did not complete
-1: error...
b: {masget} return (on-off)/off * X here..
hX is:
sclCal true : (kelvins/mascnt)/.. units are then K
sclcal false: off .. units are TsysOff
bonrecs[] {masget} return individual on records. This is an
optional argument. The units will be K (if /sclcal),
Jy if /sclJy, else masdevice units.
boffrecs[] {masget} return individual off records. This is an
optional argument. The units will be K (if /sclcal),
Jy if /sclJy, else masdevice units.
keywords that return info:
calI[] {calI} return cal info. see mascalonoff..
DESCRIPTION:
Process a position switch on,off pair. Return a single spectrum
per pol 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 header will be that from the first record of the on with the
following b.h.azimuth, b.h.elevatio set to the average for the on scan
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 masdevice units (linear in power).
structures:
b - holds the spectral data and headers:
calI: A structure holding info on the cal,caloff
calI.CALVAL[2] : Cal values used for polA,polB
calI.exposure : integration time that cntstok corresponds to.
calI.NUMPOL : number of pols. normally 2. 1--> stokes I
calI.CNTSTOK[2]: Conversion factors masDevice cnts to Kelvins
For pola, polB.
calI.npnts : number of points in spectra used for averaging.
calI.flipped : 1 normal, -1 spectra was flipped
calI. EDGEFRACT[2]: edgeFract[0,1] excluded
calI.USEMASK : 1--> use mask, 0--> use edge fraction to throw
out outliers.
calI.INDUSED[calI.npnts]: indices into spectra used to compute
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
3. If the acquired band is larger than the bandwidth of the receiver
( eg: 800 Mhz rcvr is about 100 Mhz, and you might use a 172 Mhz
band to measure it)
then you should use the edg=keyword to limit the computation to the
part of the band that has signal.
example: suppose there are 8192 channels.
masplot,b,/chn
look at where the band rises and falls. Suppose it rises at channel
1800 and falls at channel 6800.. then
edge=[1800,8192.-6800]/8192
(See /pkg/rsi/local/libao/phil/mas/masposonoff.pro)
NAME:
masrms - compute rms/Mean by channel
SYNTAX: brms=masrms(bin,rembase=rembase,nodiv=nodiv)
ARGS:
bin[] : {masget} array of masget structures
KEYWORDS:
rembase : If set then remove a linear baseline (by channel) before
computing the rms.
nodiv : If set then don't divide by the mean values of each chan.
RETURNS:
brms : {masget} where brms.bN.d[] is now the rms/mean by channel.
DESCRIPTION
The program will compute the rms/mean for each frequency channel. If the
/nodiv keyword is set then the rms is not divided by the channel mean (so
you can then see the size of the rms in spectrometer counts). The rembase
keyword will remove a linear baseline in each channel before computing
the rms.
The input data (bin) can be a single row element with multiple dumps per
row or it can be an array of row elements (with one or more dumps per row).
There must be at least 5 number per channel for the rms computation.
(See /pkg/rsi/local/libao/phil/mas/masrms.pro)
NAME:
masstokes - intensity calibrate stokes data.
SYNTAX: istat=masstokes(bdatIn,bcalIn,descDat,descCal,$
bdatOut,bcalOut,
avg=avg,han=han,edgefract=edgefract,
mask=mask,cmpmask=cmpmask, phase=phase
bpc=bpc,fitbpc=fitbpc,smobpc=smobpc,phsmo=phsmo,$
nochk=nochk,mmcor=mmcor,mmret=mmret
ARGS:
bdatIn[n]:{} data to calibrate.
bcalIn[2]:{} cal on,off. Should already be averaged to 1 on rec,
1 off rec.
descDat :{} descriptor for datafile
descCal :{} descriptor for cal file
These are needed to check the ashift pdev
parameter for the scaling of the data. It is stored
in desc.hsp1.
KEYWORDS:
avg: if set then return the averaged data record
han: if set then hanning smooth the data.
edgefract: float fraction of bandpass on each side to not use during
calibration. default .1
mask[nchan]:int if supplied then mask to use when computing things
0--> do not include channel. nchan must match
number of channels in spectra.
cmpmask : if set then compute mask from calon/caloff
and ignoring 6% of channels on each edge.
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 (1 portions of the mask excluding
the outside edges where the filter falls off.
This is only used if bpc =1,or 2.
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 this dataset.
RETURNS:
bdatOut: {masget} intensity calibrated data spectra
bcalOut: {masget} 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
:-7 desc.ashifts0s1 or ashifts2s3 not equal
DESCRIPTION:
masstokes will intensity calibrate (and optionally phase calibrate)
stokes data given data and a cal on,off. The data is passed in via
bdatInp. The cal is input vis bcalInp. The cal should already be averaged
so there is 1 cal on spectra and 1 cal off spectra.
The descDat,descCal are the descriptors from masopen() of the data
and the cal file. They are needed to get the pdev scaling done
on the data (desc.hsp1).
On output bdatOut and bcalOut will be in units of Kelvins.
The /avg keyword will cause the bdatOut to be averaged to a single
spectra.
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.
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:
I (Px*cosft(XX) + Py*cosft(YY))/(Px+Py)
Q (Px*cosft(XX) - Py*cosft(YY))/(Px+Py)
U (real(fft(cacf)*Px*Py)/(Px+Py)
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. Compute the average <calon>,<calOff> over the specified channels.
The specified channels are determined by:
a. The mask 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
I = (XX*scaleXX + YY*scaleYY)
Q = (XX*scaleXX - YY*scaleYY)
U = U*sqrt(scaleXX*scaleYY)
V = V*sqrt(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:
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,masstokes(lun,bdat,bcal); will process the first set
--print,masstokes(lun,bdat,bcal,/han); will process the 2nd set with hanning
To process the 2nd set directly with an edgefraction=.12:
--print,masstokes(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,masstokes(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,masstokes(lun,bdat1,bcal1,calinp=bcalinp,scan=210200235L)
--print,masstokes(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,masstokes(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:
WARNING:
01dec09 still in progress
.. mmcor,mmret do not yet work so don't try them
(See /pkg/rsi/local/libao/phil/mas/masstokes.pro)
NAME:
masstripch - stripchart recording of total power
SYNTAX: masstripch,desc,maxpnt=maxpnt,v1=v1,v2=v2,sub=sub,$
rec=rec,rdDelay=rdDelay,pltDelay=pltDelay,$
rdStep=rdStep,win=win,median=median,maskar=maskar,$
colar=colar
ARGS:
desc: {} from masopen()
KEYWORDS:
polList:int pols to use. 1,12,-1 (all pols a,b def)
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)
sub : if set then subtract running mean
scan : long position to scan before starting. The
default is the current position
rec : long record of file to position to before starting.
cnt from 1..def: 1
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.
win : long window number to plot in. Default is 0
median : if set, use the median to compute the power
maskar[nchn]:int if provided then use to mask off chan
with rfi. use chan where maskAr[] > 0
colar[2] :int color indices to use for plotting polA,B
def: 1-white,2-red
1-black,2-red,3-green,4=blue,5=yellow,6=
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:
istat=masopen(file,desc)
masstripch,desc
2. set fewer points per plot for better resolution. Set the top vertical
scale to .8,2. and the bottom plot to -.4,.4.
masstripch,desc,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:
istat=masopen(file,desc)
masstripch,desc,maxpnt=500,v1=[.5,3]
5. Do the same thing but wait 1 second between plots and read 200 points at
a time:
masstripch,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/mas/masstripch.pro)
NAME:
mastdim - decode fits tdim1 keyword
SYNTAX: mastdim,h,nchan=nchan,nra=nra,ndec=ndec,npol=npol,nspc=nspc
ARGS:
h{} : struct fits header (b.h)
RETURNS:
istat: 0 ok
-1 could not parse tdim1
nchan: long number of frequency channels
nra: long number of ra points
ndec: long number of dec points
npol: long number of pols 1,2,4
nspc: long number of specra
(See /pkg/rsi/local/libao/phil/mas/mastdim.pro)
NAME:
mastsyshwcal - compute tsys from files with hardware winking cal
SYNTAX: n=mastsyshwcal(fnmAr,tsysI,fnmIAr=fnmIAr,verb=verb,$
edgefract=edgefract,rmsdigar=rmsdigar)
ARGS:
fnmAr[n]: strarr list of files to process
KEYWORDS:
fnmIAr[n] :{} files to input. returned from masfilelist. If
present, ignore fnmAr[]
edgefract[2]:fltarr fraction of bandpass on each edge to ignore when computing
tsys. default [.06,.06]. if single value entered then use
same fraction on each edge.
verb : int if set then output file info as they are processed
RETURNS:
n : int > 0 number of entries in tsysI[n]. it will skip any files that do not
include the winking cal.
< 0 error
tsysI[n]: struct arr holding tsys info for each file
rmsdigar[2,2,n]: fltarr return dig rms [i/q,POlA/B,nfiles]
DESCRIPTION:
Compute tsys from files with hardware winking cal data. Each file
in fnmAr[n] (or fnmIAr[]) is input using masgethwcal(). The
cal value is input and then tsys is computed using the cal off specdtra. Tsys info
is returned in the tsysI[] structure. It contains:
YYMMDD LONG 20220711 ; ast date start of file
SECMID LONG 84093 ; ast secmid start of file
SEC1970 DOUBLE 1.6575961e+09 ; secs 1970 start of file
H STRUCT -> MASFHDR Array[1]; mas header
SCAN LONG 219200151 ; scan number
BM LONG 0 ; bm/band 0..6
CFR FLOAT 8219.00 ; center freq of band in MHz
AZ FLOAT 50.7062 ; az start of file
EL FLOAT 53.9142 ; el start of file
TSYSA FLOAT 103.902 ; tsys polA deg K avg over band
TSYSB FLOAT 97.4181 ; tsys polB deg K (if 2 pol)
CALI STRUCT-> Array[1] ; cal info, includes how much of band was used
to compute tsys (edge, rfi ignored).
EXAMPLE:
1. process tsys for a given day for proj sscontin
proj='sscontin'
yymmdd=20220711
nfiles=masfilelist('',fnmIAr,proj=proj,yymmdd=yymmdd,band=1,grp=0,/appbm)
nTsys=mastsyshwcal('',tsysI,verb=verb,fnmIar=fnmIAr)
ver, 80,125
plot,tsysI.el,tsysI.tsysA
oplot,tsysI.el,tsys.tsysB,col=colph[2]
(See /pkg/rsi/local/libao/phil/mas/mastsyshwcal.pro)