NAME:
acormapexamples - Examples using correlator mapping routines:
Starting idl.
- To use the correlator mapping routines you need to tell idl where to
look to get these procedures. You can do it manually each
time you run idl, or you can put it in an id startup file.
Manually:
idl
@phil
@corinit
Using a idl setup file:
Suppose your home directory is ~jones.
create the file ~jones/.idlstartup
add the line
!path = '/home/phil/idl/gen:' + !path
to this file.
In your .cshrc file (if you run csh) add the line
setenv IDL_STARTUP ~/.idlstartup
You can then run idl with :
idl
@corinit
You can also put any other commands in the startup
file that you want to be executed each time you
start idl. My startup file contains:
!EDIT_INPUT=500
@phil!path = '/home/phil/idl/gen:' + !path
1.Which data can use these routines.
You should be able to process data taken with cormap1, cordrift, or
cormapdec with these routines.
2. Cookbook processing of a map:
A. inputing the spectra:
- If the entire map is in a single file used:
openr,lun,'/share/olcor/corfile.17dec00.a1397.1',/get_lun ; open the file
scan=35253227L ; first scan of the map
istat=cormapinp(lun,scan,1,2,m,cals) ; input the map from brds 1,2
free_lun,lun
- If the map is spread over many files use:
stat=cormapinplist( )
- You can also input a map using the data archive. see
arch_getmap in the correlator section.
m[2,nsmp,nstrips].d is in correlator power units
B. Scale to Kelvins
istat=cormapsclk(m,cals,tsys) ; scale to kelvins using the cals
.. m[2,nsmp,nstrips].d is in kelvins
C1. baseline the data.
mb=cormapbl(m,mask,polyDeg) ; baseline the map
C2. bandpass correct the data:
numbpedge=2
mb=cormapbc(m,numbpedge)
.. This will:
1. Create an average bandpass for each strip. It can come from:
- bandpasses on the edge of each strip,
- averageing the entire strip.
2. For each spectra of the strip compute
spc[i]=spc[i]/normalized(spcBpc)
3. Remove Tsys from each spectra by:
- interpolating by the spectra used for the band pass correction
- interpolate using all of the spectra in a strip
- remove Tsys from each spectra by computing the mean of each
spectra
D. gain correct the data
cmgaincor(mb,gaincor) ; gain correct the map
.. mb.d is now gain corrected. gaincor[] holds the correction applied.
E. Things to watch out for..
- cmgaincor assumes the data is in nongaincorrected Kelvins. Once
you've called it it is in true kelvins (or Jy). Don't recall it
with the new data..
3. The data structures.
.. map data structure input by cormapinp
help,m
M STRUCT = -> Array[2, 15, 15]
[npol,nsamplesperstrip,nstrips]
help,m,/st
H STRUCT -> HDR Array[1]. header (see cor documentation)
P FLOAT 1.08431 ..total power
AZ FLOAT 284.581 ..azimuth pos deg at end of sample
ZA FLOAT 18.9583 ..za pos deg at end of sample
AZERRASEC FLOAT 0.256720 ..az tracking error asecs
ZAERRASEC FLOAT 0.289800 .. za tracking error asecs
RAHR FLOAT 3.26740 ..RA in hours,middle of sample
DECDEG FLOAT 12.5167 ..DEC in deg ,middle of sample
D FLOAT Array[2048] ..spectral data. correlator units.
.. cals data structure input by cormapinp:
help,cals1,/st
H STRUCT -> HDR Array[1] .. header for cal on record
CALVAL FLOAT Array[2] .. pola,b cal value in Kelvins
CALSCL FLOAT Array[2] .. pola,b. scale correlator
counts to kelvins.
4. Displaying things:
.. average the bandpasses for the entire map
avgBpPolA=cmavgstrips(mbl,1)
avgBpPolB=cmavgstrips(mbl,2)
plot,avgBpPolA
oplot,avgBpPolB
.. average all spectra strips 3-10 counting from 1.
avgBpPolA=cmavgstrips(mbl,1,first=3,last=10)
avgBpPolB=cmavgstrips(mbl,2,first=2,last=10)
.. display whole map overplotting all spectra in each strip
ver,35,65
cmstrips,m,1,0,2 ;offset each strip by 2 kelvin
cmstrips,m,1,0,2,/freq ;plot versus frequency
cmstrips,m,1,0,2,/vel ;plot versus velocity
.. plot individual spectra.
plot,m[0,1,3].d ; pola,2nd sample,4th strip
.. average over a range of frequency channels then display as an
image.
pol=1 ; do pol a
chn1=500 ; first channel, count from 1
chn2=550 ; last channel
avgchn=cmavgchn(mbl,pol,chn1,chn2); avgchn[nsmp,nstrips]
xloadct ; start color widget, clip on grey scale
imgdisp,avgchn,zx=10,zy=10,/histeq ; image display, zoom 10 by 10
..Do an animation of the map.
..assume a 15,15 map with 2048 frequency channels.
..display channels 951 thru 1100 (150) expanding the ra,dec dimensions
to also be 150,150
d=reform(mbl[0,*,*].d[951:1100],150,15,15); grab the data
dh=imghisteq(rebin(d,150,150,150)) ;interpolate,histogram equalize
xinteranimate,set=[150,150,150],/showload
for i=0,149 do xinteranimate,frame=i,image=reform(dh[i,*,*],150,150)
xinteranimate,/keep_pixmaps
.. click on the colors button to setup the colormap you want
.. click on the active slider to see which frame is displayed.
.. change the display rate with the top slider (frames/sec).
5. Notes on scaling the map to kelvins.
istat=cormapsclk(m,cals,tsys)
.. m[2,n,m].d will now be in kelvins.
.. tsys[2,n,m] will be the system temperature measured at each spectra.
..
.. You should examine the cals data to make sure that there are not
.. any bad cal samples. cals.calscl[0] should be a relatively constant
.. value since it is converting from correlator counts to kelvins
.. using (calInKelvins)/(calon-caloff). Variation in these values
.. can come from the cal temperature changing with time, the gain
.. of the if/lo changing, or rfi in a cal on and not in the off.
.. The cals are temperature controlled. This value does change by
.. about 1% per degree F change in the rotary floor room.
.. This is probably some amplifier gains changing.
.. look at the values with:
plot,cals.calscl[0] .. polA
plot,cals.calscl[1] .. polB
.. If 2 cals/strip were taken then everyother strip has also
.. had it's cals flipped (along with the data samples).
.. Look at the system temperatures with:
plot,tsys[0,*,*] polA all strips
(See /pkg/rsi/local/libao/phil/Cor2/cormap/acormapexamples.pro)
NAME:
cmavgchn - average freq channels of a map
SYNTAX: avgchn=cmavgchn(m,pol,chn1,ch2)
ARGS:
m[2,nsmp,nstrips] : map structure array
pol : 1 or 2.. pol to return.
chn1 : first freq chn for average. count from 1
chn2 : last freq chn for average. count from 1
RETURNS:
avgchn[nsmp,nstrips: float the averaged channels.
DESCRIPTION:
Average frequency channels chn1 through chn2 of the requested
polarization. Return a float array of the averaged channels.
EXAMPLE:
assume m[2,31,21] is a map with 2 polarizations, 31 samples per strip,
21 strips, and 1024 frequency channels. Then the call :
avgchn=cmavgchn(m,1,500,520)
returns the array avgchn[31,21] that is averaged over the frequency
channels 500=520.
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cmavgchn.pro)
NAME:
cmavgstrips- average bandpasses for a map
SYNTAX: avgbp=cmavgstrips(m,pol,first=first,last=last)
ARGS:
m[2,nsmp,nstrips] : map structure array
pol : 1 or 2.. pol to return
KEYWORDS:
first : int first strip to use. default =1
last : int last strip to use. default =last
normal : if set then normalize averaged bandpass to unity
RETURNS:
avgbp[nfrqchn] : float .. averaged bandpass.
DESCRIPTION:
Average all spectra in 1 or more strips together. Return the
averaged bandpass.
EXAMPLES:
Assume that m[2,31,21] is a map with 31 samples/strip, 21 strips,and
1024 frequency channels.
avgbp=cmavgstrips(m,1) will average all of the bandpasses for
pol A in the map
avgbp=cmavgstrips(m,1,first=2,last=5) will average all samples in
strips 2 through 5 (counting from 1).
avgbp=cmavgstrips(m,1,first=3,last=3) will average all the samples
in strip 3.
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cmavgstrips.pro)
NAME:
cmgaincor - gain correct map(s)
SYNTAX: cmgaincor,m,gaincor,jy=jy,undo=undo
ARGS:
m[2,nsmp,nstrips,(nmaps)]: array of map structure. routine will work with
a single map or an array of maps
RETURNS:
gaincor[nsmp,nstrips,nmaps]: float. The gain correction done at each point.
m.d: the data that is gain corrected.
KEYWORDS:
jy : if set then gain correct and return in Janskies. The default
is to do a relative gain correction leaving the data in the
input units (normally kelvins)
undo : if set then routine will undo a previous gain correction. In this
case the user passes in m and gaincor.
note: undo is not currently implemented..
DESCRIPTION:
The gain of the telescope in kelvins/Jy is modeled in the routine
corhgainget. This routine will remove the telescope gain variation with
azimuth and za.
By default it will do a relative correction normalizing to the average
gain begin between 5 and 10 degrees za. If the /jy keyword is set then it
will convert the input data (which should already be in Kelvins) to
janskies after correcting for the gain response of the telescope.
The routine performs the following operation on each spectra:
; get the gain value for each az,za position of the map
stat=corhgainget(m.h,gainval) ;here gainval is in K/Jy
if jy keyword is set then
gainCor=gainval ; this is Kelvain/Jy
else
gainCor=(gainval/avggainval)
endelse
m.d=m.d/gainCor
NOTE:
If you input multiple maps, they should be for the same data set
parameters (frequency and receiver).
WARNING:
cmgaincor modifies the input map. You should not call cmgaincor using
already corrected data. eg..
numbpedge=2 ; use 2 spectra each edge for bandpass
mbpc=cormapbc(m,numbpedge)
corgaincor,m,gaincor
corgaincor,m,gaincor ; this is illegal to call a second time
; but..
mbpc=cormapbc(m,numbpedge)
corgaincor,m,gaincor
mbpc=cormapbc(m,numbpedge)
corgaincor,m,gaincor ; this is Ok
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cmgaincor.pro)
NAME:
cmstostr - combine map datasets into an array.
SYNTAX: cmstostr,m,ind,marr,nmaps=nmaps,chn1=chn1,chn2=chn2
ARGS:
m[2,nsmp,nstrips]: {cormap} map data for 1 map
ind: long index for this map
marr[2,nsmp,nstrips,nmaps]:{corget} store maps here
KEYWORDS:
chn1: long . first freq chn to keep. count from 0.def:0
chn2: long . last freq chn to keep. count from 0.def:last
nmaps: long . if specified then allocate marr to hold
this many maps.
DESCRIPTION:
The data returned by cormapinp is an anonymous structure. This allows
it to be dynamically generated. A drawback is that you can't combine
different anonymous structures into an array. This routine allows
you to place cormapinp array structures from calls into an array. Each
map must be the same dimensions. The chn1,chn2 keywords allow you to
keeep only a subset of the frequency channels (to keep the memory usage
down).
EXAMPLE:
alloc=nmaps
for i=0,nmaps-1 do begin
istat=cormapinp(lun,scan,1,2,m,cals) .. get the first record
cmstostr,m,i,marr,nmaps=alloc,chn1=951,chn2=1100
alloc=0
endfor
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cmstostr.pro)
NAME:
cmstrips - overplot spectra by strip.
SYNTAX: cmstrips,m,pol,offset,step,smo=n,title=title,xtitle=xtitle,
ytitle=ytitle,first=first,last=last,freq=freq,vel=vel,
tp=tp,color=color
ARGS:
m[2,xdim,ydim] : {} map array
pol : int 1 or 2 for pol A or pol B
offset : float. the y offset to add to the first strip.
step : float. the y offset to add between strips.
KEYWORDS
smooth : int .. smooth adjacent frequency channels by n.
title : string.. title for plot
xtitle : string.. x axes label
ytitle : string.. y axes label
first : int .. first strip to plot. default is 1
last : int .. last strip to plot. default is end of map
freq : if set then plot versus frequency
vel : if set then plot versus velocity
tp : if set then plot total power instead of spectra.
color : if set then change color for each strip
(you should probably enter ldcolph before running
this so that the color table are loaded..)
DESCRIPTION:
Overplot the spectra for all the strips in a map. Each strip is
offset from the previous by offset.
If the keyword tp is set, then plot the total power by strip for all
strips in the map. Each horizontal line will be the total power for
each sample of a strip. The freq,vel keywords are not used.
You will probably need to use the ver,minver,maxver function to
for the vertical scale to be within the range of the data plotted.
EXAMPLE:
Assume the map is dimensioned m[2,31,21] with 1024 frequency channels.
cmstrips,m,1,.01,.02 will plot
m[0,*,i].d+offset+(i-1)*step for strips 1 thru 21.
cmstrips,m,1,.01,.02,/freq .. will plot versus frequency
cmstrips,m,1,.01,.02,/vel .. will plot versus velocity
To plot total power by strip with zero offset between strips:
ver,.95,1.2
cmstrips,m,1.,0,0,/tp,/color
SEE ALSO:
ver,hor (generic idl routines).
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cmstrips.pro)
NAME:
cmtotpwr - compute total power for a map
SYNTAX: mtp=cmtotpwr(m,chn1,ch2,mask=mask)
ARGS:
m[2,nsmp,nstrips] : map freq structure array
chn1 : first freq chn for average. count from 1
chn2 : last freq chn for average. count from 1
ARGS:
mask[m]: {cormask} mask to use when computing total power.
If supplied, then ignore chn1,chn2. The dimension
m can be:
1: use same mask for entire map
nstrips: use same different mask for each strip
nsmp*nstrips: use different mask for each sample
RETURNS:
mtpavg[2,nsmp,nstrips]: {cormaptp} return array of total power
DESCRIPTION:
Average frequency channels chn1 through chn2 of the map and return
in the mtp structure. It is identical to the m[] structure except that
the m.d[] frequency data is replaced by a single number which is the
averaged total power.
EXAMPLE:
assume m[2,31,21] is a map with 2 polarizations, 31 samples per strip,
21 strips, and 1024 frequency channels. Then the call :
mtp=cmtotpwr(m,100,900)
returns the struct mtp with mtp.d the total power averaged over channels
100:900.
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cmtotpwr.pro)
NAME:
cormapbc - band pass correct a map 1 strip at a time
SYNTAX: mbpc=cormapbc(m,numbpedg,pol=pol,edge=edge,$
m3sections=m3sections,m3type=m3type,$
tsys=tsys,degfit=degfit,cumfilter=cumfilter,$
mask=mask,fractmask=fractmask,coef0=coef0,coef1=coef1
ARGS:
m[2,pntstrip,numstrips] : {} input map array.
numBpEdg : int. This determins the type of bandpass correction to use:
ge 0. numBpEdg is the number of bandpasses on each edge of
a strip to average together and use for the bandpass
correction (0 is the same as 1).
The edge keyword can choose just the left or right edge.
-1 The bandpass correction is computed for each strip by
doing a linear fit for each channel over the samples of
a strip. Outliers are thrown out. The keyword edge will
be ignored.
-2 Use the entire strip to compute the bandpass by averaging
over the strip.
-3 Use the entire strip to compute the bandpass by computing
the median by freq chn along the strip. The m3sections
(Minus 3) and m3type can be used to modify how the
median is computed.
KEYWORDS:
pol : int 1 or 2. if provided, then only process this pol
edge[] : int If numBpEdg is ge 0 then edge can be used to determine
which side of the strips the bandpass correction comes from
(by default numBpEdg bandpasses from both edges are used).
The coding for edge is:
-1 use numbpedg on the starting edge of each strip
for the bandpass correction.
1 use numbpedg on the ending edge of each strip
for the bandpass correction.
0 use numbpedge from both edges.
If edge is a single number (not an array) then the specified
side of the strip will be used for all strips of the map.
If edge is an array then it must be dimensioned to be the
same length as the number of strips in the map. Each
element of edge[] then determines the side to use for
that strip. The number of bandpasses used is still
determined by numbpedg.
m3sections: int When numbpedg is set to -3, the bandpass correction is
computed using the median by channel. Normally this
is done over the entire strip. The m3sections keyword
allows you to break the strip up into m3sections and then
compute the median over each of these sections separately.
By default the minimum (for each channel) of each of these
sections is then used for the bandpass correction.
m3type lets you change how the bandpass value is selected
from the m3 sections.
m3type : int The keyword m3sections creates m3sections bandpasses. By
default the minimum (by channel) of the m3sections bandpasses
is returned as the bandpass correction. m3type lets you
change this to:
0 - The bandpass is the minimum of the m3section bandpasses.
1 - The bandpass is the maximum of the m3section bandpasses.
2 - The bandpass is the average of the m3section bandpasses.
3 - The bandpass is the median of the m3section bandpasses.
tsys: int This keyword determines how to remove the system temperature
from each spectra. The values are:
0 -(default) compute tsys from bpc spectra.
linearly interpolate between each edge if two edges.
if a single edge, use that for the entire strip
1 - linearly interpolate the total power using all of the
spectra of the strip. If numbpedg = -1 then
the average of the linear fits by channel is used.
2 - compute the mean Tsys for each spectra of the strip and
remove it.
3 - do not remove tsys
4 - compute the median Tsys for each spectra of the strip
and remove it.
5 - compute Tsys with a robust average over each spectra and
then remove it.
The mask= or edgefract= keyword can limit the number of
chnannels used to compute tsys.
degfit : int if tsys eq 1 then this is the deg of the polynomial fit to
use. The default is 1.
cumfilter : If tsys eq 1 (and bpnumedg != -1) and this keyword is
set, then a cumfilter is run on the Tsys values before
doing the linear fit. This will normally exlude
the continuum sources from the fit. see cumfilter (under
general software)
mask[nchn]: int mask to use for removing Tsys. If provided then it
should have the same number of channels as m[0,0].d
It will use the non-zero channels to compute and remove
Tsys (depending on keyword tsys). The routine will
always use all of the channels when doing the bandpass
correction (since the cal was averaged over all
channels). Note.. to use this polA and polB must
have the same number of channels.
fractmask : float fraction of channels on each side of spectra to ignore
when computing tsys for tsys removal. Use this instead
of mask= keyword. If you have 1024 channels then
fractmask=.05 would ignore 51 channels from each side.
RETURNS:
mbpc[2,pntstrip,numstrips] : map after bandpass correction
coef0[nlags,nstrips]: float bandpass fit by channel for each strip.
This is the constant.
coef1[nlags,nstrips]: float bandpass fit by channel for each strip.
This is the linear term (-nsmp/2,nsmp/2)
DESCRIPTION:
The routine will bandpass correct a map 1 strip at a time and then
subtract off the system temperature.
BANDPASS CORRECTION:
The bandpass correction is determined by the numbpedge parameter. Its
values are:
numbpedge ge 0.
use a set of bandpasses from 1 or both edges of each strip for the
bandpass correction. numbpedge specifies the number of bandpasses to
average. By default numbpedge spectra from each edge are used. If the
edge keyword is -1 or +1 then only the left or right edge is used for
the entire map.
If edge[nstrips] is an array then you can choose a different edge
for each strip.
numbpedge=-1
The bandpass correction is computed by doing a linear fit to each channel
(over the samples in a strip). The fitting iterates throwing out any
outliers over 2 sigma. The x value for the fit is symmetric about 0 so
the constant term of the fit is used for the bandpass corretion value
at each channel.
numbpedge = -2
The average bandpass over each strip is used for the bandpass correction.
numbpedge = -3
The median bandpass over each strip is used for the bandpass correction.
The m3section keyword lets you break the strip up into m3sections. The
median for each section is computed and then the minimum (by default) of
the m3sections bandpasses is used for the bandpass correction.
The m3type keyword lets you change how the value from the m3sections
bandpasses is selected for each channel. It can be max,min,average, or
median.
SYSTEM TEMPERATURE REMOVAL:
The removal of the system temperature is determined by the TSYS
keyword. The keywords MASK= and FRACTMASK= can be used to limit
the channels that are used to compute tsys (in case there is
interference). The values for TSYS are:
TSYS=0
For each strip compute Tsys from the spectra used for the bandpass
correction. If bandpasses from both edges are used then linearly
interpolate the tsys from the 2 sides across the entire strip.
TSYS=1
Compute tsys for each spectra by taking the average of each spectra. Then
do a linear fit of these values along the strip (the degfit keyword
lets you change the order of the fit). From each spectra remove the fitted
value.
If numbpedge was set to -1 then the linear fits by channel have already
been done. In this case, the average (across the channels) of the fits is
used.
TSYS=2
Compute Tsys for each spectra and remove it from each spectra. This is
handy if you want to get rid of the continuum sources.
TSYS=4
Same as Tsys=2 but use the median rather than the mean.
TSYS=5
Same as Tsys=2 but use a robust average. All points 3 sigma away from
the mean are excluded, and then the mean is recomputed. Continue looping
till no points are thrown out.
TSYS=3
Do not remove Tsys.
Take a look at: Comparing mean,median,robustAvg for a comparison of tsys=2,4,5.
For each strip the routine will:
1. Average the bandpasses for the bandpass correction.
2. For the i'th spectra in each strip compute:
spectra[i]=(spectra[i]/normalized(avgBP) -Tsys[i]
where normalized(avbBp) is the normalized average bandpass computed
in step 1. Tsys[i] is the system temperature computed from the
keywords tsys= etc..If mask= or fracmask= is set then avgBP and
Tsys are computed over the specified channels.
EXAMPLES:
1. Bandpass correct using 3 channels on each side of the map. Remove tsys
using all spectral channels interpolating between the 3 averaged bandpasses
on each edge.
m=cormapbc(m,3)
2. Bandpass correct using only the starting edge. When computing Tsys
ignore 5% of the channels on each edge.
m=cormapbc(m,3,edge=-1,fractmask=.05)
3. Suppose you want to make a continuum map..
Bandpass correct linearly interpolating along each channel. Remove
tsys by using the average (along channels) of these channel fits.
Do not use 5% of the channels on each edge of the spectra.
m=cormapbc(m,-1,fractmask=.05,tsys=1)
4. You want to remove the continuum sources (ie. you are looking for galaxies).
Use the linear fit by channel to compute the bandpass correction.
Remove Tsys by substracting the median of each spectrum. Ignore 5%
of the channels on each edge of the spectra when computing tsys.
m=cormapbc(m,-1,fractmask=.05,tsys=4)
5. Same as 4, but you get a negative image of some galaxies in your strip.
This is probably because continuum sources in the bandpass correction
are cutting down on the number of unbiased samples for the median. Try
breaking the strip up into 3 sections, computing the median of each
Then try removing Tsys using a robust average over each spectra
(ignoring 5% of the channels on each edge).
m=cormapbc(m,-3,fractmask=.05,tsys=5,m3sections=3)
RECENT HISTORY:
21nov04 - added m3sections,m3type keywords.
07jul04 - added tsys=5. use avgrobbychan for tsys removal.
06jul04 - added tsys=4. uses median rather than mean for Tsys removal.
07feb04 - bandpass by channel, after fit, recheck all the
residuals for being with nsig, not just the current good ones.
Bug in the residual compare, was not using abs().
14jan04 - allstrip -> numbpedg=-2
07jan04 - added bcbychan option. does a linear fit by chan along
each strip to determine the bandpass.
03sep03 - if mask supplied (or edge fract) then the normalizatino
of the bandpass for bandpass correction should be
over the mask, and not the entire bandpass. Before it was
normalizing over the entire bandpass, dividing, and then
using the mask to compute the tsys subtraction.
27jan03 - added cumfilter keyword
27jan03 - for tsys removal if tsys=1 then allow option to
iterate throwing out pnts above 3sigma
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cormapbc.pro)
NAME:
cormapbl - baseline all the spectra in a map.
SYNTAX: mbl=cormapbl(m,mask,polyDeg,gotmask=gotmask)
ARGS:
m[2,nsmps,nstrips] - map array of structures.
mask[nfrqchn]: float - mask to use for baselining
polyDeg : int - degree of polynomial for baselining
KEYWORDS:
gotmask : if set then user supplies mask and polyDeg as input.
RETURNS:
mask[nfrqchn] : the mask array used.
polyDeg : the degree of the polynomial used.
mbl[2,nsmps,nstrips] : baselined map info
DESCRIPTION:
cormapbl will baseline all the spectra in a map. The user inputs
the map array. The routine call bluser which allows the user to
interactively define a mask and order for baselining. After the
mask and order is chosen, the routine will fit the requested order
polynomial to every spectra using the mask.
A copy of the map array is returned with the baseline removed from the
data.
bluser is called from within the routine with the average of all
spectra in pola. The users uses this as the reference to decide which
polynomial and mask to use.
EXAMPLE:
let m[2,31,21] be the map structure with 1024 frequency channels.
the call:
mbl=cormapbl(m,mask,polyDeg) will:
1. call bluser with average of polA over the map.
the output from bluser consists of:
KEY ARGS FUNCTION
m .. define mask
f n .. fit polynomial of order n
h h1 h2 .. change horizontal scale for plot to h1,h2
v v1 v2 .. change vertical scale for plot to v1,v2
c .. print coefficients
p .. plot data - fit
q .. quit
You must first specify the mask with m.After than you can
fit various polynomials with the f deg . Use p to plot the
data - fit. Use v 1.1 1.5 or h 200 400 to change the vertical and
horizontal display.
When you are satisfied with the fit, enter q and the routine will
then fit the polynomial to each spectra and return it in mbl.
NOTE:
If you called: cormapinp, cormapsclk, and then this routine, the
spectra still need to be bandpass corrected before the temperature
scale is correct.
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cormapbl.pro)
NAME:
cormapinp - input a correlator map.
SYNTAX:istat=cormapinp(lun,scan,polABrdNum,polBBrdNum, m,cals,norev=norev,
han=han,sl=sl,maxrecs=maxrecs,avgsmp=avgsmp,
calrecs=calrecs)
ARGS:
lun: int logical unit number for file. It should already be opened.
scan: long scan number for start of map.
polABrdNum: int correlator board index (1 thru 4) to take polA data from.
polBBrdNum: int correlator board index (1 thru 4) to take polB data from.
KEYWORDS:
norev: When you drive both directions, the routine will normally
reverse the odd strips so that the data,header, and cals of
adjacent strips will line up in RA. Setting norev to true
will override this default. The data will be returned in
the order it is taken.
han: if set then hanning smooth the data on input.
sl[]: {sl} scanlist array returned from call a previous call to
getsl(). If this keyword is provided then direct access is
used rather than sequentially reading the file.
maxrecs: long If you have more than 300 records in a scan you will have
to set maxrecs to that value.
avgsmp: long Number of samples to average together. Handy if you've
oversampled each strip. If this number does not divide
evenly into the samples/strip, then samples at the end
will be dropped.
avsmp=0 or 1 is the same as no averaging
RETURNS:
istat: int 1: got all the strips
0: got part of requested map
-1: got no data.
m[2,pnts/strip,nstrips]:{} array of structures holding the returned
data and header. (see below for a description).
cals[nstrips*n] :{} array of structures containing the processed
cal on,off data (cal value and scale factor)
(see below for a description).
calrecs[nstrips*n]:{calrecs} array holding the raw cal on, cal off
input data. For the i'th calOn/off measured the
the structure contains cal spectra as float arrays:
calrecs[i].scan ..scan number forcalOnScan
calrecs[i].calPerStrip ..0,1,or 2
calrecs[i].con[nchans,2]
calrecs[i].coff[nchans,2]
where i is the i'th calon/off, nchannels are the
number of frequency channels, and [,2] has polA, polB
You can use this data to do rfi excision and then
recompute the cals[] array.
DESCRIPTION:
Input a map taken with the correlator routines cormap1,cormapdec,
cordrift, or cormapbm. You open the datafile (eg openr,lun,filename,/get_lun)
prior to calling this routine and then pass in the lun. The scan argument
tells the routine where to start reading. This should be the first scan of the
map. If the map starts with a calon/off and you gave the scan of the
first strip, the routine will search backwards looking for the cal that goes
with the first strip (so either is ok to enter). You need to also
specify which boards of the correlator contain the polA, polB data you
want to extract (it will only process 1 set of polA, polB data at a time).
The board numbers are indexed 1 through 4. If polA, polB were input on
on the same board, then use the same board number for both.
The routine will figure out from the header the size of the map, whether
it is a partial map, which directions you drove, and how the cals were taken.
The map information is returned as an array m of structures. The
dimensions of the array are:
m[pol,pntsperstrip,numstrips]. pol=2
Each element of the array contains the information for a particular
sample and polarization. The contents of m (for this example assume
2048 freq channels,and we are looking at sample 5 of strips 7
(counting from 1):
m[0,5,7].h the header for pola,
m[0,5,7].d[2048] the data for pola
m[0,5,7].p the total power value for this sample (linear scale).
m[0,5,7].az the azimuth position in degrees at the end of
each sample.
m[0,5,7].za the zenith angle position in degrees at the end of
each sample.
m[0,5,7].azErrAsec the azimuth tracking error in arcseconds at the
end of each sample (great circle).
m[0,5,7].zaErrAsec the zenith angle tracking error in arcseconds at the
end of each sample.
m[0,5,7].raHr the RA in hours at the middle of each sample.
m[0,5,7].decDeg the declination in degrees at the middle of each sample.
m[0,5,7].calscl cal scale factor for this position. calst,calend, or
(calst+calEnd)/2.
m.h holds the header information for this sample. It contains
several sub structures:
m.h.std - the standard header
m.h.pnt - the pointing header
m.h.iflo - the iflo header
m.h.dop - the doppler,freq,velocity header
m.h.cor - the correlator header
m.h.proc - the datataking procedure part of the header.
Some locations of interest are:
m.h.std.scannumber: the scan number for the strip
m.h.pnt.r.reqPosRd[2]: requested ra,dec for center of map in radians.
If you drove both directions in ra then lining up the samples numbers from
each strip would not give a map in increasing ra. By default the routine
will flip the order of the even strips (counting from 1) in the entire
map (it flips the order of the headers and cals as well as the data). eg..
suppose the map has 15 strips, counting from 1 it will flip the samples
in strips 2,4,6,8,10,12,14,12,14 so the ra will line up.
If you started in the middle of a map, then it will flip the strips relative
to the number in the entire map. eg.
suppose the map has 15 strips and your first strip was strip 4, then it
would flip 1,3,5... of the returned strips (since these would normally
be flipped in the entire map).
Setting the /norev keyword will override the default and return the data
in the order it was taken in each strip.
The cal data is returned in an array of cal structures. It is
dimensioned cals[nstrips*n]. n is 1 if you took 1 cal per strip and 2 if
you took 2. The contents of the cal structure are:
cals[i].h - the header for the on cal record.
cals[i].calval[2] - the cal value in kelvins for polA,polB interpolated
to your frequency.
cals[i].calscl[2] - the cal scaling factor to go from correlator counts
to kelvins for pol A,polB. it is:
computed as : calValKelvins/(calOn-caloff)
The correlator part of the header contains the cal on, cal off
counts:
cals[i].h.cor.calon[2] polA,polB cal on total power
cals[i].h.cor.caloff[2] polA,polB cal off total power
If you provide the calrecs keyword, then the actual cal spectra will
also be returned in an array of calrecs structure. Each element contains
the info for one cal on/off measurement. The contents are:
calrecs[i].scan ..scan number forcalOnScan
calrecs[i].calPerStrip ..0,1,or 2
calrecs[i].con[nchans,2].. cal on spectra correlator units
calrecs[i].coff[nchans,2] cal off spectra correlator units
nchans are the number of frequency channels per spectra. If there
are two sbc per board then polA will be [nchans,0], polB will be [nchans,1]
If there can be 0,1 or 2 cals per strip.
This data can be used to recompute the cal values after rfi excision.
EXAMPLE:
Suppose we are mapping with 15 samples/strip, 11 strips, and two correlator
boards with 1024 pola/polB on each board. We also took 1 cal per strip.
Since arrays are indexed from 0,lets count our samples and strips from 0.
Then:
m[0,4,3].d[*] polA sample 4,strip 3
m[1,4,3].d[*] polB sample 4,strip 3
the headers would be:
m[0,4,3].h polA header sample 4,strip 3
m[1,4,3].h polB header sample 4,strip 3
Assuming there is 1 cal rec/strip, you could scale the spectra in
the first strip to kelvins with:
m[0,*,0].d=m[0,*,0]* cals[0].calscl[0]
m[1,*,0].d=m[1,*,0]* cals[0].calscl[1]
If we took two cals per strip then the cals for the i'th strip (counting
from zero would be:
cals[2*i:2*i+1]
SEE ALSO:
cor2/corget, cor2/cormapsclK, cor2/arch_getmap,cormap/cormapinplist
NOTES:
1. For lbw the routine assumes the cal values are for the measured
polarizations. If the circular hybrid is in, then the cal values from
the file (which are linear, need to be averaged). This probably needs a
parameter to flag this (or maybe rely on the header to have the iflo info..)
2. cormapbm does a beam map in great circle azimuth,za relative to the
center of the source. The m.raHr and m.decDeg values are the az,za offsets
(great circle) in degrees rather than the ra and dec (sorry about this
but it was too late to change the variable names to m.c1,m.c2).
m.raHr -- azimuth offset relative to src center in degrees.
m.decDeg -- zenth angle offset relative to src center in degrees.
3.The header contains a proc section that holds procedure dependant information
m.h.proc.iar[10] integer data
m.h.proc.dar[10] double data
These contain some of the parameters that were used in the datataking
procedure:
---------------------
HDR.PROC.IAR definition for cormap1,cormapdec
code : bit 0: coord : 0 J2000, 1 1950
: bit 1: rawidth : 0 great cirle, 1 little circle ra width
: bit 2: direction: 0 both ways , 1 increasing ra only
: bit 3: update doppler odd strips
: bit 4: update doppler even strips
: bit 5: cal on/off start of strip
: bit 6: cal on/off end of strip
: bit 8: cal on/off start of map (once only)
: bit 9: 1--> adjust correlator power levels at start position.
iar[0] - code value used
iar[1] - number of strips requested in map (total.. not this call)
iar[2] - 1st strip to start on (count from 1)
iar[3] - integrations per strip
iar[4] - stripNumber this strip (count from 1)
(from start of map, not start of call)
iar[5] - start time seconds from midnite ast this strip
HDR.PROC.dAR definition for cormap1,cormapdec
dar[0] - seconds per integration
dar[1] - raOffAmin for map .. what user entered
dar[2] - decOffAmin for map .. what user entered
dar[3] - decStepAmin for map .. what was used
strip info
dar[4] - raRateRd for map ..
dar[5] - raRateRd this strip little circle
dar[6] - raOffRd this strip little circle
dar[7] - decOffRd this strip
---------------------
HDR.PROC.IAR for cordrift
code code where each bit signifies an option:
code code where each bit signifies an option:
1 bit0: coord : 0 J2000, 1 1950
2 bit1: cal on/off end of strip: 1 yes, 0 no
4 bit2: 0:raOffset=stripTime/2.,
1:raOffset=stripTime*1.0027../2 .. convert to sidereal secs
8 bit 3: update doppler firstStrip (only)
0x10 bit 4: update doppler each strip (including first)
0x100 bit 8: coradjpwr at start of 1st strip
0x200 bit 9: coradjpwr at start of every strip
0x400 bit10: cal onoff start of each strip
iar[0] - code value used
iar[1] - requested number of times to loop (strips if decStp != 0)
iar[2] - requested seconds to drift
iar[3] - integrations per strip
iar[4] - current strip count from 1
iar[5] - start time seconds from midnite ast this strip
iar[6] - extra settle time secs..
iar[7] - stripsPerStep in dec
HDR.PROC.dAR definition for cordrift
dar[0] - seconds per integration
dar[1] - haOffRd for map .. what user entered current coord..
dar[2] - decOffAmin for map .. what user entered
dar[3] - decStepAmin for map .. what user entered
strip info
dar[4] - decOffsetThis strip amin
---------------------
HDR.PROC.IAR for cormapbm (beam maps)
code code where each bit signifies an option:
: bit 0: coord : 0 J2000, 1 1950
: bit 1: azwidth : always great circle
: bit 2: direction: 0 both ways , 1 increasing az only
: bit 3:
: bit 4:
: bit 5: cal on/off start of strip
: bit 6: cal on/off end of strip
: bit 8: cal on/off start of map (once only)
: bit 9: 1--> adjust correlator power levels at start position.
iar[0] - code value used
iar[1] - number of strips requested in map (total.. not this call)
iar[2] - 1st strip to start on (count from 1)
iar[3] - integrations per strip
iar[4] - stripNumber this strip (count from 1)
(from start of map, not start of call)
iar[5] - start time seconds from midnite ast this strip
HDR.PROC.dAR definition for cormapbm
dar[0] - seconds per integration
dar[1] - azOffAmin for map .. great circle
dar[2] - zaOffAmin for map ..
dar[3] - zaStepAmin for map .. w
strip info
dar[4] - azRateRd for map ..(great circle
dar[5] - azRateRd this strip great circle
dar[6] - azOffRd this strip great circle
dar[7] - zaOffRd this strip
modhistory
31jun00 update for new corget format
20dec00 updated to use h.proc information
09jun01 added support for sl keyword
04nov01 when computing positions, the rate seconds are sidereal. need
to compute duration solar->sidereal
07dec01 fixed up m.rahr for cordrift. Had constant position for each strip.
16dec02 added maxrecs
16dec02 added avgsmp. updated:
all
h.cor.dumpsPerInteg
h.proc.iar[3] - integrations per strip
h.proc.dar[0] secs per integration
cordrift
iar[2] secs req to drift
07apr03 check procname and source name..
01jun03 idlversion >= 5.5 now treat embedded structures
as arrays of 1. This breaks the reform code. Need to test
version of idl to do reform correctly.
05dec03 start adding support for cormapbm
10jan04 it now searches backwards for the first cal if the user
enters the scan number of the first strip..
03feb04 was not computing number of cals to return if cals on both
edge of strip
20feb04 newer versions of idl were creating different dimensions of the
arrays .. eg d[128,1,30] and d[128,30]. This was causing problems
on the reversing.. I tried to fix it here..
07feb04 case statement line 580.. first test was cormap1 or cordrift
2nd test was cordrift, should have been cormapdec. This caused
cormapdec patterns to fail.
26jul04 need to check rcvr number. Some clever people were interleaving
maps of same name, but different receivers..
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cormapinp.pro)
NAME:
cormapinpcal - process a cal on/off for cormapinp
SYNTAX: stat=cormapinpcal(lun,b,calAr,brdIndA,brdIndB,calrec,totcals,stripNum,
rdrec=rdrec,han=han)
ARGS:
lun : int unit number assigned to file to read.
b : {corget} optional. calOn rec if already read in.
brdIndA : int board index to process for polA count(0..3)
brdIndb : int board index to process for polB count(0..3)
calrec : int index into calAr[] to place this cal data
If sucessful, this variable gets incremented on return.
totcals : int total number of cals that will be placed in calAr.
On the first cal with calrec=0, totcals is used
to initially allocate calAr[totcals].
stripNum: int the strip number for this cal (count from 1)
KEYWORDS: rdrec if set, then the cal on rec should be read from
the current location of lun and placed in b. If rdrec
is not set, then the user has already read the
the calon rec and it is passed in via b.
han: if set and returning calrecs, hanning smooth the
cal spectra.
RETURNS:
stat : int 1 successful, 0 error.
calAr[n]: {calstr} array holding calinfo for each strip of map.
The data for the current cal will be stored in
calar[calrec] (before calrec is incremented).
calrec : int will be incremented on successful return. This keeps
track of where the next cal will go in calAr.
b :{corget} If rdrec is not set, then the calOn rec is returned
DESCRIPTION:
This routine is normally called by cormapinp. It will process the
next calonoff pair of scans. If rdrec is set,
then lun should point to the calOnRec. If rdrec is not set, then the
cal on rec should already have been read into the variable b and lun
should point at the cal off rec.
The routine calls corcalonoffrec(). It loads the cal results in
the strcuture calAr[calrec]. Calrec is then incremented.
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cormapinpcal.pro)
NAME:
cormapinplist - input a map from a list of files.
SYNTAX:istat=cormapinplist(flist,scanlist,polABrdNum,polBBrdNum, m,cals,
maxstrips=maxstrips,norev=norev,han=han,
maxrecs=maxrecs,avgsmp=avgsmp)
ARGS:
cormapinplist arguments:
flist[n]: string array of files to read data from
scanlist[n]: long starting scannumber of the map in each file
cormapinp arguments:
polABrdNum: int correlator board index (1 thru 4) to take polA data from.
polBBrdNum: int correlator board index (1 thru 4) to take polB data from.
KEYWORDS:
maxstrips: long The maximum number of strips for the completed map
The program will allocate an array of this size when
inputting the data. The default max size is 41 strips
The following keywords are from cormapinp.
norev: When you drive both directions, the routine will normally
reverse the odd strips so that the data,header, and cals of
adjacent strips will line up in RA. Setting norev to true
will override this default. The data will be returned in
the order it is taken.
han: if set then hanning smooth the data on input.
sl[]: {sl} scanlist array returned from call a previous call to
getsl(). If this keyword is provided then direct access is
used rather than sequentially reading the file.
maxrecs: long If you have more than 300 records in a scan you will have
to set maxrecs to that value.
avgsmp: long Number of samples to average together. Handy if you've
oversampled each strip. If this number does not divide
evenly into the samples/strip, then samples at the end
will be dropped.
avsmp=0 or 1 is the same as no averaging
RETURNS:
istat: int 1: got map
-1: trouble getting the data
m[2,pnts/strip,nstrips]:{} array of structures holding the returned
data and header. (see below for a description).
cals[nstrips*n] :{} array of structures containing the cal on,off data
(see below for a description).
DESCRIPTION:
cormapinp is used to input an on the fly map from a data file. Large
maps may be taken a few strips per day and be spread over many datafiles.
cormapinplist will input a map that is spread over many files.
See cormapinp for a description of the m,cals structures that are
returned.
EXAMPLE:
Suppose the map is spread over 3 days and that the files are
located in /share/olcor/ with the specified starting scan numbers. The
data is in brdpolA=1 and brdPolB=2 with a total of 49 strips.
You could input the map using:
flist=['corfile.14dec02.a1632.1','corfile.15dec02.A1632.1',$
'corfile.16dec02.A1632.1',$
dir='/share/olcor/'
scanlist=[234800055L,234900056L,235000053L]
polABrd=1
polBBrd=2
istat=cormapinplist(dir+flist,scanlist,polABrd,polBBrd,m,cals,maxstrips=49)
SEE ALSO:
cormapinp(), cor2/arch_getmap
modhistory
19dec02 started
07mar04 added calrecs keyword
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cormapinplist.pro)
NAME:
cormapsclK - scale a map to Kelvins.
SYNTAX: istat=cormapsclK(m,cals,tsys,use1cal=use1cal)
ARGS:
m[pol,smpPerStrip,nstrips] : map array input from cormapinp.
cals[ncals] : cals data input with cormapinp.
KEYWORDS:
use1cal: int 1 or 2. if two cals per strips you can should to use
only the first or 2nd cal by setting use1cal to 1 or 2.
RETURNS:
m.d : will be scaled to kelvins
tsys[2,smpPerStrip,nstrips] : average system temp each sample
istat : 1 ok, 0.. trouble
DESCRIPTION:
Scale the map data to kelvins. To use this routine the data must
have been taken with 1 or 2 cals per strip. For each strip in the map
the cals.calscl scale factor is used to convert from correlator counts
to kelvins. If two cals per strip were taken then the cal values are
interpolated across the strip.
The original data in m.d is overwritten with the new values.
The location m.calscl is also updated if we have 2 cals/strip. The
interpolated value is used instead of the average that was loaded in
cormapinp.
EXAMPLE:
istat=cormapinp(lun,scan,brda,brdB,m,cals)
istat=cormapsclk(m,cals,tsys)
When calling this routine, make sure that the data array m is still
in correlator units (as returned by cormapinp). If you call cormapbc()
to bandpass correct the data, then the returned data array is no
longer in correlator units and this routine should not be used with
that data array.
NOTE:
1. Since the original data is overwritten, you do not want to run this
routine more than once on the same array. The conversion only works if
the input data is in correlator counts.
2. This is normally the first routine called. After calling the routine,
you will notice that the system temperatures in tsys are smaller than
the central portion of the spectra. This is because the bandpass correction
has not yet been done. Tsys is the average of the entire spectra. Since
the spectra falls off at the edges (because of the bandpass) the
center or the spectra will be higher than the average. You will have
the correct temperatures only after the bandpass correction is done (
by dividing through by a normalized bandpass).
This implies that scaling to kelvin and then fitting a baseline
will not give the correct temperatures. If this is really what is
wanted, then the scaling factor from correlator units to kelvins
should be computed over the central part of the spectra.
(See /pkg/rsi/local/libao/phil/Cor2/cormap/cormapsclk.pro)