NAME: abssq - compute absolute value then square SYNTAX: pwr=abssq(v) ARGS: v[n]: complex voltage to compute power RETURNS: pwr[n]: float real(v)*real(v) + img(v)*img(v) DESCRIPTION: Compute the square of the absolute value
(See /pkg/rsi/local/libao/phil/gen/abssq.pro)
NAME:
addpath - add a directory to the start of the path variable
SYNTAX: addpath,pathname
ARGS:
pathname : string variable with the path to add. If the pathname
does not begin with / or ~ then the path will be
relative to the path returned by aodefdir()
DESCRIPTION:
Add a directory to the beginning of the !path variable. If the
first character does not begin with ~, or /, then make the path
relative to the directory returned by aodefdir(). If the pathname is already
in the path, then move it to the front of the path.
EXAMPLE:
addpath, '/home/aosun/u4/bozo/idl' ..
addpath, 'Cor2' .. adds /pkg/rsi/local/libao/phil/Cor2
(See /pkg/rsi/local/libao/phil/gen/addpath.pro)
NAME:
alfatsysget - return alfa Tsys for the requested positions
SYNTAX: stat=alfatsysget(az,za,freqMhz,rotAngle,tsysValAr,date=date,fitI=fitI,
fname=fname)
ARGS:
az[n]: float azimuth in degrees
za[n]: float zenith angle in degrees
freqMhz[n]: float freq in Mhz.
rotAngle[n]: float rotation angle of alfa in degrees.
KEYWORDS:
date[2]: int [year,daynumber] to use for gain computation. The default
is to use the current day. If the Tsys curves change
with time, this allows you to access the curve that
was in use when the data was taken.
fname: string alternate file holding fit coefficients.
RETURNS:
tsysValAr[2,7,n]: float tsys for the 2pols, 7 pixels and n positions in Kelvins
fitI[2,7] : {alfatsysfitI} optionanly return the fit information (see alfatsysgetdat for
a description of the structure).
stat: int -1 --> error, no data returned
1 --> data returned ok
DESCRIPTION:
Return the system temperature (K) for the 2*7 pols*pixels of alfa using a
model of az,za,rotation Angle, and frequency. The date keyword allows you
to access a tsys curve that was valid at some other epoch (the default is the
most recent curve). If the input variables are an array of dimension N then
the retured data will be an array of Tsys(2,7,N)
Fits have been done for Tsys(az,za,freq,rotationAngle). This routine will
input the fit information (calling alfatsysinpdata) and then compute the Tsys for the
given input parameters.
The input fit data is stored in a common block so the data does not have to
be input from disc on any succeeding calls.
NOTE:
For a description of the Tsys fits see:
http://www.naic.edu/~phil/mbeam/Tsys/Fits/fittingTsys.html .
EXAMPLES:
; 1 value at az,za,rotangle,freq
az=180.
za=15.
rotAngle=19.
freq=1385.
stat=alfatsysget(az,za,freq,rotAngle,tsysVal)
;
; 18 values at za 2 thru 19 in 1deg steps, az=180,freq=1385.
n=18
az=fltarr(n) + 180.
za=findgen(n)+2 ; za=2..19
rotA=fltarr(n) + 19.
freq=fltarr(n) + 1385.
stat=alfatsysget(az,za,freq,rotA,tsysVal)
SEE ALSO:alfatsysinpdata
(See /pkg/rsi/local/libao/phil/gen/alfatsysget.pro)
NAME:
alfatsysinpData - input alfa tsys fits.
SYNTAX:
istat=alfatsysinpdata(tsysfitI,fname=fname,date=date)
ARGS:
KEYWORDS:
fname: to specify an alternate data file with fit values.
The default file is aodefdir() + 'data/tsys.datR17
date : [year,daynum] .. if specified the data when you want the
tsys for.. the default is most recent.
RETURNS:
istat: 1 ok, -1 bad filename or data in file.
tsysfitI[14]:{alfatsysFitI} return fit info. 1 structure for each
pol/pix.
DESCRIPTION:
Input the alfa tsys fit data. The default datafile is aodefdir() +
'data/tsys.datR17 (aodefdir() is a function that returns the root of the ;
aoroutines). The keyword fname allows you to specify an alternate file.
The file format is:
- col 1 ; or # is a column
- col 1 !yyyy dayno starts a date section. yyyy dayno is the
year daynumber for the start of this data set.
- data is free format , column oriented
The structure format for {alfaTsystFitI} is:
fitI.pol 0,1 ; polA or polB
fitI.pix 0..6 ; pixel number
fitI.fitType 1 ; code for type of fit used
fitI.pntsused 1 ; number of points used to compute the fit
fitI.ncoef 11 ; number of coef in the fit
fitI.sigmafit 0. ; hold the fit sigma (in deg K)
fitI.numCoef fltarr(ncoef); hold the fit coefs.
fitI.sigmacoef fltarr(ncoef); hold the fit coef errors.
fitI.startYr for the fit
fitI.startDaynum for the fit
How the different cal routines vary:
alfatsysinpdata() inputs the data from disc.
It defaults to the current date. It loads a table in common
holding the fit info for all of the pixels/pols
alfatsysget() Pass in the az,za,freq,rotAr. It will input the
data using alfatsysinpdata if necessary, do the computation
and return the tsys for the reqeusted values.
SEE ALSO: alfatsysget
(See /pkg/rsi/local/libao/phil/gen/alfatsysinpdata.pro)
NAME:
alfawappharm - compute alfa,wapp harmonics
SYNTAX: n=alfawappharm(rfCfr,freqList,harmList,maxHarm=maxHarm,rfidef=rfidef,$
print=print, bpf=bpf)
ARGS:
rfCfr : float sky center frequency (Mhz) of alfa band
freqlist[m]: float sky frequency of birdies to check (Mhz) (also see
rfidef keyword below).
KEYWORDS:
maxHarm:int largest harmonic to use (default is 5).
rfidef:int if this keyword is set then the known large interfering
frequencies are loaded into freqList and then used for the
computation. See below for a list of the frequencies included.
print : if set then print out the results when done
bpf[2] :float redefine the min,max frequency of the 250 Mhz IF
bandpass filter. The default is 200 to 300. Since the filter
is not infinitely sharp, you can try broadening it a bit.
RETURNS:
n: int number of harmonics found in the final band. This includes
any fundamentals
harmList[n]: struct structure holding the harmonics info.
DESCRIPTION:
Compute where harmonics fall in the alfa/wapp configuration. The
user inputs the sky center frequency of the alfa band as well as
a list of frequencies that may create harmonics. The program computes
where these frequencies and their harmonics end up in the final output band.
By default the program searches up to the 5th harmonic. You can change this
with the maxHarm keyword.
The program sequence is:
0. initialize the list to include all of the fundamental frequencies.
1. Add all harmonics created by the dewar.
2. Remove all frequencies outside the 1225-1525 bandpass filter.
3. Compute harmonics created at the IF in the mixer.
4. Remove all frequencies outside the 700 Mhz low pass filter.
5. Create harmonics in the IF created by the fiber/distribution amp.
6. remove all frequencies outside the 250Mhz bandpass filter
7. Compute harmonics created by the amp after the bandpass filter and the
a/d.
8. Digitally downconvert to 0-100
9. Compute where these frequencies map to on the sky.
The bpc keyword lets you change the default 100 Mhz bandwidth of the
250Mhz bandpass filter. The program assumes the filter is infinitely
steep. You can widen the 100 Mhz to see if there are any nearby birdies
that might alias in.
If the rfidef keyword is set (/rfidef), the program will load a list
of known interferers into the rfiList variable and then proceed with the
computation. The list includes:
aerostat radar: 1241.75,1244.6,1256.5,1261.25
remy radar: 1270,1290
faa radar: 1330,1350.
puntaSalinas: 1232.7, 1247.7 (The 2 modeA freqs that don't overlap
with the aerostat)
The return structure harmList[] contains:
RFCFR FLOAT 1385.0 ; center frequency of alfa band
BW250 FLOAT 100.0 ; bandwidth of 250Mhz bandpass filter
MAXHARM INT 5 ; maximum harmonic used for search
FREQRF FLOAT 1350.0 ; The rf frequency of the initial sky birdie
FREQRFH FLOAT 1350.0 ; The rf frequecy of the resulting harmonic
FREQI FLOAT 285.0 ; the IF frequency of the harmonic
FREQH FLOAT 85.0 ; the harmonic frequency at base band
HARMNUM INT 0 ; the order of the harmonic 0=1= the fundamental
STAGENAMH STRING 'sky'; The stage where the harmonic was created.
STAGENAMI STRING 'sky'; The stage that preceded the harmonic creation
You can print out the harmonic list useing alfawappharmpr,harmList
EXAMPLE:
1. pass in faa radar frequencies.
rfiFreq=[1330.,1350]
rfCfr=1385
n=alfawappharm(rfcfr,rfiFreq,harmL)
alfawappharmpr,harmL
Alfa harmonics. SkyFreq:1385.0 250BPwidth:100.0 MaxHarm: 5
freqRfi freqSrc harmNum LocationCreated
1350.0 1350.0 0 sky
1405.0 1350.0 2 AmpAfter100MhzFilter_AtoD
1380.0 1350.0 3 AmpAfter100MhzFilter_AtoD
1375.0 1350.0 4 AmpAfter100MhzFilter_AtoD
1410.0 1350.0 5 AmpAfter100MhzFilter_AtoD
2. use the default frequencies
n=alfawappharm(1375,freqList,harmL,/rfidef,/print)
Alfa harmonics. SkyFreq:1375.0 250BPwidth:100.0 MaxHarm: 5
freqRfi freqSrc harmNum LocationCreated
1330.0 1330.0 0 sky
1350.0 1350.0 0 sky
1415.0 1330.0 2 AmpAfterIFbpFilter_AtoD
1375.0 1350.0 2 AmpAfterIFbpFilter_AtoD
1340.0 1330.0 3 AmpAfterIFbpFilter_AtoD
1400.0 1350.0 3 AmpAfterIFbpFilter_AtoD
1405.0 1330.0 4 AmpAfterIFbpFilter_AtoD
1325.0 1350.0 4 AmpAfterIFbpFilter_AtoD
1350.0 1330.0 5 AmpAfterIFbpFilter_AtoD
1400.0 1350.0 5 AmpAfterIFbpFilter_AtoD
NOTES:
1. This program maps the sky frequency into the final bandpass. It does
not map the location of a birdie in the final bandpass back to the
sky frequency.
2. The harmonic number tells how fast the birdie will move as you move
the center frequency of the sky band.
(See /pkg/rsi/local/libao/phil/gen/alfawappharm.pro)
NAME:
aobeampatcmp - compute ao beam pattern from bessel function
SYNTAX n=aobeampatcmp(freqMhz,diam=diam,nsdlb=nsdlb,thMin=thMin,thMax=thMax,$
valMax=valMax,fwhmA=fwhmA,np=np,p=p,thD=thD)
ARGS:
freqMhz: float frequency in Mhz to use.
KEYWORDS:
diam: float dish diameter in meters to used. The default is 225 meters.
nsdlb: int number of sidelobes,nulls to compute. The default is 10.
np : long number of points to use in the beam pattern. The default
is 10000 points.
RETURNS:
n : int number of sidelobes computed.
thMin[n]: float angle (arcminutes) where nulls were found
thMax[n]: float angle (arcminutes) where peaks were found (includes
main beam at 0 degrees).
valMax[n]: float value at each peak (mainbeam and sidelobes).
normalized to mainbeam = 1.
fwhmA : float full width at half maximum in arcminutes.
p[np] : float beam pattern linear scale.
thD[np] : float the angle (from center of beam pattern) for each point
in p[n]. Units are Degrees.
DESCRIPTION:
Compute a beam pattern for the arecibo telescope. See
tools of radio astronomy, rohlfs and wilson page 139.
The computation uses a circular aperture with uniform illumination (
not quite the ao taper but...).
You specify the frequency in Mhz. The routine will compute and plot
the beam pattern listing the first nsdlb nulls and sidelobes on the plot.
It can pass back to the caller the computed information.
(See /pkg/rsi/local/libao/phil/gen/aobeampatcmp.pro)
NAME:
aodefdir - AO base directory for idl routines.
SYNTAX: defdir=aodefdir(doc=doc,url=url)
ARGS:
doc: if the keyword is set then return the directory for the
html documentation.
url: if the keyword is set then return the url for the
html documentation.
DESCRIPTION:
Return the directory where the ao idl routines are stored. At AO it
returns '/pkg/rsi/local/libao/phil/'. The addpath() routine will use this
directory if no pathname is given. This routine makes it easier
to export the ao idl procedures to other sites.
(See /pkg/rsi/local/libao/phil/gen/aodefdir.pro)
NAME:
arrpltazzaerr - arrow plot of azerr,zaerr vs az,za
SYNTAX: arrpltazzaerr,az,za,azerr,zaerr,tit,ticklen=tickLen,hard=hard
indlist=indlist,names=names,notab=notab,tabtit=tabtit,
colar=colar,ln=ln
ARGS:
az[npts]: float azimuth deg
za[npts]: float za deg
azerr[npts]: float azerr asecs
zaerr[npts]: float zaerr asecs
tit: string. title.. keg 'Gain K/Jy'
KEYWORDS:
ticklen : float ticklen in asecs. def: 5
hard : if set then use hardcopy settings
names[npts] : string list of source names used. print left of image.
indlist[npts]: long index into names array for each point in az,za
totab : if set then donot print avg/rms table by za
tabtit : title for table
colar[] : long color indices to use
ln : long line number (1..) to start the notes..
nmoff : long offset line to start printing the source names
DESCRIPTION:
make a 2d arrow plot of azimuth, za error versus az,za. This is
normally used to plot the residual fits of the pointing model. The
length of each arrow will be proportional to the pointing error at
that az,za. The direction of the arrow will be that of the error.
By default a table of the errors in 5 degrees steps is printed at
the bottom of the plot. (it may not show up on the screen version but
it will be there in the postscript file).
(See /pkg/rsi/local/libao/phil/gen/arrpltazzaerr.pro)
NAME:
avgrms - avgerage then compute rms for data set
SYNTAX: avgrms,toAvg,d,rmsAr,allen=allen
ARGS:
toAvg[navg]: long samples to average before computing rms.
If i, greater than 1 then the rms will be computed
on each dataset.
d[m,ncol]:float average along m and then compute rms. Do this for each ncol
KEYWORDS:
allen: if set then return the allen deviation rather than
the standard deviation.
it is y[i]=sqrt(.5*mean((d[i]-d[i+1])^2))
where d[x] has first been averaged by toavg[i].
RETURNS:
rmsAr[navg,ncol]:float return rms info for each dataset
DESCRIPTION:
Average the first index of d and then compute the rms. Do this for
all ncol. If the keyword /allen is set then the sqrt of the allen variance
is computed rather than the standard deviation.
ToAvg determines the number of adjacent samples to average before
computing the deviation. If toavg is an array then the deviation will
be computed for each set of averages.
The info is returned in rmsAr[navg,ncol] where navg is the number of
entries in toavg[].
This routine can be used to check if the rms noise decreases as
1/sqrt(bw*time). To do this, you should normalize the input data to
the mean (or median) value of each column
(since it is DeltAT/T=1./sqrt(b*tau)).
(See /pkg/rsi/local/libao/phil/gen/avgrms.pro)
NAME:
avgrobbychan - compute the robust average by chan for 2d array.
SYNTAX: result=avgrobbychan(d,nsig=nsig,ncnts=ncnts,rms=rms)
ARGS:
d[m,n] : array to compute rms
KEYWORDS:
nsig:float any points nsig beyond mean are ignored in the average.
RETURNS:
result[m]: result[i]= robustmean(d[i,*])
ncnts[m]: number samples used for each mean
rms[m]: rms of each channel (it is not divided by the mean)
DESCTRIPTION:
compute a robust mean by channel of a 2d array. For each channel
compute the rms. Throw out all points above nsig, continue doing
this until no points are thrown out. For each channel return the mean
of the remaining points.
(See /pkg/rsi/local/libao/phil/gen/avgrobbychan.pro)
NAME:
basename - return directory and basename of file
SYNTAX: baseNm=basename(filename,dirNm=dirNm,nmLen=nmLen)
ARGS:
filename: string filename to split into directory and basename
RETURNS:
baseNm: string basename (without the directory). '' will be returned if
there is no base name.
KEYWORDS:
dirNm: string if present then return the directory name here. It will
contain the trailing "\". If there is no directory
name then return ''.
nmLen[2]: long length of directory name:(nmLen[0]) and basename:(nmLen[1]).
if either name is not present then nmLen[i] will have 0.
EXAMPLE:
file='/share/olcor/corfile.01jun06.x101.1'
bnm=basename(file,dirnm=dirnm,nmLen=nmLen)
;
bnm: corfile.01jun06.x101.1
dir: /share/olcor/
nmLen: 13 22
(See /pkg/rsi/local/libao/phil/gen/basename.pro)
NAME: bitreverse - bit reverse data SYNTAX datBr=bitreverse(dat,numBits) ARGS: dat[n]: long data to bit reverse numbits: long number of bits to reverse RETURNS: datBr[n]: long the bit reversed data. DESCRIPTION: Bit reverse the data in array dat. The number of bits to use for the reversal is specified in numBits: EXAMPLE: if a=00000001b then bitreverse(a,4) would return 0001000.
(See /pkg/rsi/local/libao/phil/gen/bitreverse.pro)
NAME:
blmask - interactively create a mask for baseline fit
SYNTAX: istat=blmask(x,y,maskArr,y2=y2)
ARGS:
x[npts]: xaxis array
y[npts]: yaxis array
KEYWORDS:
y2[npts]: overplot second array
RETURNS:
maskArr[npts]: returned mask array with valuse 0,1
istat: 1: created mask, 0: no mask specified by user
DESCRIPTION:
Let the user interactively define a mask to use for fitting. The
x,y data is plotted and then the user is prompted to interactively
build the mask using the mouse. The user is prompted to position the
the mouse to the start and end of each segment that is to be filled
with ones in the mask. The left mouse button is used to specify the
point. The right mouse button will get you out of this loop. The
returned mask will have the value 1 for all the specified segments.
All other values will be set to 0.
If the y2 keyword is provided, then the data in y2 will be over
plotted in red. You might use this when you are trying to create
1 mask for 2 polarizations of data.
EXAMPLE:
istat=blmask(x,y,maskArr)
buttons used: left mark, right quit
1--P1?-- -47.2958
1--P2?-- 9.64367 first section of mask done
2--P1?-- 19.4657
2--P2?-- 22.8821 2nd section of mask done
3--P1?-- 25.1597 right button clicked so it returns.
NOTE:
It is the users responsibility to set the horizontal and vertical scale
prior to calling this routine. The normal interface to this routine is
from the bluser() routine.
SEE ALSO: bluser, cursorsubset.
(See /pkg/rsi/local/libao/phil/gen/blmask.pro)
NAME:
bluser - interactively baseline a function.
SYNTAX: istat=bluser(x,y,coef,mask,yfit,maskst=maskst)
ARGS:
x[npts]: xaxis array
y[npts]: yaxis array
KEYWORDS:
maskst[npts]: if provided, then use this mask to start with.
RETURNS:
coef[] : coefficients from fit
mask[npts]: mask used for fit that was defined by the user.
yfit[npts]: fit evaluatuted at the data points
istat: 1: fit ok, 0 no fit
DESCRIPTION:
bluser lets the user interactively baseline a function. The user passes
in the x,y arrays of data. The routine will pass back the coef's from the
fit, the mask used, and the fit evaluated at the data points. On entry
the main menu is displayed:
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
The user should first adjust the horizontal and vertical scale with
h h1 h2 and v v1 v2. Each time one of these is entered the plot will
be redisplayed with the new limits (do not use commas here..).
When the plot shows the correct limits, enter m to define the mask. The
leftmost mouse button is used to define the starting and ending portions
of the data that will be used for the fit. You can have as many of these
sections as you want. When you are done with the mask, click the right
mouse button.
After defining the mask, you can fit any order polynomial you want.
f 4
would fit a fourth order polynomial and overplot the fit.
p .. plots data - fit. (use v v1 v2 to reset the vertical scale).
You can go back and redefine a new mask if you want and the redo the fits.
q is how you exit the routine.
SEE ALSO: blmask
(See /pkg/rsi/local/libao/phil/gen/bluser.pro)
NAME:
bpcmpbychan - compute a band pass from set of spectra
SYNTAX: bpc=bpcmpbychan(d,
deg=deg,nsig=nsig,nlsig=nlsig,chnNsig=chnNsig,$
flatTsys=flatTsys,
maxFitloop=maxFitLoop,$
chanRms=chanRms,resAll=resAll,dfit=dfit,$
gdpntperchan=gdpntperchan,nloop=nloopbychan,$
indBadChan=indBadChan,indgdchan=indgdchan,zeromean=zeromean,
tsysScl=tsysScl
ARGS:
d[nchn,npts]: float input data to process
KEYWORDS:
deg: int degree of polynomial fit along each channel. default=1
nsig: float The clipping level (in sigmas) to use for the fit
residuals. Any value greater than this is not included
in the fit. The default is 3 sigma.
nlsig: float if nlsig is provided, then channels that have large values
of sigma/median will have there points to exclude recomputed
with nlsig rather then nsig. A good value would be 2.
By default this is not done.
chnNsig: float : If indBadChan, or indGdChan is requested, use
chnNsig to determine which channels are good:
(rms/sigma le chnNsig), and which channels are bad:
(rms/sigma gt chnNsig). The default for chnNsig is 3 sigma.
flatTsys: if set, then divide each spectra by the median
value. This will allow continuum sources to be included
in drift scan mode.
maxFitLoop:int The maximum number of times to iterate on the fit for
a single channel. The default is 20.
zeromean : if set, the the input data has zero mean, donot divide
byte the mean;
RETURNS:
bpc[nchn] :float the bandpass to use for the bandpass correction.
chanRms[nchn] :float the rms/mean computed for each channel.
resAll[nchn,npts] :float the fit residuals for each point (y-yfit) in
units of sigma of that channel.
dfit[nchn,npts]:float the 2d fit of the data to the polynomials by chan.
gdPntPerChan[nchn]:long the number of good points used in each channel.
indGdChan[] :long The indices (0 based) for all channels who had
an rms/mean less than chnNSig
indBadChan[nchn] :long The indices (0 based) for any channels who had
an rms/mean greater than chnNsig.
nloop[nchn]:long The number of fitting iterations that where done
on each channel.
tsysScl[npts] :float If flatTsys was set, then this will contain
1/medianTsysBySpectra for each spectra. have the
tsys scaling factor for each spectra.
DESCRIPTION:
Compute a bandpass correction for a set of data.
It works with an array of npts spectra each having nchn channels.
The algorithm is:
1. If tsysFlat is selected then
compute the median over freq for each spectra and call it Tsys[npts]
Use this to flatten each time point so that continuum sources
do not skew the data and can be included in the statistics.
If TsysFlat is not selected, then set Tsys[npts] to 1.
2. for each channel take all of the time points for that channel:
y=d[ichn,*]
a. Normalize this channel to the average tsys over time:
y[npts]=y[npts]/Tsys[npts]
b. define all npts to be good points.
c. fit a polynomial of order deg to the good points,
coef=poly_fit(x[gdpnts],y[gdpnts],deg)
d. compute the residuals over all points using the coef from c.
res=y-coef(x,coef)
e. find all points that are less than nsig times the fit sigma.
f. If there are fewer points than we started with in c, goto c, if
not, then we are done with the channel
g. store the following:
- chanRms[ichn]= rms/mean(last fit)
- dfit[ichn,*] = last fit along this channel
- gdpntPerChan get number of points we were left with after the fit.
- nloop[ichn] is the number of times we looped on the fitting.
3. When done with all channels, if indBadChan or indGdChan is provided,
call meanrob( robust mean) with the chanrms[] array. It will
find all of the channel sigmas that stick out less than chnNsig.
4. Compute the average bandpass correction by averaging the fit array
over all time samples: bpc=total(dfit[nchn,npts],2)/npts
6. The function returns bpc[nchn]
Some Notes on the returned data:
1. The data is flattened in the time direction. The results:
bpc, dfit,chanRms,resAll are relative to this dataset. With this
processing, continuum sources will not have large residuals (since
they were flattened by the mean power of that sample). If you overplot
the bpc with the input data, they will differ by a scaling factor. This
will be corrected for when the cal is applied.
2. resAll is in units of rms of the flattened channel: (y-yfit)/sigma.
3. chanRms[nchn] has been divided by the mean of each channel. This
flattens the bandpass edges. You should be able to predict what
this value should be from the 1./sqrt(bw*tau)
EXAMPLE:
Suppose you have spectra spc[1024,300] and a calOnOff[1024,2]
Processing would be:
bpc=bpcmpbychan(spc,chanRms=chanRms,indgdChan=indgdchan)
;
; don't use 100 channels on each edge, plus all the chanRms gt 3 sigma
; for scaling to kelvins.
;
mask=lonarr(1024)
mask[indgdchan]=1.
mask[0:99]=0
mask[1024-99:*]=0
ind=where(mask eq 1)
ngood=n_elements(ind)
bpc=bpc/(total(bpc[ind])/ngood) ; bp now normalized to unity
for i=0,nsmp-1 do spc=spc[*,i]/bpc
;
; now scale to kelvins
;
caldeflTP=total(calOnoff[ind,0]-calonOff[ind,1])/ngood ; use same good chan
corToK=calKelvins/calDeflTp
spc=spc*corToK spectra now in kelvins.
; I've left the steps separate for illustration. You could combine the
; calscaling and bandpass corretion in one multiply.
(See /pkg/rsi/local/libao/phil/gen/bpcmpbychan.pro)
NAME:
bytesleftfile - return the unread bytes left in a file
SYNTAX: bytesLeft=bytesleftfile(lun,bytereq=bytereq,chktime=chktime,$
maxloop=maxloop,cursize=cursize)
ARGS:
lun: int lun to already open file that we should check.
KEYOWRDS:
bytereq: long Wait until at least this number of bytes is available.
use the chktime keyword to determine how often to check
chktime : float number of seconds to delay between file size checks. This
is used if bytereq is specified. The default is 1 second.
maxloop : long max number of times to loop if bytereq option used.
The default is 999999
RETURNS:
bytesLeft: LONG returns the unread bytes in the file
cursize : LONG current number of bytes in file
DESCRIPTION:
The routine will return the number of unread bytes in a file.
It takes 5 to 10 milliseconds to check the file size.
The bytereq option will wait until at least bytereq bytes are available
in the file. The chktime option specifies how often to check the file size.
The maxloop option tells how many times to check the file before quitting.
EXAMPLE:
openr,lun,'/share/olcor/corfile.15jun02.x101.1',/get_lun
nbytes=bytesleftfile(lun)
Wait till there are at least 5k bytes. Check every 10 secs and loop
a maximum of 60 times:
nbytes=bytesleftfile(lun,bytereq=5000,chktime=10,maxloop=60)
(See /pkg/rsi/local/libao/phil/gen/bytesleftfile.pro)
NAME:
calget - return the cal value given hdr,freq.
SYNTAX: stat=calget(hdr,freq,calval,date=date,swappol=swappol)
ARGS:
hdr: {hdr} header holding at least hdr.iflo
freq[n]: float freq in Mhz for cal value.
KEYWORDS:
date[2]: int [year,daynumber] epoch for cal value.
Default is the date in the header.
swappol: If set then swap the polA polB calvalues. This
can be used to correct the 1320 hipass polarization
problem, or to compensate for a xfer switch being
used.
RETURNS:
calval[2,n]: float cal values in deg K for polA,polB. if alfa then it returns
calval[2,7,n]
stat: int -1 error, 0 got the values ok.
DESCRIPTION:
Return the cal values in degrees K for the requested freq. The hdr can be
a correlator or ri header (as long as it includes hdr.iflo). This
routine always returns 2 values: pola, and polB.
The calvalues for the receiver in use are looked up and then the
values are interpolated to the observing frequency.
NOTE:
Some cals have measurements at a limited range of frequencies (in some
cases only 1 frequency). If the requested frequency is outside the range
of measured freqeuncies, then the closest measured calvalue is used
(no extrapolation is done).
hdr should be a single element rather than an array.
This routine extracts info from the header and then calls calget1().
SEE ALSO:calget1, calval, calinpdata, corhcalval
(See /pkg/rsi/local/libao/phil/gen/calget.pro)
NAME:
calget1 - return the cal value given rcvr,frq,type.
SYNTAX: stat=calget1(rcvrNum,caltype,freq,calval,date=date,hybrid=hybrid,
fname=fname,swappol=swappol)
ARGS:
rcvrNum: int receiver number 1..16 (see helpdt feeds).
calType: int type of cal used 0 through 7.
0-lcorcal,1-hcorcal,2-lxcal,3-hxcal,
4-luncorcal,5-huncorcal,6-l90cal,7-h90cal
freq[n]: float freq in Mhz for cal value.
KEYWORDS:
date[2]: int [year,daynumber] data for calvalue. Default is the
current date.
hybrid: if set then a hybrid was in use and the cal values should
be averaged together.
fname:string use an alternative filename for cal data.
swappol: if set then swap the polA, polB calvalues on return.
This can be used to correct for the 1320 hipass
polarization cable switch or a xfer switch.
RETURNS:
calval[2,n]: float .. calValues in deg K for polA,polB
Note if this is alfa, then [2,7] value are returned
for the 7 pixels..
stat: int .. -1 error, 0 got the values ok.
DESCRIPTION:
Return the cal values in degrees K for the requested reciever, caltype,
and frequency. This routine always returns 2 (or 2x7) values: pola, and polB.
The calvalues for the receiver in use are looked up and then the
values are interpolated to the observing frequency.
EXAMPLES:
Get the cal values for lbw (rcvrNum=5) using the high correlated cal
(caltype=1) at 1400. Mhz.
stat=calget1(5,1,1400.,calval)
NOTE:
Some cals have measurements at a limited range of frequencies (in some
cases only 1 frequency). If the requested frequency is outside the range
of measured freqeuncies, then the closest measured calvalue is used
(no extrapolation is done).
If you have a datatking header, you can use calget(). It will
take the rcvrNum, and caltype, from the header.
SEE ALSO:calget, calval, calinpdata, corhcalval
(See /pkg/rsi/local/libao/phil/gen/calget1.pro)
NAME:
calInpData - input cal data for rcvr/calType.
SYNTAX:
istat=calInpData(rcvNum,calNum,calData,fname=fname,date=date)
ARGS:
rcvNum: 1 thru 16. receiver to use (same as hdr.iflo.stat1.rfnum).
calNum: 0 thru 7. cal Type to use (same as hdr.iflo.stat2.calType).
the values are:0-lcorcal,1-hcorcal,2-lxcal,3-hxcal,
4-luncorcal,5-huncorcal,6-l90cal,7-h90cal
KEYWORDS:
fname: to specify an alternate data file with cal values.
The default file is aodefdir() + 'data/cal.dat{rcvNum}
date : [year,daynum] .. if specified the data when you want the
cals for..default is most recent.
RETURNS:
istat: 1 ok, -1 bad file/rcvnum,data , -2 bad calnum
calData: return data in structure:
calData.calNum - calNum
calData.numFreq - number of freq entries
calData.freq[numFreq] - for each cal value (in Mhz)
calData.calA[numFreq] - polA cal value in Kelvins
calData.calB[numFreq] - polB cal value in Kelvins
DESCRIPTION:
Input the cal data for all the frequncies of a particular receiver (rcvNum)
and cal type (calNum). The calNum and rcvNum can be extracted from the
headers with iflohrfnum() and iflohcaltype ().
The default datafile is aodefdir() + 'data/cal.dat{rcvNum} (aodefdir() is
a function that returns the root of the aoroutines). The keyword
fname allows you to specify an alternate file. The file format is:
- col 1 # is a column
- data is free format , column oriented
- nd1 is noise diode 1, nd2 is noise diode2, H-highcal,L-lowcal,
- A is polA, B id polB, ndxx-->A/B implies that diode N feeds pol X
freq nd1H->A nd1L->A nd1H->B nd1L->B nd2H->A nd2L->A nd2H->B nd2L->B
The mapping of cal type to diodes used is:
corCal diode 1-> polA, diode 1-> polB
uncorCal diode 1-> polA, diode 2-> polB
xCal diode 1-> polB, diode 2-> polA
90cal diode 2-> polA (with 90deg phase shift), diode 2-> polB
This routine is called automatically by corhcalval and calget().
How the different cal routines vary:
calinpdata() inputs the data from disc. You must specify the rcvrnum,calnum.
It defaults to the current date. It loads a table in common
but it does not interpolate or compute a calvalue.
calval() Pass in the frequency and the caldata array input via
calinpdata(). It will interpolate the frequency and compute the
cal value.
calget() You supply a header and a frequency. The routine figures out the
caltype,rcvrNum, and date from the header and calls calinpdata()
and calval(). It returns the cal values.
corhcalval() You specify the correlator sbc header (eg b.b1.h). It will
compute the frequency and then call calget(). It returns the
cal values.
NOTE:
The following receivers will always return a single calNum independent
of what is requested (since they only have 1 single cal).
rcvnum nam calnumreturned
12 sbn 0 low cor cal
3 610 0 low cor cal
6 lbn 0 low cor cal
100 ch 5 h uncorcal
SEE ALSO: corhcalval, calget
(See /pkg/rsi/local/libao/phil/gen/calinpdata.pro)
NAME:
calmfit - fit curves to cal on/off data
SYNTAX: istat=calmfit(d,nsteps,nloops,fitI,fitIAvg,fitISum,verb=verb,$
bw=bw,indAvgA=indAvgA,indAvgB=indAvgB,$
masktoUseA=masktouseA,masktouseB,avgstop=avgstop,
dopol=dopol, _extra=e,
ARGS:
d[n]: {mcal} : data input from mcalinp. Uses
(spcCalOn - spcCalOff)/spcaloff.
nsteps: int ; number of 100 Mhz steps to move through the entire band
nloops: int ; number of entire band was measured
fitI[nloops] :{ } ; results of the fitting each loop
fitIAvg :{ } ; fit to the average of the loops
fitISum :{ } ; summary of the average fit
KEYWORDS:
verb :int ; passed to corblauto for plotting
bw :float ; bandwidth in Mhz each measurement (25Mhz)
indAvgA[ma]:long ; indices of nloops to use when averaging polA
indAvgB[mb]:long ; indices of nloops to use when averaging polB
maskToUseA[mm,j]:long ;mask array polA for each fit. if j=0 then use
the same mask for each fit.1==> use,0=-> ignore
maskToUseB[mm,k]:long ;mask array polB for each fit. if k=0 then use
the same mask for each fit.1==> use,0=-> ignore
dopol : int if provided then only process the pols specified
by dopol. 0=polA,1=polB. default is both pols
avgStop: if set then stop after plotting each avg
so the user can look at the plot output.
_extra: ; passed to corblauto().
deg=deg (polynomial fit order)
fsin=fins (harmonic fit order)
DESCRIPTION:
mcalinp() inputs a set of calon/calOff -1 spectra. CalMfit will
fit a function to the spectra covering the entire receiver bandpass.
Each entry in d holds 100 Mhz worth of data. nsteps of these will
cover the entire bandpass and there are nloops copies. Since there can
be multiple passes through the entrie band, d has dimensions of
d[nlags,nsteps,nloops]).
corblauto is used to fit to the spectra calon/caloff-1 for a complete
pass through the receiver band. The fit is repeated nloops time (once for
each pass through the band). For each fit the following data is returned:
help,fitI,/st
** Structure <8445c7c>, 7 tags, length=329324, data length=329324, refs=2:
NP LONG 20480 ; number pnts across band
FRQSTP FLOAT 0.0976562 ; freq step Mhz
FRQMIN FLOAT 4000.00 ; minFreq
FRQMAX FLOAT 6000.00 ; max freq
FITI STRUCT -> Array[1] ; fit coef.
YFIT FLOAT Array[20480, 2]; fit evaluated at freq points
MASK FLOAT Array[20480, 2]; mask used for fit (non
zero values).
A robust average by channel is taken for the nloops through the band. A
fit to this average spectrum is then done and returned in fitIAvg. By default
all of the passes through the data are used when computing the average. The
keywords indAvgA, indAvgB let you specify which passes through the data should
be included when making the averages.
fitISum holds summary info for the average:
** Structure <8438ab4>, 11 tags, length=163896, data length=163896, refs=1:
NP LONG 20480
NLAGS LONG 256
NSBC INT 4
NSTEPS INT 20
FRQSTP FLOAT 0.0976562
FRQMIN FLOAT 4000.00
FRQMAX FLOAT 6000.00
NLOOPS LONG 7
AVGD FLOAT Array[20480, 2]
USEAVGA INT Array[7]
USEAVGB INT Array[7]
It also includes the averaged data in .avgD
NOTE: if nloops eq 1 then fitiavg is not returned, since it is the same
as fitI. fitISum will be returned.
(See /pkg/rsi/local/libao/phil/gen/calmfit.pro)
NAME:
calonofffind - find all of the cal on/offs in file
SYNTAX: numfound=calonoffind(lun,sl,indon)
ARGS:
lun: int lun for file to search
sl[]: {getsl} scan list structure (returned from getsl(sl).
RETURNS:
indon[numfound]: int indices into sl for the start of each cal on of
a pair.
DESCRIPTION:
Find all of the calonoff pairs in a file. Return the indices into sl
for the calon scans. For a pair to be included the calon scan must be
immediately followed by a cal off scan.
EXAMPLES:
sl=getsl(lun)
numfound=calonoffind(lun,sl,indon)
for i=0,numfound-1 do begin
print,posscan(lun,sl[ind[i]].scan,1,sl=sl) ; position to start calon
istat=corcalonoff(lun,retdat) ; process each onoff pair
endfor
SEE ALSO:
getsl,corcalonoff
(See /pkg/rsi/local/libao/phil/gen/calonofffind.pro)
NAME:
calval - return the cal value for a given freq.
SYNTAX:
istat=calval(freqReq,calData,calV,hybrid=hybrid,swappol=swappol)
ARGS:
freqReq[n]: float frequency in Mhz for cal
calData: {calData} already input via calInpData()
KEYWORDS:
hybrid: if keyword set, then average the polA,polb values
together (used when linear receivers are converted to circular
by a hybrid after the dewar).
swappol: if set then swap the polA, polb calvalues in the calV array
on return. This can be used to correct for the
1320 hipass polarization cable switch or the use of
one of the xfer switches.
alfapixnum: if this is alfa data, then the pixel number (0 thru 6) of
the alfa pixel to use. If the cal values are alfa data
and this is not specified, use pixel 0.
RETURNS:
istat: 1 ok within range, 0 outside range used edge,-1 all zeros
calV[2,n]: float array of [2] floats holding interpolated cal values for
polA and polB. If this is alfa then the calV is
dimensioned 2,7 for the 7 pixels.
DESCRIPTION:
Interpolate the cal value to the requested frequency. The calData
should have already been read in with calInpData. If the requested
frequency is outside the data range, return the cal values at the
edge (no extrapolation is done). The data is returned in an array
of two by n values.
The normal way to get cal values is via corhcalval() or calget().
They call this routine.
SEE ALSO:corhcalval, calget, calinpdata.
(See /pkg/rsi/local/libao/phil/gen/calval.pro)
NAME:
cataloginp - input a pointing catalog
SYNTAX: nsrc=cataloginp(file,format,retdata,comment=comment)
ARGS:
file :string filename of catatlog
format :int format for catalog:
1: srcname hh mm ss dd mm ss
2: srcname hhmmss ddmmss
retdata[]:{srccat} return data here
in retdata.
KEYWORDS:
comment : string comment characters for catatlog.def:#
DESCRIPTION:
Read in all of the source names and positions catalog specified by
file
The returned srccat array will contain:
help,retdat,/st
** Structure CATENTRY, 6 tags, length=52:
NAME STRING '' source name
RA FLOAT Array[3] hh mm ss.ss
DEC FLOAT Array[3] dd mm dd.dd (alway positive)
DECSGN INT 0 +/- 1 sign of declination
RAH DOUBLE 0.0 ra in hours (includes sign)
DECD DOUBLE 0.0 dec in hours (includes sign)
EOL STRING string dec to end of line
(See /pkg/rsi/local/libao/phil/gen/cataloginp.pro)
NAME:
chebeval - evaluate chebyshev polynomial
SYNTAX: y=chebeval(a,b,coef,x)
ARGS:
a : double min value of xrange used for fit.
b : double max value of xrange used for fit.
coef[m]: double coefficients from fit.
x[n: double xvalues where polynomial should be evaluated.
RETURNS:
y[n]: fit evaluated at the requested x values.
DESCRIPTION:
Evaluate a chebyshev polynomial at the requested x values. These
values should be within the x values used for the fit. The a,b
parameters are the min,max x values used in the fit.
SEE ALSO: chebfit()
(See /pkg/rsi/local/libao/phil/gen/chebeval.pro)
NAME:
chebfit - chebyshev polynomial fit to data
SYNTAX:coef=chebfit(x,y,deg,merr=merr,$
yfit=yfit,rangex=rangex,singular=singular,chisq=chisq,$
covar=covar)
ARGS:
x[n]: float/double independent var.
y[n]: float/double measured dependent variable
deg: int deg of fit (ge 1)
KEYWORDS:
merr[n]: float/double measurement errors for y.default is uniform
RETURNS:
coef[deg+1]: coefs from fit
KEYWORDS RETURNS:
yfit[n]: float or double . fit evaluated at x locations.
rangex[2]: float/double min max values of x used for fit.
these were used to map the x axis
into [-1,1] for the fit.
chisq : float/double chisqr from svdfit
covar[] : float/double covariance matrix from svdfit().
singular: int number of singular points found (see svdfit()).
DESCRIPTION:
Do a chebyshev polynomial fit of order deg to the x,y data. Merr
are the measurement errors (see idl svdfit routine).
Return the coefs for the fit as well as the mapping of the xrange
into [-1,1].
SEE ALSO:
chebeval() to evaluate the coef.
NOTE:
The fitting function svdcheb() is contained in this file. If the
routine gives an error that it cannot find svdcheb() just compile
this routine explicitly (.compile chebfit).
(See /pkg/rsi/local/libao/phil/gen/chebfit.pro)
NAME:
checkkey - check if any keys have been pressed
SYNTAX: key=checkkey(first=first,all=all,wait=wait,noflush=noflush)
ARGS : NONE
KEYWORDS:
first: if set then return the first character. default is the last char.
all: if set then return all characters. default is the last char.
wait: if set then wait for at least one character
noflush: if set and first is set, then don't flush the other characters
in the buf.
RETURNS:
key : character or characters entered.
DESCRIPTION:
checkkey will check to see if any keys have been pressed. It will return
with the last key waiting in the input buffer or '' if nothing was there.
The wait keyword will cause the routine to wait for at least 1 keypress.
The first keyword will return the first rather last key of any string.
The all keyword will return all keys entered. By default all characters
waiting in the input buffer are read. If the noflush and first keyword are
set then only the first char is returned. The other chars can be read at
a later time.
(See /pkg/rsi/local/libao/phil/gen/checkkey.pro)
NAME:
chkswaprec - check hdrlen to see if the record needs to be swapped
SYNTAX: swap=chkswaprec(stdhdr)
ARGS:
stdhdr: {hdrStd} standard portion of data header
RETURNS:
swap : int 1 --> need to swap,0 --> no need to swap
DESCRIPTION
Check to see if the record needs its data swapped because of
big/little endian differences. The user passes in the standard header
and the routine checks that the abs(hdrstd.hdrlen) < 65535. If the
number i larger than this then the record needs to be swapped (headers
are never larger than 65535.
(See /pkg/rsi/local/libao/phil/gen/chkswaprec.pro)
NAME: cmpbflytwiddle - compute butterfly twiddle factors SYNTAX - cmptwiddle,lenfft,nbfly,ntwiddleAr,bflytwiddle ARGS: lenfft: long length of fft RETURNS: nbfly: int number of butter fly stages. ntwiddleAr[nbfly]: step size, number of twiddle values this butterfly bflytwiddle[2,lenfft,nbfly]: float hold the twiddle values for this butterfly for each butterfly i, the first ntwiddleAr[i] entries of bflytwiddle[lenffit,i] will hold the values. DESCRIPTION: Compute the fft twiddle factors for lenfft. This routine is for debugging, not for speeding up the fft computation. For each butterfly stage the m uniq twiddlefactors for that stage are returned in bflyTwiddle[0:m-1,ibutterfly]. NtwiddleAr[nflys] holds how many uniq twiddle factors there are at each stage.
(See /pkg/rsi/local/libao/phil/gen/cmptwiddle.pro)
NAME: colorinfo - return info on current color setup. SYNTAX: colorinfo ARGS: DESCRIPTION: Queries idl on the current color setup. Output sent to stdout.
(See /pkg/rsi/local/libao/phil/gen/colorinfo.pro)
NAME:
condblevels - compute db contouring levels for a map
SYNTAX: levels=condblevels(map,nlevels,dbstep,maxval)
ARGS:
map[m,n]: float data to to compute levels for.
nlevels : int number of levels requested
dbstep; float dbstep between levels
RETURNS:
levels[nlevels]: float holding the values to use to mark the contours
maxval : float the maximum value in map. The contours are relative
to this value.
DESCRIPTION:
Compute nlevels space dbstep apart from the maximum value in map.
(See /pkg/rsi/local/libao/phil/gen/condblevels.pro)
NAME:
contourph - phil's interface to idl contour routine.
SYNTAX: contourph,d,numcontours,dbstep,maxval,levels,axes=axes,$
fill=fill,_extra=e
ARGS:
d[nx,ny] : float data to contour.
numcontours: int number of contours to plot.
dbstep: float db step between contour levels.
RETURNS:
maxval: float db scaling is relative to this value.
levels[numcontours]: float the contour levels in dbs relative to maxval.
KEYWORDS:
axes[4]: float axes label values: [minx,maxx,miny,maxy]. If not
supplied then 0:nx-1,0:ny-1 are used.
fill: int if set then fill the contours with colors.
e: extra keywords sent to the contour routine.
DESCRIPTION:
contourph interfaces to the idl contour routine. It will scale
a dataset to dbs relative to the maximum value. The number of levels and
db step size is input by the user. Annotation can be entered via the
_extra=e keyword.
EXAMPLE:
let bmmap[120,50] be a baselined beammap dataset of 20 arcminutes in az and
15 arcminutes in za. To contour the data with 12;contours at 2db steps:
axes=[-10.,10.,-7.5,7.5]
contourph,bmmap,12,2,maxval,levels,axes=axes,xtitle='az [amins]',$
ytitle='za [amins]'
(See /pkg/rsi/local/libao/phil/gen/contourph.pro)
NAME: cor3lvlstokes - to 3 level correction for stokes data SYNTAX: acfCor=cor3lvlstokes(acfIn,nlags,bias,ntodo,double=double,ads=ads ARGS: acfIn[nlags,4,ntodo]: uncorrected acf's nlags: int number of lags in each auto or cross correlation. ntodo: int number of sets of aa,bb,ab,ba to do. KEYWORDS: double: if set then do computation in double precision RETURNS: acfcor[nlags,4,ntodo]: float/double corrected correlation functions ads[2,ntodo]: float/double digThreshold/sigma for the auto correlations. DESCRIPTION: 3 level correct the auto and cross correlations. The data should be ordered AA[nlags],BB[nlags],ab[nlags],ba[nlags] NOTE!!! THIS IS STILL BEING TESTED...
(See /pkg/rsi/local/libao/phil/gen/cor3lvlstokes.pro)
NAME:
corbychan - auto/xcorrelate dynamic spc by chan.
SYNTAX: cormat=corbychan(spc1,spc2,avg1=avg1,avg2=avg2,rms1=rms1,rms2=rms2)
ARGS:
spc1[nchn,mtime] : float first dynamic spectra
spc2[nchn,mtime] : float 2nd dynamic spectra (not needed if auto
correlation done
KEYWORDS:
RETURNS:
cormat[nchn,nchn]: float correlation matrix returned
avg1[nchn] : float robust average over time of spc1
avg2[nchn] : float robust average over time of spc2
rms1[nchn] : float rms by channel of spc1
rms2[nchn] : float rms by channel of spc2
DESCTRIPTION:
Compute the correlation matrix for a set of dynamic spectra. It computes
the correlation between every two pairs of channels in one set of
dynamic spectra (if only spc1 entered) or between channels in two
sets of dynamic spectra (if spc1,spc2 entered).
The computation is:
Let:
nchn=the number of freq channels
mtm =the number of time samples
for spc 1 or 2 let:
savg[nchn]= mean(spc[,]) averaging over the time samples
srms[nchn]= rms(spc) computed along the time axis of each chan
s[i,j] = (spc[i,j] - sAvg[i]))/srms[i] .. remove mean normalize to sigma
cormat[i,j]= sum_k(s1[i,k]*s2[j,k])/nsmp
The routine uses robust averaging and rms (outliers are not used).
(See /pkg/rsi/local/libao/phil/gen/corbychan.pro)
NAME: corinit - initialize to use the idl correlator routines. SYNTAX: @corinit DESCRIPTION: call this routine before using any of the correlator idl routines. It sets up the path for the idl correlator directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/corinit.pro)
NAME: corinit1 - initialize to use the idl correlator routines (no lut load). SYNTAX: @corinit DESCRIPTION: call this routine before using any of the correlator idl routines. It sets up the path for the idl correlator directory and defines the necessary structures. It calls geninit1 instead of geninit. The color table ldcolph call is not made. This speeds things up for remote observers.
(See /pkg/rsi/local/libao/phil/gen/corinit1.pro)
NAME: covarnorm - normalize the covariance matrix SYNTAX: covar=corvarnorm(covar) ARGS: covar[]: float covariance matrix. DESCRIPTION: Normalize a covariance matrix to have unit diagonols.
(See /pkg/rsi/local/libao/phil/gen/covarnorm.pro)
NAME: cp - read cursor position after button press. SYNTAX: cp,x=x,y=y,/up,/down,/change,/nowait,/data,/device,/normal ARGS : none KEYWORDS: up : if set then wait for button up event down : if set then wait for button down event change: if set then wait for button change event nowait: if set then return immediately with the current x,y values data : if set then return x,y in data coordinates (default) device: if set then return x,y in device (screen) coordinates normal: if set then return x,y in normalized device coordinates (0.,1.) x : float return x value here y : float return y value here DESCRIPTION: This routine calls the idl cursor routine. The keywords are the same as cursor. By default the routine waits for the cursor to be depressed (or the cursor is already down). It then reads the x,y coordinates, prints them out, and then returns. If you want to loop reading the cursor n times with a button push on each step of the loop , then you need to use /up or /down. EXAMPLE: plot,findgen(100) cp .. user clicks left button at desired position on plot. 24.0208 23.2295 .. x,y positions printed out. NOTE: If the window system is set so that the window focus follows the cursor, then you must make sure that the cursor is in the idl input window before you enter the command cp. SEE ALSO: The idl routine cursor for a descriptoin of the keywords.
(See /pkg/rsi/local/libao/phil/gen/cp.pro)
NAME:
cumfilter - cumfilter routine from carl heiles.
SYNTAX: cumfilter,data,range,limit,indxgood,indxbad,countbad,
median=median,correct=correct
ARGS:
data[n]: data to filter
range : long number of channels in center of distribution to use
to compute limits.
limit : float used to compute the min, max values to keep.see below
KEYWORDS:
median : if set then remove a median filtered version of data
before cumfiltering.
correct: if set then replace the bad data. If median not set then
replace bad data with the median of the entire dataset. If
median set then replace bad data with the median filtered
version of the original data (filter len=16).
RETURNS:
indxgood[]: long indices into data for ok data .
indxbad[] : long indices into data for data that should be filtered out.
countbad : long number of elements in indxbad.
DESCRIPTION:
This is carl heiles' cumfilter routine. The basic idea is how to
define reasonable limits for clipping bad data when there may be
large outliers and you don't know the distribution. The algorithm sorts
the data and then uses a range about the center to define the limits:
original data:
data[n]
sorted data:
sdata=sort(data)
find the hi , low value of the data set "range" elements about the
center.
low =sdata[n/2-range/2]
high=sdata[n/2+range/2]
set the limits of "ok" DATA to be "limit" times this low,high.
min=limit*low
max=limit*high
all points with values between min,max are ok..
If the median keyword is set, we median filter (filter length=16) the
data set and remove this filtered version from the original data before
cumfiltering ( eg: data=data - median(data,16)).
If the correct keyword is set, then the bad data points are replaced by
the median of the data set (no median keyword) or by the median filtered
value for that point (if median keyword is set).
EXAMPLE:
Compute the total power of a spectrum using cumfiltering.
d[1024] are the spectral channels.
range=1024/4
limit=3
cumfilter,d,range,limit,indxgood,indxbad,countbad
tp=mean(d[indxgood])
(See /pkg/rsi/local/libao/phil/gen/cumfilter.pro)
NAME:
CURVEFITPP-phils version of curvefit with a few bug fixes.
PURPOSE:
Non-linear least squares fit to a function of an arbitrary
number of parameters. The function may be any non-linear
function. If available, partial derivatives can be calculated by
the user function, else this routine will estimate partial derivatives
with a forward difference approximation.
CATEGORY:
E2 - Curve and Surface Fitting.
CALLING SEQUENCE:
Result = CURVEFIT(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE,$
CHISQ=CHISQ,FLAMBDASTEP=FLAMBDASTEP,covar=covar,
cfplot=cfplot,cfparms,trouble=trouble,halfass=halfass)
INPUTS:
X: A row vector of independent variables. This routine does
not manipulate or use values in X, it simply passes X
to the user-written function.
Y: A row vector containing the dependent variable.
Weights: A row vector of weights, the same length as Y.
For no weighting,
Weights(i) = 1.0.
For instrumental (Gaussian) weighting,
Weights(i)=1.0/sigma(i)^2
For statistical (Poisson) weighting,
Weights(i) = 1.0/y(i), etc.
A: A vector, with as many elements as the number of terms, that
contains the initial estimate for each parameter. IF A is double-
precision, calculations are performed in double precision,
otherwise they are performed in single precision. Fitted parameters
are returned in A.
KEYWORDS:
FUNCTION_NAME: The name of the function (actually, a procedure) to
fit. IF omitted, "FUNCT" is used. The procedure must be written as
described under RESTRICTIONS, below.
ITMAX: Maximum number of iterations. Default = 20.
ITER: The actual number of iterations which were performed
TOL: The convergence tolerance. The routine returns when the
relative decrease in chi-squared is less than TOL in an
interation. Default = 1.e-3.
CHI2: The value of chi-squared on exit (obselete)
CHISQ: The value of reduced chi-squared on exit
NODERIVATIVE: IF this keyword is set THEN the user procedure will not
be requested to provide partial derivatives. The partial
derivatives will be estimated in CURVEFIT using forward
differences. IF analytical derivatives are available they
should always be used.
pjp flambdastep: This method moves between a steepest descent and the
inverse hessian method (using the curvature matrix to
compute the answer when we are close to the solution).
pjp trouble : 0 converged ok
-1 chisq infinite
-2 flambdacount > 30 * 10/flambdastep
-3 iteration > itermax .default 20
-4 alpha/c not finite
pjp nostop : if set then don't stop if alpha/c not finite, just return
with trouble set
pjp halfass : 0..1. multiply step by this amount to slow down motion
pjp cfparms : if set then print parameters input to fit
pjp cfplot : 0 no plot,
1 plot no wait
2 plot wait at last on fit
3 plot wait each one
OUTPUTS:
Returns a vector of calculated values.
A: A vector of parameters containing fit.
OPTIONAL OUTPUT PARAMETERS:
Sigma: A vector of standard deviations for the parameters in A.
COMMON BLOCKS:
NONE.
SIDE EFFECTS:
None.
RESTRICTIONS:
The function to be fit must be defined and called FUNCT,
unless the FUNCTION_NAME keyword is supplied. This function,
(actually written as a procedure) must accept values of
X (the independent variable), and A (the fitted function's
parameter values), and return F (the function's value at
X), and PDER (a 2D array of partial derivatives).
For an example, see FUNCT in the IDL User's Libaray.
A call to FUNCT is entered as:
FUNCT, X, A, F, PDER
where:
X = Variable passed into CURVEFIT. It is the job of the user-written
function to interpret this variable.
A = Vector of NTERMS function parameters, input.
F = Vector of NPOINT values of function, y(i) = funct(x), output.
PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
PDER(I,J) = DErivative of function at ith point with
respect to jth parameter. Optional output parameter.
PDER should not be calculated IF the parameter is not
supplied in call. IF the /NODERIVATIVE keyword is set in the
call to CURVEFIT THEN the user routine will never need to
calculate PDER.
PROCEDURE:
Copied from "CURFIT", least squares fit to a non-linear
function, pages 237-239, Bevington, Data Reduction and Error
Analysis for the Physical Sciences. This is adapted from:
Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear
Parameters", J. Soc. Ind. Appl. Math., Vol 11, no. 2, pp. 431-441,
June, 1963.
"This method is the Gradient-expansion algorithm which
combines the best features of the gradient search with
the method of linearizing the fitting function."
Iterations are performed until the chi square changes by
only TOL or until ITMAX iterations have been performed.
The initial guess of the parameter values should be
as close to the actual values as possible or the solution
may not converge.
pjp Notes: This method moves between a steepest descent (move in the
direction of decreasing chisq using the gradient of chisq
with respect to the ai) and the
inverse hessian method (using the curvature matrix to
compute the answer when we are close to the solution).
The value Lambda added to the diagonol elements of the
hessian matrix moves you between these 2 modes. When
lambda is large then the diagonal elements become dominant
and the motion matrix solution is just the gradient motion.
When lambda becomes small then the of diagonal terms become
important and you are solving the curvature matrix.
flambdastep determines how fast you move between these two
solution methods. When things get worse, lambda is increased by
a factor of flambdastep and you move towards the linear descent.
when things get better, then lambda decreases by flambdastep and
you move towards solving the curvature matrix.
note that in both instances the input data array is the
difference between the input data and the current fit value.
NOTATION;
let i=0-npts-1
let m=0-nparams
delta:i = y(i) - yfit(i) i=0,npts-1
pder:i,m= dByda:m at y(i)
alpha = transpose(pder) # ((Weights # (fltarr(nterms)+1))*pder)
alpha:m,n= sum:i( pder:m,i * ( 1/sig^2:i * pder:i,n)
so this is the hessian matrix with sig^2 included..
beta:m = sum:i(delta:i # pder:i,m) linear step for correction
sum deriv of each a:m over all y(i)i
alpha is the hessian matrix
b=beta is the current error in the fit.
EXAMPLE: Fit a function of the form f(x) = a * exp(b*x) + c to
sample pairs contained in x and y.
In this example, a=a(0), b=a(1) and c=a(2).
The partials are easily computed symbolicaly:
df/da = exp(b*x), df/db = a * x * exp(b*x), and df/dc = 1.0
Here is the user-written procedure to return F(x) and
the partials, given x:
pro gfunct, x, a, f, pder ; Function + partials
bx = exp(a(1) * x)
f= a(0) * bx + a(2) ;Evaluate the function
IF N_PARAMS() ge 4 THEN $ ;Return partials?
pder= [[bx], [a(0) * x * bx], [replicate(1.0, N_ELEMENTS(f))]]
end
x=findgen(10) ;Define indep dep variables.
y=[12.0, 11.0,10.2,9.4,8.7,8.1,7.5,6.9,6.5,6.1]
Weights=1.0/y ;Weights
a=[10.0,-0.1,2.0] ;Initial guess
yfit=curvefit(x,y,Weights,a,sigma,function_name='gfunct')
print, 'Function parameters: ', a
print, yfit
end
MODIFICATION HISTORY:
Written, DMS, RSI, September, 1982.
Does not iterate IF the first guess is good. DMS, Oct, 1990.
Added CALL_PROCEDURE to make the function's name a parameter.
(Nov 1990)
12/14/92 - modified to reflect the changes in the 1991
edition of Bevington (eq. II-27) (jiy-suggested by CreaSo)
Mark Rivers, U of Chicago, Feb. 12, 1995
- Added following keywords: ITMAX, ITER, TOL, CHI2, NODERIVATIVE
These make the routine much more generally useful.
- Removed Oct. 1990 modification so the routine does one iteration
even IF first guess is good. Required to get meaningful output
for errors.
- Added forward difference derivative calculations required for
NODERIVATIVE keyword.
- Fixed a bug: PDER was passed to user's procedure on first call,
but was not defined. Thus, user's procedure might not calculate
it, but the result was THEN used.
Steve Penton, RSI, June 1996.
- Changed SIGMAA to SIGMA to be consistant with other fitting
routines.
- Changed CHI2 to CHISQ to be consistant with other fitting
routines.
- Changed W to Weights to be consistant with other fitting
routines.
_ Updated docs regarding weighing.
(See /pkg/rsi/local/libao/phil/gen/curvefitpp.pro)
NAME:
daynotodm - convert daynumber to day,month
SYNTAX: [day,month]=daynotodm(daynum,year)
ARGS:
daynum: int/long daynumber of year 1..365or 366
year : int/long 4 digit year
RETURNS:
[day,month] as a vector.
DESCRIPTION
convert daynumber and year to day of month (1..31) and
month of year (1.l12).
(See /pkg/rsi/local/libao/phil/gen/daynotodm.pro)
NAME:
daynotojul - convert daynumber,year to julday
SYNTAX: julday=daynotojul(dayno,year)
ARGS:
dayno[n]: int/long/double daynumber of year 1..365or 366
year[n]: int/long 4 digit year
KEYWORDS:
gmtoffHr: double offset from gmt for dayno,year. gmtOffHr/24. will
be added to the computed julian days. Probably best
used when dayno id double or float.
RETURNS:
julday[n]:double julian day. This starts at noon
DESCRIPTION
convert daynumber and year to julian daynumber with fraction of day.
Method:
1. loop for each year.
a. take the first day of the year (day1), convert it to long (iday).
b. convert iday1 to julianday
c. for all the data of the year juldayYr=julday1 + day-iday1
(See /pkg/rsi/local/libao/phil/gen/daynotojul.pro)
NAME:
daysinmon - return number of days in this month
SYNTAX: days=daysinmon(mon,year)
ARGS:
mon : int/long month of year 1..12
year : int/long 2/4 digit year
RETURNS:
days : int/long number of days in this month
DESCRIPTION:
Return the number of days in the specified month.
If two digit year is entered then value gt 50 are 1999
(See /pkg/rsi/local/libao/phil/gen/daysinmon.pro)
NAME: dbit - convert to db's SYNTAX: a=dbit(b,minval=minval) ARGS: b[]:input value KEYWORDS: minval: float.. all values < this set to minval before log RETURNS: a[]: b in db's
(See /pkg/rsi/local/libao/phil/gen/dbit.pro)
NAME:
delpath - remove pathname from the path variable.
SYNTAX: delpath, path
ARGS:
path : string complete pathname to delete from the !path variable.
It must match how it appears in the !path variable.
(See /pkg/rsi/local/libao/phil/gen/delpath.pro)
NAME: dms1_deg - convert ddmmss.sss as a double to degrees. SYNTAX: angDeg=dms1_deg(ddmmss) ARGS: ddmmss: double value to convert RETURNS: angdeg: double the angle converted to degrees. DESCRIPTION Convert packed degrees, minutes, seconds to degrees. The input is a single double with ddmmss.ss with dd degrees, mm minutes, ss.s seconds.
(See /pkg/rsi/local/libao/phil/gen/dms1_deg.pro)
NAME: dms1_dms3 - convert deg,min,secs 1 word to deg,min,sec separate words SYNTAX - ret=dms1_dms3(hhmmss) ARGS: ddmmss : double angle to convert RETURNS: ret[4] : double deg,min,sec, and sign
(See /pkg/rsi/local/libao/phil/gen/dms1_dms3.pro)
NAME: dms1_rad - convert ddmmss.sss as a double to radians. SYNTAX: angRad=dms1_rad(ddmmss) ARGS: ddmmss: double value to convert RETURNS: angRad: double the angle converted to radians. DESCRIPTION Convert packed deg, minutes, seconds to radians. The input is a single double with ddmmss.ss with dd deg, mm minutes, ss.s seconds.
(See /pkg/rsi/local/libao/phil/gen/dms1_rad.pro)
NAME:
dmtodayno - convert day,mon,year to daynumber
SYNTAX: daynum=dmtodayno(day,mon,year)
ARGS:
day[] : int/long day of month
mon [] : int/long month of year 1..12
year[] : int/long 4 digit year
RETURNS:
daynum[]: int/long daynumber of year. First day of year is 1.
DESCRIPTION:
Convert from dayofmonth, month , and year to daynumber of year.
It also works with arrays.
(See /pkg/rsi/local/libao/phil/gen/dmtodayno.pro)
NAME:
dmytoyymmdd - convert ddMonyy to yymmdd
SYNTAX: yymmdd=ddmytoyymmdd(ddMONyy)
ARGS:
ddMONyy: string convert 10mar02 to 020210 etc..
RETURNS:
yymmdd: long return 0 if bad format..
DESCRIPTION
The datataking files use ddMONyy in the name where dd is the day
of the month, MON is a 3 letter abbreviation for the name,
and yy is the last two digits of the year. This routine will
convert the value into a long yymmdd
(See /pkg/rsi/local/libao/phil/gen/dmytoyymmdd.pro)
NAME:
dotprod - compute the dot product of two vectors
SYNTAX: val=dotprod(v1,v2)
ARGS:
v1[m] : vector
v2[m] : vector
returns:
val
DESCRIPTION:
return val=total(v1*v2)
(See /pkg/rsi/local/libao/phil/gen/dotprod.pro)
NAME:
explain - list documentation
SYNTAX: explain,{subject_document or routine_name}
ARGS:
None : If no arguments, display the list of subject documents.
namedoc:string If namedoc is provided then display the documentation
for this subject. It will be a list of routine names
with 1 line descriptions (eg cordoc) or a just a list
of the routine names (cordocnames).
routine :string If the name of a routine is entered, display the
complete documentation for the routine.
EXAMPLES:
explain - list the topics available
explain,cordoc - 1 line description of the correlator routines.
explain,cordocnames - list the names of all of the correlator routines.
explain,corplot - list the documentation for corplot.
(See /pkg/rsi/local/libao/phil/gen/explain.pro)
NAME:
fftint - integer based fft routine
SYNTAX: fftint,xr,xi,lenfft,bshift=bshift,doplot=doplot,coefBits=coefBits,$
cosAr=cosAr,sinAr=sinAr
ARGS:
xr[lenfft]: long real data
xi[lenfft]: long imaginary data
lenfft : long length of fft
KEYWORDS:
bshift: long bitmap telling whether or not to downshift on each butterfly
stage. For an fft of lenfft there are nbut=alog2(lenfft) butterfly
stages. Bit nbut-1 (counting from 0) is the first butterfly stage.
A one in the bit bitposition --> downshift, a 0 is no downshift.
The default value is 0x555555 downshift everyother butterfly.
coefBits:long number of bits to use in the coefficients. The default is 16
doplot: if set then plot each stage of the butterfly and wait for
users response. (s --> stop in routine, e--> exit routine, any
other key is continue).
cosAr[lenfft]:l64 input/output the cos table.
sinAr[lenfft]:l64 input/output the sin table.
RETURNS:
xr[lenfft] long the fft values are returned in place.
yr[lenfft] long
cosAr[lenfft]:l64 input/output the cos table.
sinAr[lenfft]:l64 input/output the sin table.
DESCRIPTION:
Perform and integer based fft. The user inputs the data as real and imaginary
arrays. Lenfft should be a power of 2. The coef's are scaled to numCoef bits (16 by
default). The multiplications are done in long long format to prevent overflows. The
return data is converted back to long.
If you supply cosAr,sinAr (with length ne lenfft) then the routine will return the
cosar, sinar used in the computation. On the next call you can input these arrays
and save the time to compute the sin,cos's.
You can use this routine to investigate integer round off errors, bit size on
different transformlengths.
Each butterfly step has been vectorized. For and 8k transform it runs about
9 times slower than the idl fft routine (50 vs 6 milliseconds).
EXAMPLE:
;
; make a sine wave with peak value of 35 counts
pksin=35
lenfft=1024
xr=long(mksin(lenfft,10)*pksin)
xi=long(mksin(lenfft,10,phase=-.25)*pksin)
fftinf,xr,xi,lenfft,bshift=0x2aa
freq=lindgen(lenfft) - lenfft/2
spcpwr=shift((xr^2L + xi^2L),lenfft/2L)
plot,freq,spcpwr
(See /pkg/rsi/local/libao/phil/gen/fftint.pro)
NAME:
fignum - put figure number on the page
SYNTAX: nextnum=fignum(fnum,xp=xp,ln=ln)
ARGS :
fnum : int figure number to put on plot. 1..
KEYWORDS:
xp : float 0..1 xposition for start of FIG N
ln : int 3..33 linenumber for FIG N
RETURNS:
nextnum : int input value incremented by 1
DESCRIPTION:
fignum() will place the string FIG fnum on the plot for you. The
horizontal position defaults to .92 of the screen (where the screen goes
0..1 horizontally). The vertical position is set to 3 where the vertical
screen runs 0 through 33.
NOTE:
the vertical line numbers are for the entire screen. If you use !p.multi
then you will have to decrease ln= by the corresponding amount.
SEE ALSO: note
(See /pkg/rsi/local/libao/phil/gen/fignum.pro)
NAME:
file_exists - check if a file name exists
SYNTAX: stat=file_exists(filename,fullname,dir=dir,size=size)
ARGS :
filename: string filename to search for
KEYWORDS:
dir[]: string if supplied then search through these directories
size : if supplied then return the file size in bytes
if it exists.
RETURNS:
stat : 1 file found, 0 not found
fullname: full directory/filename where file was found
size: if keyword size supplied then return the number of bytes in the
file
DESCRIPTION:
Check if a file exists. Return 1 if it does, 0 if it doesn't. Also return
the fully qualified name where the file was found.
If keyword dir is supplied then search through all of the directories
in the string array dir. In this case filename should not contain a
directory path.
This routine is handy to search for the location of an online datafile.
They start in /share/olcor/ but get moved to /proj/projid/ directories
at some later point.
EXAMPLE:
istat=file_exists('/share/olcor/corfile.13aug02.x101.1',fullname)
dir=['/share/olcor/','/proj/x101cor/']
istat=file_exists('corfile.13aug02.x101.1',fullname,dir=dir)
NOTE: this routine will only find regular files, a
directory name will return a non-existant file.
(See /pkg/rsi/local/libao/phil/gen/file_exists.pro)
NAME:
fisecmidhms3 - secs from midnite to hh:mm:ss
SYNTAX: label=fisecmidhms3(secsMidnite,hour,min,sec,float=float,$
nocolon=nocolon)
ARGS:
secsMidnite: long seconds from midnite to format.
KEYWORD:
float: if set then return secs at float
nocolon: if set then do not return colons
RETURNS:
hour: long hour of day.
min: long minute of hour.
sec: long sec of hour.
lab: string formatted string: hh:mm:ss
DESCRIPTION:
Convert seconds from midnight to hours, minutes, seconds and then
return a formatted string hh:mm:ss. The 2 digit numbers are 0 filled to the
left. If the input data is float/double and the float keyword is
not provided,the the data is truncated to long. The float keyword rounds
to 2 digits beyond the decimal point.
(See /pkg/rsi/local/libao/phil/gen/fisecmidhms3.pro)
NAME:
fitazza - fit function to azimuth and zenith angle
SYNTAX: fitazza,az,za,y,fitI,weights=weights,fittype=fittype,
yfit=yfit,covar=covar,variance=variance,sigmaA=sigmaA,
singular=singular,chisq=chisq
ARGS:
az[npts] : float azimuth in deg
za[npts] : float za in deg
y[npts] : float data to fit
KEYWORDS:
weights[npts] : weights to use with fit. default = 1.
fittype : int. 1..4 see below
RETURNS:
fitI : {azzafit} structure. return fitinfo here.
yfit[npts] : fit values evaluated at input az,za returne here
covar[m,m] : normalized covariance matrix. m=# of coef.
variance[m] : of the diagnonal elements (not normalized)
sigmaA[m] : sigmas of the coef (not normalized)
singular : int, number of singular coef found.
chisq : float. chisq
DESCRIPTION:
There are 4 or 10 coefficients in the fit. The functional form depends
on the value of fittype.
fittype: coef: 0-3
1: c0 + c1*za + c2(za-14)^2 + c3(za-14)^3 + azterms.. use za-14 when za gt 14
2: c0 + c1*(za-10) + c2(za-10)^2+c3(za-10)^3 + azterms. for all za
3: y=(za-10)/10 then
c0 +c1*(y)+ c2*(2*y*y-1.) + c3*(4*y^3-3*y) + azterms. for all za
4: c0 +c1*za +c2*(za-14)^2 +c3*(za-14)^3 no azterms.. use za-14 when za gt 14
5: c0 + c1*(za-10) + c2(za-10)^2+c3*(za-10)^3+
costerm have 1az,3az and sinza*3az.. no 2az term
6: c0 + c1*(za-10) + c2(za-10)^2+c3(za-10)^3 no az terms
coef: 4-9
az terms:
c4*cos(az) +c5*sin(az) + c6*cos(2az) + c7*sin(2az)
c8*cos(3az) +c9*sin(3az)
The fit info is returned in the fit structure {azzafit}. It contains:
numCoef: 10L , for fit
fittype: 1 , 1-def,2-about za10,3-chebyshev 3rd order
freq : 0. , Mhz
coef : dblarr(10),coef.
sigmaCoef : dblarr(10), sigmas on each coef.
covar : dblarr(10,10), covariance matrix
chisq : 0.D , of fit
sigma : 0.D , of fit - data
zaSet : 14. , za cutoff for higher order,or pivot
rfNum : 0 , rcv num
pol : ' ' , a , b , i stokes I
type : ' ' ,gain,sefd,tsys,etc..
title : ' ' , for any plots top
ytitle : ' ' , for any plots left
date : ' ' fit run
Since the fittype is returned in the structure, later routines (fitazzaeval
etc..) know which functional form to use.
You must call @geninit once before calling this routine to define the
{azzafit} structure.
SEE ALSO:
fitazzaeval, fitazzapr, fitazzaplres, fitazzaprcov
(See /pkg/rsi/local/libao/phil/gen/fitazza.pro)
NAME:
fitazzaeval - evaluate the fitazza fit at az,za positions.
SYNTAX: val=fitazzaeval(az,za,fitI,azonly=azonly,zaonly=zaonly)
ARGS:
az[n] : float azimuth positions to evaluate fit.
za[n] : float zenith angle positions to evaluate fit.
fitI :{fitazza} fit info returned from fitazza.
KEYWORDS:
azonly: if set then only evaluate the az terms of the fit.
zaonly: if set then only evaluate the za terms of the fit.
DESCRIPTION:
fitazza() will do a fit to data points as a function of azimuth and
zenith angle. This routine will allow you to evaluate that fit at
any az,za you want. The azonly, zaonly keywords limit the evaluation
to only the az or za terms of the fit.
SEE ALSO:
fitazza
(See /pkg/rsi/local/libao/phil/gen/fitazzaeval.pro)
NAME:
fitazzalog - write fitazza info in tabular form to a file
SYNTAX: fitazzalog,fitI,lun,calval,calval2,title=title
ARGS :
fitI: {azzafit} structure holding info returned from fitazza.
lun: int the lun open to the file for output.
calval: float the cal value that was used in the fit.
calval2: float 2nd cal value if fit was to stokes I
KEYWORDS:
title: string, title to write in file.
DESCRIPTION
Print out the fit values returned from fitazza in tabular form to a file.
You should open the file for write access (maybe append if you don't
want to overwrite the file).
SEE ALSO:
fitazza, fitazzapr, fitazzaeval
(See /pkg/rsi/local/libao/phil/gen/fitazzalog.pro)
NAME:
fitazzaplres - plot residuals from az,za fit
SYNTAX: fitazzaplres,az,za,y,fitI,tit=tit,key=key,sbc=sbc,sym=sym
ARGS:
az[npts] : float azimuth positions used for fit
za[npts] : float za positions used for fit
y[2,nsbc,npts]: float raw data used for fit pola,polB
fitI[2,nsbc] : {azzafit} returned from fitazza
KEYWORDS:
key : int key=1 residuals by sample
2 residuals by za
default both
tit : string title of plot
sbc : int sbc to plot 0-3. default all
sym[2] : int symbols for pola, b def.. none
DESCRIPTION:
Plot the residuals from the az,za fit. The routine assumes that
you have fit the polA and polB data for a number of subcorrelators.
It is the callers responsibility to move the fitI single structure into
the array of structures passed into this routine.
(See /pkg/rsi/local/libao/phil/gen/fitazzaplres.pro)
NAME:
fitazzapr - print/plot info on the az,za fit
SYNTAX: fitazzapr,fitI,over=over,ln=ln,sclln=sclln,nocoef=nocoef,tit=tit,$
nosigma=nosigma,noplot=noplot,plterms=plterms,$
plcomb=plcomb,xoff=xoff,_extra=e,azreq=azreq,zareq=zareq,
xp=xp
ARGS :
fitI: {azzafit} structure holding info returned from fitazza.
KEYWORDS:
over: int, if set then overplot
ln: int, line # to start printing fit values . def= 3
sclln: float scale spacing between lines. default=1.
nocoef: int, if set, then don't bother to print fit coef.
nosigma: if set then don't output sigma For coeff.
they aren't meaningfull unless you specified weights
noplot : if set then don't do the line plot, just print the values.
tit: string, title for plot
plterms: 0 - plot az terms
1 - plot za terms
2 - plot all terms versus az
3 - plot all terms versus za
plcomb: if set then plot the combined fit rather than overplotting
each term (only makes a difference for az terms).
xoff : float offset to add to default az or za before evaluating
the fit (overplot does not lay on top of previous plot).
_extra: keyword=value. passed to plot,oplot routines.
(eg.. psym=2 to just plot symbols)
azreq[] : float azimuth values (deg) to compute fit at.
zareq[] : float za values (deg) to compute fit at.
xp : float xposition (0 to 1) for text (def .05)
DESCRIPTION
Print out the fit values returned from fitazza. Output to the terminal
and plot to a plotfile.
The default output is to plot the individual azimuth terms versus
az. The program will by default evaluate the fit at fixed azimuth and
za values. You can use azreq zareq to evaluate it at different values
(if you specify plterms=2,3 and you want to use azreq,zareq then you
must provide both values).
EXAMPLE:
f(za): 9.16089 +( 0.77126)*za +(-0.06284)*(za-14)^2 +( 0.06741)(za-14)^3
f(1az): 0.54982*cos(1az) + (-1.85612)*sin(1az)
f(2az): 1.30097*cos(2az) + ( 0.69671)*sin(2az)
f(3az): 0.86003*cos(3az) + ( 1.33455)*sin(3az)
SigCof: za 0.25757 0.02781 0.08049 0.01504
az 0.18664 0.14984 0.17228 0.17435 0.16076 0.14014
f(za): ffffffff +(ffffffff)*za +(ffffffff)*(za-14)^2 +(ffffffff)(za-14)^3
f(1az): ffffffff*cos(1az) + (ffffffff)*sin(1az)
f(2az): ffffffff*cos(2az) + (ffffffff)*sin(2az)
f(3az): ffffffff*cos(3az) + (ffffffff)*sin(3az)
SigCof:az f8.5 f8.5 f8.5 f8.5
za f8.5 f8.5 f8.5 f8.5 f8.5 f8.5
(See /pkg/rsi/local/libao/phil/gen/fitazzapr.pro)
NAME:
fitazzaprcov - print out the covariance matrix
SYNTAX, fitazzaprcov,fiti,fd
ARGS:
fitI: {fitazza} returned from fitazza
fd: int write the data to this file descriptor. Default is stdout.
(See /pkg/rsi/local/libao/phil/gen/fitazzaprcov.pro)
NAME:
fitazzarob - robust fit of function to azimuth and zenith angle
SYNTAX: fitazzarob,az,za,y,fitI,weights=weights,fittype=fittype,
yfit=yfit,covar=covar,variance=variance,sigmaA=sigmaA,