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.. 01jun03idlversion >= 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)