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: alfamoninit - initialize to alfa dewar monitor routines SYNTAX: @corinit DESCRIPTION: call this routine before using any of the alfa deware monitor (amxxx) idl routines. It sets up the path for the idl alfamon directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/alfamoninit.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: angletosterad - angle to steradian SYNTAS: sterad=angletosterad(theta,deg=deg) ARGS: theta: float/double angle in radiabns to convert to steradians KEYWORDS: deg : if set then theta is in degrees. RETURNS: sterad: double steradians that angle subtends DESCRIPTION: Convert theta to steraradians. This assumes that theta is the full opening angle of a cone (not the half angle)l You can derive this by integrating the area of a sphere from 0 to th/2 with a radius of 1 for below let th be the half angle... a= integral(2*pi*rPrime*dthp) integrate thp from 0 to th but rPrime=sin(thp) if r =1 so.. a =2*pi *integral(sin(thp)*dthp) from 0 th = 2*pi*(-cos(thp)) eval at 0,th = 2*pi*(1-cos(th)) except this routine wants the full angle, not the half angle ...
(See /pkg/rsi/local/libao/phil/gen/angletosterad.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:
aocatinp - input an ao source catalog
SYNTAX: nsrc=aocatinp(file,srcI,calib=calib,galhI=galhI,gcl=gcl,lewis1=lewis1,$
lewis2=lewis2,oh=oh,psr=psr,brightpsr=brightpsr,snr=snr,$; sort=sort,novel=novel)
ARGS:
file :string filename of catatlog
KEYWORDS:
The following keywords select one of the standard catalogs available from
Online at AO. These keywords will not work if you are not running at AO
(since the files are not available).
calib : if set then use calib.cat source calibrator file
galhI : if set then use galaxy_HI.cat
Arecibo HI standard galaxy list
Objects with: no neighbor w/i 10' & 1000km/s radii;
Optical Diameter < 2.0'
gcl : if set then use gcl.cat
Globular Clusters in Harris' catalog visible from AO
lewis1 : if set then use lewis_gals.cat
Galaxy catalog (WITH velocities)
Lewis et al. (1985 ApJS 59 161)
lewis2 : if set then use lewis_gals2.cat
Galaxy catalog (WITHOUT velocities)
Lewis et al. (1985 ApJS 59 161)
oh : if set then use oh.cat
Bright OH/IR stars in the Arecibo Sky
Catalog for 1612 MHz velocity calibration sources
psr : if set then use the princeton pulsar catalog (our version
of it used by cima
brightpsr : if set then use bright pulsars from the princeton pulsar
catalog.
snr : if set then use snr.cat
Green's supernova remnant list
sort : if set then sort the sources by increasing ra
novel : if set then the catalog has no velocity info.
stop after the coordsys
RETURNS:
nsrc : int number of sources found
-1 could not read file.
srci[nsrc]:{srccat} return data here
DESCRIPTION:
Read in all of the sources in an ao catalog file. These are the standard
and usr source list files that are used by CIMA (the telescope gui control
program).
The file format is:
# col 1 is a comment
0 1 2 3 4 5 6
srcNm ra dec CoordSys vel velFrame velType # comments
ra : hh:mm:ss.s or hhmmss.s
dec: dd:mm:ss.s or ddmmss.s
Coordsys: J,B
vel : velocity of source relcative to vel coordinate system
velFrame: Topo - topocentric
Helio - helio centric
lsr - local standard of reset
velType : v - velocity km/sec
zo - z optical
zr - z radio
The returned srcI array will contain:
help,srcI,/st
history:
str:
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)
Coord string 'b','j'
vel double 0.
velCrdsys string '' t-topo,h-helio,l-lsr,
velType string '' v,zo,zr
EOL STRING '' string after : following velType.
(See /pkg/rsi/local/libao/phil/gen/aocatinp.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,nmOff=nmOff,tblln=tblln,xpsrc=xpsrc
cs1=cs1,csn=csn,ltit=ltit,lticklen=lticklen,p12m=p12m
decl=decl,fixedlen=fixedlen,font=font
ARGS:
az[npts]: float azimuth deg
za[npts]: float za deg
azerr[npts]: float azerr asecs(def)
zaerr[npts]: float zaerr asecs(def)
tit: string. title.. eg 'pointing errors'
KEYWORDS:
ticklen : float ticklen in asecs. def: 5
lticklen : str label for tick len. def=asec
hard : if set then use hardcopy settings
names[nsrc] : string list of source names used. print left of image.
indlist[npts]: long index into names array for each point in az,za
notab : if set then donot print avg/rms table by za
tabtit : title for err table by za. def:errors
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
tbllen : long if supplied then line number to start table
xpsrc : float xpostion to start srcnames. [0.,1]
default -.1 screen,-.2 plotting
ltit : int line to start the plot title
p12m : if set then for 12meter. rise, set are reversed
decl : if set then overplot dec lines (for 305m)
every 5 degrees
fixedlen : if set then the input data is fixed length, we are just using
the rotation angle (eg.fit residuals can to go zero..)
This will leave off the ticklen message.
over : if set then overplot. should be a 2nd call to routine
csn : float charsize for notes.
font : font to use 1-> truetype
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:
avgchanbyrms - use rmsbychan to throw out outliers when averaging
SYNTAX: avgAr=avgchanbyrms(d,median=median,rmsar=rmsar,nsig=nsig,cntar=cntar,$
nodiv=nodiv,double=double)
ARGS:
d[m,n] : array to average over 2nd dimension
KEYWORDS:
median: if set then use median rather than mean when doing the avg , and dividing
the rms by the mean.
rmsAr[m]: rms values for each chan used for thresholding
nsig : when averaging, use all points <- nsig. def=3
double : if set the compute using doubles
RETURNS:
avgAr[m]: average array
cntar[m] : long number averaged in each channel
rmsAr[m] : long number averaged in each channel
DESCTRIPTION:
compute the standard deviation/mean by channel.
(See /pkg/rsi/local/libao/phil/gen/avgchanbyrms.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:
basketweave - basket weave two maps
SYNTAX: basketweave,map1,map2,map1w,map2w,deg=deg,fract=fract
ARGS:
map1[m,n]: float/double input map 1. it was driven in the x coordinate
map2[m,n]: float/double input map 2. it was driven in the y coordinate
KEYWORDS:
deg: int if not zero then fit a polynomial of deg deg to the subtracted data
if zero, use a robust mean.
def: robust mean
fract: float fraction of fit to remove on each iteration, (0> to 1.).
verb : if set then plot results as we go along
RETURNS:
map1w[m,n]: float/double backet weave map1 using map1 - map2
map2w[m,n] float/double basket weave map2 using map2 - map1
DESCRIPTION:
Basketweave two maps. This will help remove stripes in a map.
When making the map you can:
- drive in ra and step in dec
or- drive in da and step in ra.
Stripes typically between the strips.
To basketweave to maps:
1. the maps should already be gridded to the same coord grid
2. copy map1 to map1w, copy map2 to map2w
3. compute the map difference mdif12= map1 - map2
4. for each xstrip
- fit a robust polynomial to the mdif12[*,i] (or just use the robust mean)
- remove a fraction of this difference from the map1w[*,i] (def=.7)
5. do the same to map 2
mdif21=map2 - map1
for each y strip
- fit a robust polynomial of deg to mdif21[i,*] (or just use the robust mean)
- remove a fraction of the fit from map2w[i,*]-=fract*yfit
Return map1w,map2w ...
Normal processing would iterate:
mapc1=map1
mapc2=map2
loop n
rms1=rms(mapc1-mapc2,/quiet) ; rms of difference before
basksketweave(mapc1,mapc2,map1w,map2w,deg=1)
rms2=rms(map1w-map2w,/quiet) ; rms of difference before
print.. rms before,afte
mapc1=map1w
mapc2=map2wk
endloop
you could then averaged the 2 basket weaved maps
(See /pkg/rsi/local/libao/phil/gen/basketweave.pro)
NAME: bdwfinit - initialize to use the idl bdwf (brown dwarf routines). SYNTAX: @bdwfinit DESCRIPTION: call this routine before using any of the bdwf_ idl routines. It sets up the path for the idl bdwf directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/bdwfinit.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,xfit=xfit)
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.. you should use xfit=xfit
as the xvalues for the fit;
yy=poly(xfit,coef)
mask[npts]: mask used for fit that was defined by the user.
xfit[npts]: x values used for fitting [0,1) or [0,-1] . these are
the values you should used when evaluating coef
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: cadiinit - initialize the cadi ionosonde software SYNTAX: @cadiinit DESCRIPTION: call this routine before using any of the cadi_... idl routines.
(See /pkg/rsi/local/libao/phil/gen/cadiinit.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, 1 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,p12m=p12m)
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.
alfaBmNum: 0..6. If alfa receiver then only return cal values
for this beam number. The default is all 6 beams.
p12m if set then this is 12meter data..
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, 1 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:
calget1fit - return linear fit to cal values
SYNTAX: stat=calget1fit(rcvrNum,caltype,f1,f2,coefAr,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
f1: float Mhz min freq
f2: float Mhz max freq
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.
alfaBmNum: 0..6. If alfa receiver then only return cal values
for this beam number. The default is all 6 beams.
RETURNS:
coefAr[2,npol=2]: float .. linear fit to polA,polB
Note if this is alfa, then coefAr=[2,npol,7]
value are returned for the 7 pixels unless
alfabmnum is selected.
stat: int .. -1 error, 1 got the values ok.
DESCRIPTION:
Return a linear fit to the calvalues between f1 and f2 for the requested
calType and receiver. The cal value in kelvins could then be computed as:
CalAK[f]=coefAr[0,0] + coefAr[1,0]*FreqMhz
CalBK[f]=coefAr[0,1] + coefAr[1,1]*FreqMhz
If the receiver is alfa then the returned coefs are:
coefAr[2,npol,nbeams] unless alfaBmNum keyword is used. In that case
the coef for the single beam are returned.
EXAMPLES:
Get the cal values for lbw (rcvrNum=5) using the high correlated cal
(caltype=1) between 1300 and 1500 Mhz.
stat=calget1fit(5,1,1300,1500,coefAr)
calAK=poly(findgen(200)+1300),coefAr[*,0])
NOTE:
SEE ALSO:calget, calval, calinpdata, corhcalval
(See /pkg/rsi/local/libao/phil/gen/calget1fit.pro)
NAME:
calInpData - input cal data for rcvr/calType.
SYNTAX:
istat=calinpdata(rcvNum,calNum,calData,fname=fname,date=date,p12m=p12m)
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.
p12m : if set then get 12meter cals. 11=xb,7=sb
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
3 610 0 low cor cal
6 lbn 0 low cor cal
100 ch 5 h uncorcal
sbn now has hcorcal,lcorcal.
12 sbn 0 low cor cal,hcorcal
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,$
alfaBmNum=alfaBmNum)
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.
alfaBmNum: 0..6 If specified then only return alfa data for this
beam. By default data for all 7 beams are returned.
If not alfa data, this is ignored
RETURNS:
istat: 1 ok within range, 0 outside range used edge,-1 all zeros,
-2 error
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:
calvalfit - fit calvalues to a freq range
SYNTAX:
istat=calvalfit(f1,f2,npnts,calData,coefAr,fitOrder=fitOrder,$
,hybrid=hybrid,swappol=swappol,alfaBmNum=alfaBmNum)
ARGS:
f1: float Mhz. first frequency
f2: float Mhz. last frequency
npnts: long number of points to return with f1,f2 being first
last freq values
calData{}: calinfo input via calinpdata
KEYWORDS:
fitOrder: int order for fit. Default is linear
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.
alfaBmNum: 0..6 If specified then only return alfa data for this
beam. By default data for all 7 beams are returned.
If not alfa data, this is ignored
RETURNS:
istat: 1 ok within range, 0 outside range used edge,-1 all zeros,
-2 error
coefAr[ncoef,npol]: float linear fit coef. c0 + c1*freq
If this is alfa then the calV is [ncoef,2,7] unless
alfaBmNum is supplied
DESCRIPTION:
Compute a linear fit to the cal data between f1 and f2
The calData should have already been read in with calInpData.
The data is returned in an array of (ncoef,2)
the fit is calVal= calval[0,*] + calval[1,*]*freqMhz
It is ok for f1 to be gt f2..
The normal way to get cal values is via mascalval()
They call this routine.
SEE ALSO:corhcalval, calget, calinpdata.
(See /pkg/rsi/local/libao/phil/gen/calvalfit.pro)
NAME:
cataloginp - input a pointing catalog
SYNTAX: nsrc=cataloginp(file,format,retdata,comment=comment,crdsys=crdsys)
ARGS:
file :string filename of catatlog
format :int format for catalog:
1: srcname hh mm ss dd mm ss
2: srcname hhmmss ddmmss
3: srcname hh:mm:ss dd:mm:ss
4: fmt 2 or 3. program checks for :
retdata[]:{srccat} return data here
in retdata.
KEYWORDS:
comment : string comment characters for catatlog.def:#
crdsys : string coord system for coordinates. by defaul
assume 'j': j2000. other is 'b': b1950
DESCRIPTION:
Read in all of the source names and positions catalog specified by
file
The returned srccat array will contain:
help,retdat,/st
history:
** 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: cgeninit - initialize to use cummings generator routines SYNTAX: @cgeninit DESCRIPTION: call this routine before using any of the cgen routines to access the cummings generator datafiles.
(See /pkg/rsi/local/libao/phil/gen/cgeninit.pro)
NAME: cginit - add path for coyote graphics SYNTAX: @cginit DESCRIPTION: call this routine before using any of the coyote graphics routines
(See /pkg/rsi/local/libao/phil/gen/cginit.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
pjp dbg : if set then print error messages when fit doesn't work
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,scl=scl)
ARGS: b[]:input value
KEYWORDS:
minval: float.. all values < this set to minval before log
scl : float.. rather than scale to max value,scale to scl..
(for no scaling set scl=1.)
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:
digmixtmd - digitally mix a complex timedomain waveform
SYNTAX: bufOut=digmixtmd(mixI,bufInp)
ARGS:
mixI : {} mixing info. Returned by digmixtmdinit()
bufInp[n]: complex complex data input (before mixing).
RETURNS:
bufOut[n]: complex data after mixing
DESCRIPTION:
digmixtmd() will mix complex time domain baseband data. The user first calls
digmixtmdinit() to define the sampling frequency and the new center frequency
after mixing.
You then call digmixtmd() as many times as needed with contiguous time domain
buffers. The lophase is stored in mixI and is updated on each call to digmixtmd().
Example:
Assume:
1. complex time domain data sampled at 100 Mhz
2. you want to mix 30 Mhz to the center of the band
3. Do thie for an entire files worth of data.
- open file
- mixI=digmixtmdinit(100e6,30e6)
bout=complexarr(maxsamples)
icur=0
- while ((nsamples=readbuf(bufInp)) != eof)
dout[icur:icur+nsamples-1]=digmixtmd(bufInp,mixI)
icur+=nsamples
endwhile
SEE ALSO:digmixtmdinit()
(See /pkg/rsi/local/libao/phil/gen/digmixtmd.pro)
NAME:
digmixtmdinit - initialize structure for digital mixing
SYNTAX: mixI=digmixtmdinit(smpFrq,newCfr)
ARGS:
smpFrq: double sample frequency of complex data.
newCfr: double new Center frequency for baseband data.
RETURNS:
mixI: {} structure holding mixing info
DESCRIPTION:
digmixtmdinit() initializes the mixI structure for calls to digmixtmd.
The data to be mixed is assumed to be complex baseband data (centered at 0 freq).
The user specifies smpFrq (in whatever units you want). This will define the time
spacing between individual complex samples (1/smpFrq).
The newCfr is the new center freq after mixing (in the same units as smpFrq).
It should be between +/- smpFrq/2. After mixing, this will be the new center of the
band.
The mixI structure contains:
** Structure <8c1ca8>, 4 tags, length=32, data length=32, refs=1:
SMPFRQ DOUBLE 1.5000000e+08
NEWCFR DOUBLE 68000000.
LOFRQ DOUBLE -68000000.
LOPHASE DOUBLE -0.52941176
Each time digmixtmd is called, the lophase is updated to the correct phase for
the next call to digmixtmd. This lets you make multiple calls to digmixtmd.
Example:
Assume:
1. data complex sampled at 100 Mhz
2. you want to mix 30 Mhz to the center of the band
3. you want to do this for an entire files worth of data.
- open file
- mixI=digmixtmdinit(100e6,30e6)
bout=complexarr(maxsamples)
icur=0
- while ((nsamples=readbuf(bufInp)) != eof)
dout[icur:icur+nsamples-1]=digmixtmd(bufInp,mixI)
icur+=nsamples
endwhile
(See /pkg/rsi/local/libao/phil/gen/digmixtmdinit.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
KEYWORDS:
YYMMDD: if supplied then use this rather than day,mon,year
SYNTAX: daynum=dmtodayno(day,mon,year,yymmdd=yymmdd)
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:
dopcorbuf - init for dopcorbuf()
SYNTAX:dopCorI=dopcorbufinit(dopFrqHz,dopTmU,bufStTmU,smptmStepU,$
tmUToSec,bandflipped=bandflipped)
ARGS:
dopFrqHz[n]: double measured doppler frequencies
dopTmU[n] : double times for each dop Frq
bufSttmU : double start time for first sample of buf
smpTmStepU : double time step for each sample in buf
tmUToSec : double value that converts user time units to
seconds.
KEYWORDS:
bandflipped: if true then freq band had been flipped
this will flip the sign of the
doppler used for the correction
RETURNS:
dopcorI :{} structure to pass to each call of dopcorbuf()
DESCRIPTION:
Initialize for calls to dopcorbuf().
The user passes an array of doppler frequencies as
well as the timestamp of each of these. The frequency is
positive if the signal frequency is greater than the rest frequency
(blue shifted.. object coming toward us)
The user also provides the timestamp for the first
sample in the buffer as well as the timestep between
samples in a buffer.
The time value can be anything the user want, but
it must be the same for all of the provided times.
All the variables ending in U are in these user time units.
The variable tmUToSec converts from the users time units
to seconds (eg. if units are hours, then tmUToSec=3600D).
The user must also verify there is enough resolution
in the data type (best to use double).
Set the bandflipped keyword if the sampled data band has been
flipped in frequency (by an odd number of high side mixing stages).
After calling this routine, the user should
make multiple calls to dopcorbuf() starting with the
buffer that starts at bufStTm.
Notes:
- this routine does not correct for time dilation.
- resolution:
using hours and double, the time resolution is 3.6*10-13 sec
using mjd and double, the time resolution is .86 usecs
using jd and double, the time resolution is .86 millisecs
For sample times on the order of usecs, it's probably best to
use hours for the time unit.
- these two routines probably need some more debugging.
EXAMPLE:
Suppose the doppler values have units of 1day, and the complex
signal has bandwidth BW (in hz) then
dopFrq[n] - doppler frequencies
dopTmU[n] - these will have units of 1 day.
bufStTmU - starting daynumber.fract for first sample in buffer
in dayno units
smpTmUStep = (1/bandwidhz)/(86400D) . this is the sample time in
units of 1day.
tmUToSec = 86400D .. to convert user time unit (1day) to seconds
(See /pkg/rsi/local/libao/phil/gen/dopcorbufinit.pro)
NAME:
dopcorbuf - offline doppler correct voltage data
SYNTAX:bufCor=dopcorbuf(dopI,buf,firstBuf=firstBuf)
ARGS:
dopI: {} structure holding dopper info. initialize
with dopcoroflinit(0
buf[n]: complexar data to doppler correct
KEYWORDS:
firstBuf: if true then assume this is the first buf
It will reset the buf time to the time of
first buf, and set the phase offset to 0.
remdc : if true, then remove dc from voltages before
doppler shifting
RETURNS:
bufCor[n]: complexarr doppler corrected data
dopI.dopFrq[m]: doppler frequency for m tims
dopI.dopTms[m]: time for each doppler value
dopI.tmStpSmp : time step each sample
these get updated on each call
dopI.phoff : phase off for start of buf
dopI.tmStBuf : updated on each call
DESCRIPTION:
Offline doppler correct voltage data. You must
call dopcorbufinit() before starting to call dopcorbuf().
The init routine will store the doppler info, startime
of 1st buffer, sampletime, and current phase offset.
Each call to this routine will doppler correct the
buffer and increment the start time and phase of the next
buffer.
The sequence should be:
dopI=dopcorbufinit()
loop from first buff
input next buffer
bufC=dopcorbuf(dopI,buf)
endloop
If you want to start over, you could recall dopcorbufinit()
or call dopcorbuf(...,/firstbuf). /firstbuf will reset the
current buf time to the time of the first buffer, and 0 the
accumulated phase.
(See /pkg/rsi/local/libao/phil/gen/dopcorbuf.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:
drawcircle - draw a circle
SYNTAX: drawcircle,x,y,radius,npnts=npnts,color=color,thick=thick,$
arc=arc,xret=xret,yret=yret,$
_extra=_e
x,y float center of circle
radius: float radius of circle
KEYWORDS:
npnts: int number of points to use in drawing the circle. The default
is 9 points. The circle will be drawn as an n-1 sided
polyhedron.
color: int number 0 to 10. Color to use when plotting symbol. The
default is the same color as the lines beging drawn. See
shcolsym for a mapping of numbers to colors.
thick:float the line thicknes to use when drawing the lines. 1. is the
default.
arc[2]:fltarr [deg] if supplied then draw an arc rather than a full circle
arc[0]= min theta
arc[1]= max theta.. theta is measured CCW from the x axis
xret[npnts]:fltarr return the x points of the circle/arc
yret[npnts]:fltarr return the y points of the circle/arc
over : If set then oplot rather than plot
_extra=_e any other params to plot,oplot
DESCRIPTION:
draw a circle about x,y with given radius
(See /pkg/rsi/local/libao/phil/gen/drawcircle.pro)
NAME: ecalinit - initialize to use the idl e callisto routines SYNTAX: @ecalinit DESCRIPTION: call this routine before using any of the e callisto idl routines. It sets up the path for the idl ecal directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/ecalinit.pro)
NAME:
ephmaord - read in an ao ephmeris file.
SYNTAX: n=ephmaord(fname,ephmI,rangedop=rangedop,$
yyyymmdd1=yyyymmdd1,yyyymmdd2=yyyymmdd2)
ARGS:
file: string filename
KEYWORDS:
rangedop: if true then range, doppler ar last two fields
yyyymmdd1: double first record to return is >= to this
utc date. includes fraction of day.
yyyymmdd2: double last record returned is <= this utc date
e. includes fraction of day.
RETURNS:
n : int -1 if can't read file
>0 number of entries
ephmI: {} structure holding info.
DESCRIPTION:
input an ao ephemeris file. This is the format
used by the pointing program. The format is:
UTC RA DEC
yyyy-mm-dd hh:mm:ss hhmmss.ss +ddmmss.s az el
the ra,dec can be j2000 or of date (current).
The data is returned in the structure
help,ephmI,/st
TARGET STRING 'Moon (301)'
J2000 BYTE 0
N LONG 212
DAT STRUCT -> Array[212]
help,ephmI.dat,/st
JD DOUBLE 2456503.9
YYYYMMDD LONG 20130730 utc
SECMID FLOAT 32460.0 utc
RAHR DOUBLE 2.8377028
DECD DOUBLE 16.014861
AZD DOUBLE 91.399200 .. Source azimuth
ZAD DOUBLE 24.834300
if keyword set /rangedop add:
range double AU
dop double km/sec .. neg is approaching
See /share/megs/phil/x101/ephm/ obs_tbl.sc and hortoao.pl to grab
horizons file and convert to ao format
(See /pkg/rsi/local/libao/phil/gen/ephmaord.pro)
NAME:
ephmaotordr - convert ao short to ao radar ephemeris
SYNTAX: istat=ephmaotordr(aofileIn, rdrfileOut)
ARGS:
aofileIn: string input file to convert
rdrfileOut: string output file for radar format
DESCRIPTION:
The ao pointing program supports tracking using an ephemeris file
(more info:http://www.naic.edu/~phil/software/datataking/vxworks/pntCmd.html#ephem).
There are two format supported:
- standard ao radar ephemeris (rcv and tx positions)
- ao short format ( just the date time ra dec az el range rangeRate
You can download ao short ephemerides from jpl horizons:
1. see /share/megs/x101/ephm/obs_tbl.sc (or one of the scripts that calls it
2. hortoao.pl < jplInp > aoOutput .. will format the jpl output to ao short format
This routine will then take the ao short format and generate a radar ephemeris file.
To do this:
1. the ao format file should have used -q 2,4,20 options to obs)tbl.sc ..
This returns apparent ra,dec,az,el and well had range and range rate
2. I reads the ao format file. copying the comments (#) to the output file
3. For each data line it generates the rcv and tx data (the tx data is a copy of the rcv data).
- the delay is computed as rangeAu*AutoKm/vellight
20may20.. update to have a tx set that is close to correct
(except that the horizon's ephm using venus is not an exact match of the actual
radar version).
SEE ALSO
an example of using this routine to get radar ephm from tle's can be found at
/share/obs4/usr/radar/r3399/mkephm.sc
(See /pkg/rsi/local/libao/phil/gen/ephmaotordr.pro)
NAME:
ephmradarrd - read radar ephemeris file
SYNTAX: n=ephmradarrd(fname,ephmI)
ARGS:
fname: string name of ephemeris file to input
RETURNS:
n : > 0 number of entries read
< 0 error code
ephmI[n]: struct holding ephm info
DESCRIPTION:
Read in an ao radar ephmfile. This is the format that is passed to the
datataking software. It is jon's format after being massaged be ??? script.
The return info is:
help,ephmI,/st
JD DOUBLE 2456505.2 juldate of entry
SECS LONG 1375289760 secs1970 of entry
STEPSECS DOUBLE 60.000000 step (secs) of entries in ephmfile
YMD LONG 20130731 yyyymmdd of entry (utc)
HMS LONG 165600 hhmmss of entry (utc)
LYMD STRING '2013-07-31' yyyy-mm-dd as string (utc)
LHMS STRING '16:56:00' hh:mm:ss as string (utc)
RXI STRUCT -> Array[1] receive info
TXI STRUCT -> Array[1] tx info
help
IDL> help,ephmI.rxI,/st
EL FLOAT 40.4000 elevation (deg) to .1 deg
AZS FLOAT 59.0000 source az (deg) to .1 deg
AZF FLOAT -121.000 feed az (deg) to .1 deg
RAHMS DOUBLE 123858.80 apparent ra hhmmss.sss
DECDMS DOUBLE 351243.50 apparent dec ddmmss.ddd
RAH DOUBLE 12.649667 apparent ra hh.hhh
DECD DOUBLE 35.212083 apparent dec dd.ddd
DELAY DOUBLE 27.067925 delay secs (rtt)
ephmI.txI is similar for the tx portion
(See /pkg/rsi/local/libao/phil/gen/ephmradarrd.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: fftinterp - fft interpolation of real data SYNTAX: fftinterp,n,yin,yout, ARGS: n : interpolation factor yin[m] : float input data yout[m*n] : flout interpolated data DESCRIPTION: interpolate the yin function by a factor of n using the fft. This probably needs a better filter. You get lots of ringing between interpolated points.
(See /pkg/rsi/local/libao/phil/gen/fftinterp.pro)
NAME:
fignum - put figure number on the page
SYNTAX: nextnum=fignum(fnum,xp=xp,ln=ln,chars=cs)
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
chars : float size of char. def=1
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:
filereadsock - read file via socket.
SYNTAX: istat=filereadsock(cpu,fname,dat,verb=verb ,tm=tm,usesh=usesh)
ARGS:
cpu: string name of remote cpu where file is located
fname: string name of the file to read (on remote machine).
dat : xxx data to read into. User should allocate prior to
call. Data will be read to fill this array.
KEYWORDS:
verb: int if set then output some debugging info (dd cmd etc..)
tm: double if supplied the return time to execute routine
usesh: int if true then use a shell. If not supplied or set to
0 then start child without shell.
returns: istat: 1 ok
-1 error probably not enough data
see !error_state.msg
DESCRIPTION:
We normally use nfs to read remote files. This routine will read
a disc file using a socket and rsh. The command used is:
- /usr/bin/rsh cpu /bin/dd bs=65536 if=fname conv=notrunc status=none
How it works:
- The idl spawn routine is used with the lun= option
- idl spawns the command and then attaches an lun to the
output of the pipe.
- The routine then does a readu,dat
to read in the data.
- when done the pipe is closed.
What you need to use this routine:
- rsh must be enabled on the machine you are running on.
- the remote cpu must be in your .rhosts file
Data will be read to fill the dat[] array the user passes in (so you must
allocate it prior to calling this routine.
- If the read hits eof before filling dat[] then istat returns with -1
and !error_state.msg will have the error message.
(See /pkg/rsi/local/libao/phil/gen/filereadsock.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/float/double 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:
fitazeleval - evaluate the 12m fit at the az,el positions.
SYNTAX: val=fitazeleval(az,el,fitI)
ARGS:
az[n] : float az positions to evaluate fit. (if fittype does not use az, just
enter a dummy value)
el[n] : float el positions to evaluate fit.
fitI :{azelfit12m} fit info to use for fitting
DESCRIPTION:
fitazeleval() will evaluate the fit stored in fitI. Currenly this
is just the 12m gain curves. eventually i'll write fitazel() to make it more general
This routine will allow you to evaluate that fit at
any az,el you want.
SEE ALSO:
(See /pkg/rsi/local/libao/phil/gen/fitazeleval.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
7: c0 + c1*cos(az) +c2*sin(az) + c3*cos(2az) + c4*sin(2az)
c5*cos(3az) +c6*sin(3az)
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.54/982*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,
singular=singular,chisq=chisq,
gindx=gindx,ngood=ngood,bindx=bindx,nbad=nbad,$
maxiter=maxiter,fpnts=fpnts,iter=iter,nsig=nsig
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
nsig : float when iterating throw out points with
residuals gt nsig*sig.
zapivot : float za for pivot. The default is 14. or 10. depending
on the fittype.
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
ROBUST INFO:
fpnts : float the fraction of points used for the final
computation
gindx: long[] indices into d for the points that were used
for the computation.
ngood long number of points in gindx.
bindx: long[] indices into d for the points that were not used
nbad long number of points in bindx.
iter long number of iterations performed.
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
7: c0 + c1*cos(az) +c2*sin(az) + c3*cos(2az) + c4*sin(2az)
c5*cos(3az) +c6*sin(3az)
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/fitazzarob.pro)
NAME:
fitngauss - fit n gaussians
SYNTAX: coef=fitngauss(x,y,ngauss,coefInit,chisq=chisq,sigmaCoef=sigmaCoef,
weights=weights,yfit=yfit,covar=covar,trouble=trouble)
ARGS:
x[npts]: inpependent data
y[npts]: dependent data
ngauss : int the number of gaussians to fit for 1 to 3
coefInit[m] initial values for coef.
KEYWORDS:
weights[npts] : for data. default is unity
trouble : int 0 converged ok
-1 chisq infinite
-2 flambdacount > max=30 *10/flambdaste
-3 iterdone >itmax(def 20)
-4 divide by zero on partials
_extra=e
cfplot : if set, then curveFit will plot each strip
cfparms : if set, then curveFit will print input parms
flambdastep : float. how much we decrease step curvature matrix when
tol : float. when chisq decreases by this amount, done..def:.001
itmax : int.. iterations max. default:40
RETURNS :
coef[m] fit coef..
keywords:
yfit[npts]: fit evaluated at input values
chisq : double (y-yfit)^2/(npts-mcoef)
sigmaCoef[m]:formal errors in coef. sqrt(diag of matrix*chisq)
covar[m,m] : covariance matrix
DESCRIPTION:
fit the function:
y= a0 +a1*x + a3*exp[(x-a4)^2/a5*(1+alpha*x) +
a6*exp((x-a7)^2/a8 +
a9*exp((x-a10)^2/a11
a0= constant
a1= linear in x
a2= alpha skew first gaussian
a3=Ampl
a4=Mean (offset error)
a5=sigma input as fwhm in units of x
a6=Ampl
a7=Mean
a8=sigma input as fwhm
a9=Ampl
a10=Mean
a11=sigma input as fwhm
--> 19mar19 .. note: coma has been disabled...
(See /pkg/rsi/local/libao/phil/gen/fitngauss.pro)
NAME:
fitngaussfunc - function for fitting n gaussians
SYNTAX:
fitngaussfunc,x,a,f,pder
ARGS:
x[npts] : independant variable
a[2+4*ngauss] : parameters to fit for
f[npts] : return value of function here
pder[2+4*ngauss] :return partial derivatives with respect to each param
DESCRIPTION:
evaluate the function and optionally it's partial derivatives
for the curvefit routine of idl:
f(x)= a0 + a1*x + a3*exp[-(x-a4)**2 * (1+a2*(x-a4))]
-------------------- + n more gaussians
a5^2
a0= constant
a1= linear in x
a2=skew first gaussian
a3=Ampl
a4=Mean
a5=sigma
a6=Ampl
a7=Mean
a8=sigma
a9=Ampl
a10=Mean
a11=sigma
note that we are fitting for sigx not sigx^2
(See /pkg/rsi/local/libao/phil/gen/fitngauss.pro)
NAME:
fitngaussnc - fit n gaussians (no coma)
SYNTAX: coef=fitngaussnc(x,y,ngauss,coefInit,chisq=chisq,sigmaCoef=sigmaCoef,
weights=weights,yfit=yfit,covar=covar,trouble=trouble,
ARGS:
x[npts]: inpependent data
y[npts]: dependent data
ngauss : int the number of gaussians to fit for 1 to 3
coefInit[m] initial values for coef.
KEYWORDS:
weights[npts] : for data. default is unity
trouble : int 0 converged ok
-1 chisq infinite
-2 flambdacount > max=30 *10/flambdaste
-3 iterdone >itmax(def 20)
-4 divide by zero on partials
_extra=e
cfplot : if set, then curveFit will plot each strip
cfparms : if set, then curveFit will print input parms
flambdastep : float. how much we decrease step curvature matrix when
tol : float. when chisq decreases by this amount, done..def:.001
itmax : int.. iterations max. default:40
iter : int.. iterations done.
noderiv : if set then have let the routine compute the derivatives
from differences rather than from the analytic formula.
RETURNS :
coef[m] fit coef..
keywords:
yfit[npts]: fit evaluated at input values
chisq : double (y-yfit)^2/(npts-mcoef)
sigmaCoef[m]:formal errors in coef. sqrt(diag of matrix*chisq)
covar[m,m] : covariance matrix
DESCRIPTION:
fit the function:
y= a0 +a1*x + a2*exp[(x-a3)^2/a4 +
a6*exp((x-a7)^2/a8 +
a9*exp((x-a10)^2/a11
where:
a0= constant
a1= linear in x
a2=Ampl 1st gaussian
a3=Mean (offset error)
a4=sigma input as fwhm in units of x
a5=Ampl 2nd guassian
a6=Mean
a7=sigma input as fwhm
a8 =Ampl 3rd gaussian
a9 =Mean
a10=sigma input as fwhm
The variable ngauss determines how many gaussians to fit for. This
routine is similar to fitngauss but it does not have the coma
parameter (which tends to diverge).
(See /pkg/rsi/local/libao/phil/gen/fitngaussnc.pro)
NAME:
fitngaussncfunc - function for fitting n gaussians
SYNTAX:
fitngaussncfunc,x,a,f,pder
ARGS:
x[npts] : independant variable
a[2+3*ngauss] : parameters to fit for
f[npts] : return value of function here
pder[2+3*ngauss] :return partial derivatives with respect to each param
DESCRIPTION:
evaluate the function and optionally it's partial derivatives
for the curvefit routine of idl:
f(x)= a0 + a1*x + a2*exp[ (x-a3)**2 )]
-------------------- + n more gaussians
a4^2
a0= constant
a1= linear in x
a2=Ampl
a3=Mean
a4=sigma
a5=Ampl
a6=Mean
a7=sigma
a8=Ampl
a9=Mean
a10=sigma
note that we are fitting for sigx not sigx^2
(See /pkg/rsi/local/libao/phil/gen/fitngaussnc.pro)
NAME:
fitsin - fit to Asin(Nx-phi) where N=1 to 6. or 123
SYNTAX:result=fitsin(x,y,N,cossin=cossin)
ARGS:
x - independent var. (x values for fit). should already be in radians
y - measured dependent variable.
N - 1..6 .. integral period to fit . 1 to 6.
- 123 . fit 1,2,3 az terms at once
KEYWORDS:
cossin: if set then return amplitudes of cos,sin rather then sin,phase
RETURN:
result[]: float
[0] constant coefficient
[1] amplitude
[2] phase in radians
If cossin is set then return:
[0] constant coefficient
[1] amplitude cosine
[2] amplitude sin
if n=123 then:
[0] constant
[1,2] amp,phase 1 az
[3,4] amp,phase 2 az
[5,6] amp,phase 3 az
DESCRIPTION:
Do a linear least squares fit (svdfit) to a sin wave with integral
values of the frequency (1 through 6 are allowable values). It also
supports simultaneous fit for 1,2,3az.
Returns the coefficients of the fit.
NOTES:
Asin(Nt-phi)= Asin(Nt)cos(phi) - Acos(Nt)sin(phi) = Bsin(Nt) + Ccos(Nt)
B=Acos(phi)
C=-Asin(phi)
phi = atan(sin(phi)/cos(phi))/ = atan(-c,b)
amplitude=sqrt(B^2+C^2)
so the fit for B,C is linear.
result from svd:
a[0] - constant
a[1] - sin coef
a[2] - cos coef
(See /pkg/rsi/local/libao/phil/gen/fitsin.pro)
NAME:
fitsineval - evaluate the fitsin() fit
SYNTAX: y=fitsineval(x,coefAr,cossin=cossin)
ARGS:
x[n]:float evaluate the fit at these x data points. Units should be
radians.
coefAr[m]:float coeff returned from finsin
KEYWORDS:
cossin : if set then fitsin was called with /cossin, so the
coef are cos sin, order rather than amp,phase
RETURNS:
y[n]:float fit evaluated at the x points.
DESCRIPTION:
fitsin() will fit a function: c(0) + c(2*I)*sin(I*az - c(2i+1)) .. i=1,2,3,4.. 6 and
return the coef of the fit in a float array.
This routine will evaluate the fit at the positions specified y x. X must be in radians.
If the call to fitsin( /cossin) had cossin set , then you must also set that keyword
in the call to fitsineval().
(See /pkg/rsi/local/libao/phil/gen/fitsineval.pro)
NAME:
fitsinn - fit to terms of Asin(nx-phi)
SYNTAX:coefI=fitsinn(x,y,N)
ARGS:
x - independent var. (x values for fit). should already be in radians
y - measured dependent variable.
N - number of terms to fit 1 .. N
RETURN:
coefI: {} hold the fit info .. see belo
yfit: float optionally evalute and return the fit
DESCRIPTION:
Do a linear least squares fit (svdfit) to a sin wave with integral
values of the frequency 1 through N. Return the coeficients in coefI.
EXAMPLE:
let az,azerr be the azimuths in deg and the azimuth errors.
Fit a0 +a1*sin(az-ph1) + a2*sin(2az-ph2) + .. an*sin(naz-phn)
N=3 ; use sin(az, 2az,3az)
coefI=fitsinn(az*!dtor,azerr,N,yfit=yfit)
help,coefI
** Structure <8bd0e6c>, 5 tags, length=124, data length=122, refs=1:
N INT 3 .. number we fit to
C0 DOUBLE 0. .. the constant term
PHDEG DOUBLE Array[3] .. phase indegrees
AMP DOUBLE Array[3] .. amplitue of fit
PHRD DOUBLE Array[3] .. phase in radians
COSSIN DOUBLE Array[2, 3] .. cos,sin Amplitudes for each order
NOTES:
A*sin(Nt-phi)= A*sin(Nt)cos(phi) - A*cos(Nt)sin(phi) = B*cos(Nt) + c*sin(Nt)
B=-A*sin(phi)
C=A*cos(phi)
phi = atan(sin(phi)/cos(phi))/ = atan(-B,c)
amplitude=sqrt(B^2+C^2)
so the fit for B,C is linear.
(See /pkg/rsi/local/libao/phil/gen/fitsinn.pro)
NAME:
fitsinneval - evaluate the fitsinn fit
SYNTAX: y=fitsinneval(x,coefI)
ARGS:
x[n]:float evaluate the fit at these x data points. Units should be
radians.
coefI:{} structure returned from fitsinn
RETURNS:
y[n]:float fit evaluated at the x points.
DESCRIPTION:
fitsinn() will fit a function: ai*cos(i*x) + bi*sin(i*x) where x=
0..n. It returns a structure containing the coef of the fit. This routine
fitsinneval() will evaluate this fit at user specified x coordinates (which
should be in radians).
(See /pkg/rsi/local/libao/phil/gen/fitsinneval.pro)
NAME:
fitsinnl - nonlinear least squares fit to a sin
SYNTAX: fitsinnl,x,y,weights,a,yfit,chisq,iter,sigma,tol=tol,itmax=itmax
ARGS:
x[npts]: float x values. The unit will determine the frequency
unit (eg. sec--> cycle/sec, min--> cycle/min).
y[npts]: float measured data values
a[4] : float initial values for fit:
a[0] constant,
a[1] amplitude
a[2] frequency cycles/(xunit)
a[3] phase (fraction of a cycle)
weights[npts]:float weights for each data point
KEYWORDS:
tol : float tolerance for conversion. default 1e-3
itmax: int max number of step iterations before non-convergence.
RETURNS:
a - fitted coefficients
yfit - fit value evalutated at x
chisq - chisq for fit
iter - number of iterations that were made
sigma[4] - sigmas for each coef.
DESCRIPTION:
Do a non-linear least squares fit to a sin wave fitting for offset,
amplitude, frequency, and phase. Use the idl routine curvefit().
(See /pkg/rsi/local/libao/phil/gen/fitsinnl.pro)
NAME:
flag - flag a set of vertical lines
SYNTAX: flag,x,yr=yr,_extra=e
ARGS:
x[n]: float x values where the vertical lines should be drawn
KEYWORDS:
yr[2]: float min,max y values for the veritical lines. The defalut is to
cover the entire y range.
_extra: Any extra keywords are passed to the oplot routine. This
can be color=, linestyle=, etc..
DESCRIPTION:
Draw 1 or more vertical lines on the plot. By default the entire
y range is used. You can limit that with the yr= keyword. Any keywords
accepted by oplot can also be passed in.
EXAMPLE:
After plotting, x vs y, draw two vertical dashed lines at x=15 and x=23.
plot,x,y
flag,[15,23],linestyle=2
(See /pkg/rsi/local/libao/phil/gen/flag.pro)
NAME:
flux12m - flux for some 12meter sources
SYNTAX: istat=flux12m(srcNm,freqar,fluxI,yr=yr,list=list,cat=cat)
ARGS:
srcNm: string source to use. will match any case.
freqAr: fltarr[n] freq to evaluate flux at (MHz)
KEYWORDS:
list : int if not zero the just list the source we
know.
yr : long year for flux. default current year.
this is used if baars et al and there is
a correction for variabilty eg.. casa
cat : string catalog to use: P -perley and butler
B -baars et al (if available)
w -wmap combined:7yrwmap cal src. weiland 2010.
RETURNS:
istat : 1 ok, -1 error
fluxI[n]: {} array of fluxes (Jy)
fluxI:{ srcNm: string,$
n : n ,$; number of fluxes we return
cat : '' ,$; 'p'-perly, 'b' baars et al.
freq:fltarr(n),$ ; freq in Mhz
flux :fltarr(n)} ; flux in Jy
DESCRIPTION:
Compute the flux (JY) at the requested freq for some
sources. Data is taken from perley and butler 2017 or baars et al.
If you request a catatlog and we only have the other, then
you get the one we have (and flux.cat will be set accordingly).
currently just casa and cygnus A
(See /pkg/rsi/local/libao/phil/gen/flux12m.pro)
NAME:
fluxfitvla - return source flux using vla formula
SYNTAX: flux=fluxfitvla(coef,freq,newfit=newfit)
ARGS:
coef[2,3,4...]: float coefficients for fit.
freq: float Frequency in Mhz.
KEYWORD:
newfit : if set then force using the newer vla fit format(aka 1999)
The new fit coef are for Ghz, not Mhz
RETURNS:
flux: float flux in Janskies.
DESCRIPTION:
Evaluate the function
log(s)= coef[0]+coef[1]*log(freq)+coef[2]*(log(freq)^2) + coef[3]*log(freq)^3...
This can be used for the standard calibrators: 3C286,3C48, and 3C147.
You must input the coefficients for each source.
It is taken from the vla calibration manual (chris salters copy
or the web).
There are no corrections to baars et scale done here.
If number of coef > 3 then the newer fit (circa 1999) will be used
that includes freq^3,4... terms).
(See /pkg/rsi/local/libao/phil/gen/fluxfitvla.pro)
NAME:
fluxkuehr - compute flux given kuehr et al. coefficients
SYNTAX: flux=fluxkuehr(coef,freqMhz)
ARGS:
coef[4]: float kuehr coefficients .. see below
freqM : frequency in Mhz
RETURNS:
flux : in janskies
DESCRIPTION:
Return the flux computed from the kuehr coef..
(Kuehr et al., A+AS, 45, 367, (1981))
log10(flux)= coef[0] + coef[1]*x + coef[2]*exp(-coef[3]*x)
and x=log10(freqMhz)
(See /pkg/rsi/local/libao/phil/gen/fluxkuehr.pro)
NAME:
fluxsrc - return source flux
SYNTAX: flux=fluxsrc(srcname,freq,radec=radec,width=width,riseset=riseset)
ARGS:
srcname : string source to search for in files.
freq : float the frequecy of interest in Mhz.
RETURNS:
flux : float The flux in janskies.
radec[2]: double return the ra,dec B1950 in the radec keyword if
provided.
radec[0] (ra) is in hours.
radec[1] (dec) is in degrees.
width[2]: float width (major,minor) axis in arcseconds.
0--> no data
< 0 --> true values are less than these values
riseset[2]: float hhmmss.s rise,set time for source
DESCRIPTION:
The file aodefdir()/data/fluxsrc.dat has fits of flux versus frequency
for a number of sources. This routine will return the flux given a
frequency if the source name is in the file. If you provide the
radec keyword, the the B1950 ra,dec position will be provided.
SEE ALSO:
fluxsrclist, fluxsrcload.
(See /pkg/rsi/local/libao/phil/gen/fluxsrc.pro)
NAME:
fluxsrclist - return list of source names in flux file
SYNTAX: fluxnames=fluxsrclist(print=print,freq=freq,size=size,bnames=bnames,$
riseset=riseset,sortRise=sortRise,$
srclist=srclist,all=all,exclvar=exclvar,$
retall=fluxstr)
ARGS: none
KEYWORDS:
print: if set then print to stdout the source info for the sources
selected. This inclues the name, coefficients, fit rms,
code and comments.
freq : float frequency in Mhz. If provided then evaluate the flux at
this frequency for every source returned. Return the
source names and flux in Janskies.
all : if set then return all available info on each source.
name (flux if Freq specified) CodeForSrc size/comments
size : if set then include any size/comments for the sources
bnames: if set only include sources that start with B.
riseset: if set then include rise,set times in list
sortRise: if set then sort by rise time.. default is name
srclist[]: string If supplied, then only supply info on the sources in
this list that also meet the other criteria (bnames..).
exclVar: if set then exclude variable sources (fit type = 2)
retall: {} fluxstr. Return the entire flux structure for all of the
selected sources.
RETURNS:
A string array containing an entry for each source.
the first entry of every source is the name:
ret[0,*]= srcname
if freq keyword is present, the flux at this frequency will
follow the srcname
ret[1,*]= flux
if size keyword is set then any comments will follow the freq (or name
if no freq requested)
ret[2,*]= comments or
ret[1,*]= comments (if no freq keyword)
DESCRPTION:
Return a list of the sources in the flux file (created by chris salter).
If the /print keyword is supplied then the routine will also write the
source info to stdout.
By default the source name for all sources is returned. You can modify
what is returned with the following keywords:
freq: Evaluates and adds the flux for each source at freq(Mhz).
size: adds the size/comments field in the returned data
all : returns srcName (flux if frq supplied) code size/comments
bnames: Only return sources that start with B. The file has multiple
names for the same source (J,B, and 3C).
riseset: Include rise/set times hhmmss.s in output
srclist: instead of searching the entire file, only look at the sources
in this list.
EXAMPLES:
1. dat=fluxsrclist() .. dat[] is a string array of all the sources
print,dat[0:3]
B0010+005 3C5 B0017+154 3C9
2. dat=fluxsrclist(/print) .. the file is also listed to std out:
prints:
codes:1-good calibrator,2-lousy,3-flux.cat ??
B0010+005 coef: 2.77 -0.82 0.00 rms: 9.00 code:1 ;3C5; Size~27"; Sp(80-5000)
3. dat=fluxsrclist(freq=1420).. dat[2,*] string array
dat[0,*] holds the source names.
dat[1,*] holds the flux (in Jy) for each src
at 1420Mhz
print,dat[*,0:1]
3C132 3.21
3C138 9.14
4. dat=fluxsrclist(freq=1420,/size).. dat is a 2-d string array.
dat[0,*] holds the source names.
dat[1,*] holds the flux for each source
dat[2,*] holds comments for this source
print,dat[*,0:1]
3C132 3.21 3C132 C27.6(756)
3C138 9.14 3C138 Fit to all S's
5. srclist=['B2223+21','3C138']
dat=fluxsrclist(freq=1420,/all,srclist=srclist).. dat[4,3] stringarray
print,dat
B1607+268 4.47 1 cut-off below 1 GHz
B1615+212 1.80 1
B1622+238 2.63 1 3C336
NOTE:
1. There are multiple names for the same source in the flux file (eg
B2314+038,3C459,J2316+040).
2. The code field is: 1 - good flux calibrator, 2-bad flux calibrator,
3 - source from kuehr et al.Chris has not yet
tried to fit this source.
The flux file is located at aodefdir()/data/fluxsrc.dat. At AO, aodefdir() is
/pkg/rsi/local/libao/phil/.
SEE ALSO:
fluxsrc, fluxsrcload.
(See /pkg/rsi/local/libao/phil/gen/fluxsrclist.pro)
NAME:
fluxsrcload - load source flux into common block
SYNTAX: fluxsrcload,retdata,file=file,salter=salter
ARGS:
retdata[]:{fluxsrc} if supplied, then return all the flux structures
in retdata.
KEYWORDS:
file : string filename that holds the source fluxes. Default
is aodefdir() + 'data/fluxsrc.dat
salter : if set then input file from chris salter's file. This
file is used to generate the fluxsrc.dat file. The option
is used when checking that the files are in sync (see chkflux
below).
DESCRIPTION:
Read in all of the source fluxes from the fluxsrc.dat file into
the common block fluxcom. The user can optionally specify
another file to use for the source fluxes.
The routine fluxsrc() computes fluxes from data in this common block.
It will automatically call fluxsrcload() if the common block is not
initialized.
The routine aodefdir()/data/chkflux.pro is used to keep the two files
in sync (chris's file and fluxsrc.dat).
SEE ALSO:
fluxsrc(), fluxsrclist()
(See /pkg/rsi/local/libao/phil/gen/fluxsrcload.pro)
NAME:
foldtmseries - fold a time series
SYNTAX: df=foldtmseries(d,smpTm,period,nbins,stphase=stphase,cnts=cnts)
ARGS:
d[m,n]: float data to fold. If 1d then just fold the time series
if 2-d then each bin of the folded data will have elements
smpTm: float seconds per sample in d
period: float period in secs to fold to.
nbins: long number bins after folding
KEYWORDS:
sTPhase: float the starting second for d[0]. The default is 0.
RETURNS:
df[m,nbins] :double folded data, normalized to the number of counts
cnts[nbins]:long number of counts in each phase bin.
DESCRIPTION:
Fold the time series of d into an array df[nbins]. If d is 2d, then the
the data is folded on the 2nd dimension (each bin of the folded output will have
m entries). The folded data is normalized to the number of counts in each bin.
Optionally return the number of counts for each bin in cnts[nbins].
The folding computations are done in double precision for more accuracy.
EXAMPLE:
; Take a wapp pulsar datafile on B1937+21 and:
;1. input the file
;2. dedisperse it.
;3. fold it
file2='/share/wapp21/p2175.B1937+21.wapp3.53889.0001'
period=1.55772007d-3 ; pulsar period
dm=71.040 ; dm
smpTm=80d-6 ; sample time of the data
npts=sp_dedisp(file2,dm,dedisp2) ; input and dedisperse
nbins=32 ; bins to fold into
dfold=foldtmseries(dedisp2,smpTm,period,nbins) ; fold it
(See /pkg/rsi/local/libao/phil/gen/foldtmseries.pro)
NAME: fromsecs1970 - convert from unixsecs 1970, to ymd,hms SYNTAX: [hhmmss,yyyymmdd]=fromsecs1970(secs1970,utc=utc,daynoyr=daynoyr,str=str) ARGS: secs1970: long seconds from 1970 KEYWORDS: utc: is set then return time as utc. def: local DESCRIPTION: Convert from unix seconds 1970 to hhmmss, yyyymmdd. If /utc is set then return as utc (the default is local time). if the keyword daynoyr= is supplied then also return dayno,year in the daynoyr variable if the str= keyword is supplied then also return time in standard systime string format RETURNS: [hhmmss,yyyymmdd]: lonarr time [dayno,year] : dblarr dayno will include the fraction of day. str : string standard systime string of time Note:secs1970 are always utc. the idl systime routine will convert this to local time.
(See /pkg/rsi/local/libao/phil/gen/fromsecs1970.pro)
NAME: fwhm2tosig2f - FWHM^2 to sigma^2 factor. SYNTAX: fac=fwhm2tosig2f() ARGS: none RETURNS: fac: float conversion factor. DESCRIPTION: for a gaussian defined as: y=A0*exp( -[(x-x0)/sig]^2) convert full width half max squared to sigma squared. If the gaussian is defined as: y=A0*exp( -.5*[(x-x0)/sig]^2) you need to multiply the returned value by .5
(See /pkg/rsi/local/libao/phil/gen/fwhm2tosig2f.pro)
NAME:
fwhmtosigf - convert factor fwhm to sigma
SYNTAX: scl=fwhmtosigf(div2=div2)
KEYWORDS:
div2: if set then return 1/(sqrt(8*log(2))
DESCRIPTION:
Returns factor to go from fwhm to sigma for a
gaussian of: y=a0*epx(x^2/sig^2)
If div2 keyword is set then return factor for:
y=a0*exp(.5* x^2/sig^2)
Example:
let y=exp( - (x/sig)^2)
for fwhm y=.5
if fwhm=x0 then x=x0/2 is the one sided value
so log(.5)= - (x/sig)^2
x=sig * sqrt(-log(.5)) .. but log(1/2)=-log(2)..
so
x =sqrt(log(2))*sig
x0=C*sig where c = .5*sqrt(log(2))
To convert fwhm to sig units you multiply by 1/C
This can be used in guass fitting..
After fitting, you get the widths in sigmas.. got get back
fwhm, you divide by scl
sigma/scl=fhwm
(See /pkg/rsi/local/libao/phil/gen/fwhmtosigf.pro)
NAME:
gainget - return telescope gain(az,za,freq) for rcvr.
SYNTAX: stat=gainget(az,za,freq,rcvrNum,gainVal,date=date,zaonly=zaonly)
ARGS:
az[n]: float azimuth in degrees
za[n]: float zenith angle in degrees
freq: float freq in Mhz.
rcvrNum: int receiver number 1..16 (see helpdt feeds))
KEYWORDS:
date[2]: int [year,daynumber] to use for gain computation. The default
is to use the current day. If they gain curves change
with time, this allows you to access the gain curve that
was in use when the data was taken.
zaonly: if set then only return the zenith angle dependence (averages
over azimuth)
RETURNS:
gainval[n]: float gain in Kelvins/Jy
stat: int -1 --> error, not data returned
0 --> requested outside freq range of fits. Return gain
of the closed frequency value.
the closest value.
1 --> frequency interpolated gain value returned.
DESCRIPTION:
Return the gain (K/Jy) for the requested receiver at the specified
frequency and az,za. The default is to use the current gain values. The
date keyword allows you to access a gain curve that was valid at some
other epoch.
Fits have been done for g(az,za) at different frequencies for various
receivers. This routine will input the fit information and compute the
gain for the two closest frequencies and then interpolate to the
requested freqeuncy. The input fit data is stored in a common block so
the data does not have to be input from disc a second time unless you
pick a different receiver.
NOTE:
Some receivers have no gain fits. They will return -1 in the status.
If a requested frequency is outside the fitted values, then the value
at the closest frequency is returned (no extrapolation is done).
If you have correlator data, you can use corhgainget() to get the
gain value. It will figure out the receiver number and date from the
header and then call gainget.
For a description of the gain calibration see:
http://www.naic.edu/~phil. Look under calibration for the receiver
of interest. The lines with the remark (gain curves) were used to compute
the gain curves.
EXAMPLES:
lbw=5
get gain at 1400Mhz az=120,za=10 for lbw
stat=gainget(120.,10.,1400.,lbw,gainval)
az=fltarr(20) ; az = 0 degrees.
za=findgen(20)+1 ; za=1..20
date=[2001,200] ; for 2001, daynumber:200
stat=gainget(az,za,1321,lbw,gainval,date=date)
gainval will be an array of 20 values for za 1 to 20 degrees and azimuth
of 0 degrees.
to convert from daynumber to day,month,year
daynum=dmtodayno(d,mon,year)
dm =daynotodm(daynum,year)
where dm=[day,month]
SEE ALSO:gaininpdata, calinpdata, corhcalval
(See /pkg/rsi/local/libao/phil/gen/gainget.pro)
NAME:
gainget12m - return 12m telescope gain(el,freq) for rcvr.
SYNTAX: stat=gainget12m(el,freq,rcvrNum,gainVal,date=date)
ARGS:
el[n]: float elevation angle in degrees for gains
freq: float freq in Mhz.
rcvrNum: int receiver number 7 or 11
KEYWORDS:
date[2]: int [year,daynumber] to use for gain computation. The default
is to use the current day. If they gain curves change
with time, this allows you to access the gain curve that
was in use when the data was taken.
RETURNS:
gainval[n]: float gain in Kelvins/Jy
stat: int -1 --> error, not data returned
0 --> requested outside freq range of fits. Return gain
of the closed frequency value.
the closest value.
1 --> frequency interpolated gain value returned.
DESCRIPTION:
Return the gain (K/Jy) for the requested receiver at the specified
frequency and elevations. The default is to use the current gain values. The
date keyword allows you to access a gain curve that was valid at some
other epoch.
Fits have been done for g(el) at different frequencies for various
receivers (currently only xb). This routine will input the fit information and compute the
gain for the two closest frequencies and then interpolate to the
requested freqeuncy. The input fit data is stored in a common block so
the data does not have to be input from disc a second time unless you
pick a different receiver.
NOTE:
Some receivers have no gain fits. They will return -1 in the status.
If a requested frequency is outside the fitted values, then the value
at the closest frequency is returned (no extrapolation is done).
For a description of the gain calibration see:
http://www.naic.edu/~phil. Look under calibration for the receiver
of interest. The lines with the remark (gain curves) were used to compute
the gain curves.
EXAMPLES:
SEE ALSO:gaininpdata12m,
(See /pkg/rsi/local/libao/phil/gen/gainget12m.pro)
NAME:
gainInpData - input gain data for rcvr.
SYNTAX:
istat=gainInpData(rcvNum,gainData,fname=fname,date=date)
ARGS:
rcvNum: 1 thru 16. receiver to use (same as hdr.iflo.stat1.rfnum).
KEYWORDS:
fname: to specify an alternate data file with gain values.
The default file is aodefdir() + 'data/gain.datR{rcvNum}
date : [year,daynum] .. if specified the data when you want the
gain for..default is most recent.
RETURNS:
istat: 1 ok, -1 bad file/rcvnum,data probably no fit data available.
gainData[n]:{gainData} return fit info. 1 structure for each frequency
DESCRIPTION:
Input the gain fit data for all the frequncies of a particular receiver
(rcvNum). The rcvNum can be extracted from the headers with iflohrfnum().
The default datafile is aodefdir() + 'data/gain.datR{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,
- 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
freq fitType c0 c1 c2 .... cN pol calVala calvalB
c0..cN are the fit coefficients, pol is I or A or B, CalValA,B are the
cal values used for this fit.
The structure format for {gaindata} is:
gainData.rcvNum receiver number
gainData.numFreq number of frequencies found
gainData.startYr for the fit
gainData.startDaynum for the fit
gainData.endYr for the fit
gainData.endDaynum for the fit
gainData.fitI[numFreq] {azzafit} structure holding the coef and other
info for each fit
gainData.calVal[2,numFreq] cal values used when each fit was made.
See azzafit for a description of the {azzafit} structure.
This routine is called automatically by corhgainval and gainget().
How the different cal routines vary:
gaininpdata() inputs the data from disc. You must specify the rcvrnum.
It defaults to the current date. It loads a table in common
holding the fit info for all of the frequencies measured.
gainval() Pass in the frequency and the rcvrnum. It will input the
data using gaininpdata if necessary, do the interpolation
and return the gain.
corhgainval() You specify the correlator sbc header (eg b.b1.h). It will
compute the frequency and then call gainval(). It returns the
gain value.
You can call the azzafitpr,eval routines fit the gainData.fitI[]
to plot out the fit fitvalues.
SEE ALSO: corhgainval, gainget ,azzafit, azzafiteval,azzafitpr
(See /pkg/rsi/local/libao/phil/gen/gaininpdata.pro)
NAME:
gainInpData12m - input gain data for 12meter rcvrs.
SYNTAX:
istat=gainInpData12m(rcvNum,gainData,fname=fname,date=date)
ARGS:
rcvNum: 7 or 11 sband or xband (prior to wbrcvr)
KEYWORDS:
fname: to specify an alternate data file with gain values.
The default file is aodefdir() + 'data/gain12m.datR{rcvNum}
date : [year,daynum] .. if specified the data when you want the
gain for..default is most recent.
RETURNS:
istat: 1 ok, -1 bad file/rcvnum,data probably no fit data available.
gainData[n]:{gainData12m} return fit info. 1 structure for each frequency
DESCRIPTION:
Input the gain fit data for all the frequncies of a particular receiver
(rcvNum). The rcvNum can be extracted from the headers with iflohrfnum().
The default datafile is aodefdir() + 'data/gain12m.datR{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,
- 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
freq fitType A0 A1 .. sigFit
c0..cN are the fit coefficients, pol is I or A or B,
The structure format for {gaindata} is:
gainData.rcvNum receiver number
gainData.numFreq number of frequencies found
gainData.startYr for the fit
gainData.startDaynum for the fit
gainData.endYr for the fit
gainData.endDaynum for the fit
gainData.fitI[numFreq] {azelfit12m} structure holding the coef and other
info for each fit
{azelfit12m}:
{numcoef : 10,$; for fittype
fittype : 1,$;
freq : 0.,$;Mhz
coef : dblarr(10),$; 10 is max coef
sigma : 0.d ,$; of fit
rfnum : 0 ,$; 7 or 11
pol : ' ' ,$; I, A, or B
type : '' ,$; ascii string for type: gain,tsys,sefd..
date :' ' } ;date fit run
This routine is called automatically by gainget12m().
How the different cal routines vary:
SEE ALSO: gainget12m ,masgainget12m
(See /pkg/rsi/local/libao/phil/gen/gaininpdata12m.pro)
NAME:
getscanind - get indices for start of each scan
SYNTAX: getscanind,scanlist,scanind,scanlen
ARGS :
scanlist[len]: long array of scan numbers
RETURNS:
scanind[len] : long indices into scanlist for start
of each scan.
scanlen[len] : long number of entries in each scan.
DESCRIPTION:
The routine corpwr() returns power information for each record in
a file. This includes the scan number of each record. Typically there
will be many records in a scan. This routine will search the array of
scan numbers and return the starting index for the start of each scan
and the number of records in the scan. It does this by scanning the array
and looking for where the scan number changes.
EXAMPLE:
print,corpwr(lun,9999,p) ... up to 9999 recs
getscanind,p.scan,scanind,scanlen
; now loop thru each scan returned
nscans=(size(scanind))[1]
for i=0,nscans-1 do begin
; grab those belonging to 1 scan
p1=p[scanind[i]:scanind[i]+scanlen[i]-1]
...process
endfor
(See /pkg/rsi/local/libao/phil/gen/getscanind.pro)
NAME:
getscanindx - extract scan from array.
SYNTAX: subarr=getscanindx(datarr,scanind,scanlen)
ARGS :
datarr[]: any type extract subarray from here
index : long .. index to extract
scanind[len]: long .. index array from getscanind
scanlen[len]: long .. len each scan from getscanind
RETURNS:
datsubarr[] : anytype.. ith scans data
of each scan.
DESCRIPTION:
Some routines return multiple scans data as one large array with
one of the elements in the array being the scan number. getscanind()
will find the start and length of each scan in this array.
getscanindx will extract the i'th scans data from the array.
EXAMPLE:
print,corpwr(lun,9999,p) ... up to 9999 recs
getscanind,p.scan,scanind,scanlen
; now loop thru extracting each scans data
nscans=(size(scanind))[1]
for i=0,nscans-1 do begin
; grab those belonging to the ith scan
p1=getscanindx(p,i,scanind,scanlen)
...process
endfor
(See /pkg/rsi/local/libao/phil/gen/getscanindx.pro)
NAME:
getsl - scan a corfile and return the scan list.
SYNTAX: sl=getsl(lun,scan=scan,maxscans=maxscans,search=search)
ARGS:
lun: int assigned to open file
KEYWORDS:
scan: long if supplied then start scanning from this scan number.
default is the beginning of the file
maxscans: long if supplied then quit after this many scans are found
default: 5000L
search: if set then search for 1st header
RETURNS:
sl[nscans]:{sl} holds scan list info
DESCRIPTION
This routine reads a corfile and returns an array of scanlist structures.
This array contains summary information for each scan of the file:
sl.scan - long scan number for this scan
sl.bytepos - unsigned long for start of this scan
sl.stat - .. not used yet..
sl.rcvnum - byte receiver number 1-16
note: ch is rcvnum 100
sl.numfrq - byte .. number of freq,cor boards used this scan
sl.rectype - byte 1 -calon
2 -caloff
3 -onoff on pos
4 -onoff off pos
5 -coron (track just on position)
6 -cormap1
7 -cormapdec
8 -cordrift
9 -corcrossch
10 -x111auto (rfi monitoring)
11 -one murrays on with calon/off at the null
12 -onoffbml murrays on, off and cal at null
13 -corcrosschL 8 strips instead of 4
sl.numrecs - long .. number of groups(records in scan)
sl.freq[4] float- topocentric frequency center each subband
sl.julday double- julian date start of scan
sl.srcname string - source name (max 12 long)
sl.procname string - procedure name used.
Some routines can use the sl structure to perform random access to
files (bypassing the need to search for a scan). The sl[] array can
also be used with the where() command to rapidly extract subsets of a
file.
EXAMPLE:
openr,lun,'/share/olcor/corfile.02nov00.x101.1',/get_lun
sl=getsl(lun)
1. process all of the lband wide data in a file:
ind=where(sl.rcvnum eq 5,count)
for i=0,n_elements(ind)-1
print,corinpscan(lun,b,scan=sl[ind[i]].scan,sl=sl)
.. process
endfor
2. Find all of the on/off patterns in a file. Make sure that the
number of records in the on equals the number in the off.
indon=where(sl.rectype eq 3 ,count)
if count le 0 then goto,nopairs
; make sure an off follows the on and has the same number of records..
; (actually this will fail if the last rec of the file is an on since
; indon+1 will go beyond the end of the sl array..)
ind= where((sl[indon+1].rectype eq 4) and $
(sl[indon].numrecs eq sl[indon+1].numrecs),count)
if count le 0 then goto,nopairs
indon=indon[ind]
3. plot all of the cal on/off records in a file with cormon().
ind=where(sl.rectype le 2)
cormon,lun,sl=sl[ind]
Note this will not work with files > 2gigabytes since it is
using a 32 bit integer.
(See /pkg/rsi/local/libao/phil/gen/getsl.pro)
NAME:
gs - generate a gaussian
SYNTAX: f=gs(len,height,fwhm,position)
ARGS:
len: int .. length to make f. x values will be 0 thru len-1
height: float.. height of the gaussian.
fwhm : float.. full width at half maximum. full range 0 to len-1
position: float.. position for the peak (0..len-1)
RETURNS:
computed gaussian as a double.
(See /pkg/rsi/local/libao/phil/gen/gs.pro)
NAME:
gs2d - generate a 2d gaussian
SYNTAX: f=gs2d(len,height,fwhm)
ARGS:
len: int .. length one side. x values will be 0 thru len-1
make an odd value for center to be located
at center of array.
height: float.. height of the gaussian.
fwhm : float.. full width at half maximum. full range 0 to len-1
RETURNS:
xy(2,len,len): double x,y coorinate for each point
computed 2d gaussian as a double.
(See /pkg/rsi/local/libao/phil/gen/gs2d.pro)
NAME:
gseval - evaluate a gausian at the requested positions.
SYNTAX: val=gseval(fwhm,position)
ARGS:
fwhm : float.. full width at half maximum of gaussian.
position: float. to evaluate (same units as fwhm
RETURNS:
vals : evaluated at position.
Assume gaussian is unit height
(See /pkg/rsi/local/libao/phil/gen/gseval.pro)
NAME:
gsfit2d - cross pattern 2d fit to total power az,za stripsI
SYNTAX: fitCoef=gsfit2d,az,za,z,initCoef,linearza=linearza,zfit=zfit,$
covar=covar,sigCoef=sigCoef,chisq=chisq,sigma=sigma,$
trouble=trouble,weights=weights
ARGS:
az[npts]: az pos (offsets from center) in gcdeg all pnts in pat
za[npts]: za pos (offsets from center) in gcdeg all pnts in pat
z[npts]: measured data points
initCoef[m]: float. coef to fit for. first guess
fitCoef[m]: float coef from fit
KEYWORDS:
linearza : if true then include a linear term in za (a7)
zfit[npts]: float return fit evaluated at input points
sigma=sigma : sigma for fit
covar[10,10]: covariance matrix
sigCoef[10] : sigmas for coef.
chisq : float ..
weights[npts]: float weights
trouble : 0 converged,
-1 chisq infinite,
-2 flamdacount>30*10/flstep
-3 iteration > iterationmax def. 20
-4 alpha/c not finite.. probably 0 in partial deriv
nostop : if set then if trouble=-4 then don't stop, just return
with trouble=-4
DESCRIPTION:
Fit to the function:
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2 - yp^2/sigyp^2]
If linearza is set then fit to:
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2 - yp^2/sigyp^2] + a7*za
You need to give the coef initial values when the routine is called.
The units are:
a0 - constant
a1 : amplitude
a2 : xoffset ,az coord direction offset is az units
a3 : yoffset ,za coord direction offset in za units
a4 : sigx^2 ,in prime coordinate system: fwhm az units
a5 : sigy^2 ,in prime coordinate system: fwhm za units
a6 : theta ,rotate unprimed to primed aligned along ellipsoid of beam,Deg
a7 : zaslope ,The za slope in amplitude units per za unit
xm=(x-x0)
ym=(y-y0)
xp = xm*cos(th) + ym*sin(th)
yp =-xm*sin(th) + ym*cos(th)
angle theta rotates from az,za to axes aligned with the major axis
of the beam elipse
The x,y values are passed via a common block. The call passes in
an index to this common block.
The fit is done in the units passed in (az,za). They should be
greatcircle units with the same scale (don't mix hrs and deg)
(See /pkg/rsi/local/libao/phil/gen/gsfit2d.pro)
NAME:
gsfit2dc - cross pattern 2d fit to az,za stripsI with coma
SYNTAX: gsfit2dc,az,za,z,a,zfit=zfit,covar=covar,sigCoef=sigCoef,$
chisq=chisq,sigma=sigma,trouble=trouble
ARGS:
az[npts]: az pos (offsets from center) in gcdeg all pnts in pat
za[npts]: za pos (offsets from center) in gcdeg all pnts in pat
z[npts]: measured data points
a[m] : float. coef to fit for. first guess
KEYWORDS:
zfit[npts]: float return fit evaluated at input points
sigma=sigma : sigma for fit
covar[10,10]: covariance matrix
sigCoef[10] : sigmas for coef.
chisq : float ..
trouble : 0 converged,
-1 chisq infinite,
-2 flamdacount>30*10/flstep
-3 iteration > iterationmax def. 20
DESCRIPTION:
Fit to the function:
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2*(1.+alphax*xpp) -
yp^2/sigyp^2*(1.+alphay*ypp)]
You need to give the coef initial values when the routine is called.
The units are:
a0 - constant
a1 : amplitude
a2 : xoffset ,az coord direction units:amin
a3 : yoffset ,za coord direction :amin
a4 : sigx^2 ,in prime coordinate system: fwhm amin
a5 : sigy^2 ,in prime coordinate system: fwhm amin
a6 : alphax ,comma in pp system aligned along coma direction
a7 : alphay ,comma in pp system aligned along coma direction
a8 : theta ,rotate unprimed to primed aligned along ellipsoid of beam,Deg
a9 : thetap ,rotate primed to coma aligned, Deg
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2*(1.+alphax*xpp) -
yp^2/sigyp^2*(1.+alphay*ypp)]
xm=(x-x0)
ym=(y-y0)
xp = xm*cos(th) + ym*sin(th)
yp =-xm*sin(th) + ym*cos(th)
xpp= xp*cos(thp) + xp*sin(thp)
ypp=-xp*sin(thp) + yp*cos(thp)
angle theta rotates from az,za to axes aligned with the major axis
of the beam elipse
angle thetap rotates from xp,yp axes to xpp,ypp which are aligned with the
coma direction of the fit.
The x,y values are passed via a common block. The call passes in
an index to this common block.
we convert postions to arcminutes and angles to radians before calling
fit, we then back convert when done
(See /pkg/rsi/local/libao/phil/gen/gsfit2dc.pro)
NAME:
gsfit2deval - evaluate coef returned from gsfit2d
SYNTAX: z=gsfit2deval,x,y,fitcoef
ARGS:
x[npts]: x pos
y[npts]: y pos
fitcoef[m]: float. coef from gsfit2d
RETURNS:
z[npts]: float fit evalutated at x,y
DESCRIPTION:
Evaluate the Fit returned from gsfit2d at the
requested x,y positions
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2 - yp^2/sigyp^2] + a7*za
The units are:
a0 - constant
a1 : amplitude
a2 : xoffset ,az coord direction units:amin
a3 : yoffset ,za coord direction :amin
a4 : sigx^2 ,in prime coordinate system: fwhm amin
a5 : sigy^2 ,in prime coordinate system: fwhm amin
a6 : theta ,rotate unprimed to primed aligned along ellipsoid of beam,Deg
a7 : zaslope ,The za slope in amplitude units per deg za
xm=(x-x0)
ym=(y-y0)
xp = xm*cos(th) + ym*sin(th)
yp =-xm*sin(th) + ym*cos(th)
angle theta rotates from az,za to axes aligned with the major axis
of the beam elipse
The x,y values are passed via a common block. The call passes in
an index to this common block.
we convert postions to arcminutes and angles to radians before calling
fit, we then back convert when done
get the x,y values in the prime system (aligned with the ellipitcal beam).
(See /pkg/rsi/local/libao/phil/gen/gsfit2deval.pro)
NAME:
gsfit2deval_12m - evaluate fitI returned from mascrossfit2d()
SYNTAX: z=gsfit2deval_12m(fitI)
ARGS:
fitI[0]: singleo struct (fitI) returned from mascrossfit2d().
RETURNS:
z[npts,2,nbms]: float fit evalutated at fit points (npnts,az/el)
DESCRIPTION:
mascrossfit2d() fits a 2d gaussian to crosses made on a source with the 12meter.
It uses gsfit2d_12m() do do the actual fitting. mascrossfit2d then returns an array of
structures; fitI[n] (one for each cross done).
This routine gsfit2deval_12m() will evaluate one of these structures at the fit values
note:mascrossfit2d() can return and array of structs if more than one cross is passed in.
this routine will only evaluate one fit at a time.
eg mascrossfit2d(..fitiar)
z=gsfit2deval_12m{fitiar[4]) to evaluate the 5th crosss.
the fit. holding info from the fits. This routine will take a s
The 2d fit is:
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2 - yp^2/sigyp^2] + a7*elOffset
see mascrossfit2d() or gsfit2d_12m() for a description of the coef and xp,yp
The x,y values for the evaluation are taken from fitI that is passed in.
The fiti can have more than one bm (freq band). The returned evaluation
will be dimensioned z[npnts,2,nbms) where 2 is az,el strip.
Example:
1. assume you called mascrossfit2d() with xband data that had
1200 samples/strip 7 bms and it returned fitI and 1 cross.
2. to evaluate the fit:
zret=gsfit2deval_12m(fitI)
zret would be dimensioned: zret[1200,2,7]
3. to overplot the data with the fits:
npnts=fiti.pntstrip ; normally 1200.
nbms=n_elements(fit.fitibm)
xaz=fiti.offard[*,0] ; az offsets
xel=fiti.offard[*,1] ; el offsets
zaz=reform(fiti.fitibm.z[*,0],npnts,nbms) ; [1200,7] az data values
zel=reform(fiti.fitibm.z[*,1],npnts,nbms ; [1200,7] el data values
!p.multi[0,1,2]
stripsxy,xaz,fiti.fitibm.z[*,0],0,0,/step ; az strips data
stripsxy,xaz,zaz,0,0,/step,/over ; overplot the az fits
stripsxy,xaz,fiti.fitibm.z[*,1],0,0,/step ; el strips data
stripsxy,xel,zel,0,0,/step,/over ; overplot el fits
Notes:
1. 14mar23 - an option was added to remove linear baselines prior to the 2d fitting.
This changed the number of elements in fitI[] (to include the baseline info)
if this was selected:
- the z data values in fitI.fitibm[ibm].z have had these linear baselines removed
before the 2d fitting is done.
- if you want the actual offsets (tsys), you need to look at the results from the
linear fits
: fitI.fitibm.blfit.az.coef[2]: coef[0] +coef[1]*azoff . gives the entrire fit
: fitI.fitibm.blfit.el.coef[2]: coef[0] +coef[1]*eloff . ditto
or
: fitI.fitibm.blfit.az.avgVal avg val from fit (tsys)
: fitI.fitibm.blfit.el.avgVal ditto
(See /pkg/rsi/local/libao/phil/gen/gsfit2deval_12m.pro)
NAME:
gsfit2d_12m - cross pattern 2d fit to total power az,el stripsI
SYNTAX: fitCoef=gsfit2d_12m,azOff,elOff,z,el,initCoef,zfit=zfit,$
covar=covar,sigCoef=sigCoef,chisq=chisq,sigma=sigma,$
rot=rot,trouble=trouble,weights=weights
ARGS:
azOff[npts]: az pos (offsets from center of strip) in gcdeg all pnts in pat
elOff[npts]: el pos (offsets from center of strip) in gcdeg all pnts in pat
z[npts]: measured data points
el[npnts]: encoder offset of elevation from central elevation value.
(needed since the el offset along the az strip relative to the source
is zero, but the change in actual elevation is non zero (since the source
is moving in az,el during the cross.
initCoef[m]: float. coef to fit for. first guess
fitCoef[m]: float coef from fit
rot : int if set then do the beam rotation. the new default is no rotation.
the initCoef should still contain 8 elements. initCoef[6] rotation angle
is ignored.
For the 305m the beam was elliptical and rotation was needed
For the 12meter the beam in more spherical, and there is not enough
ellipticity for this to work well.
KEYWORDS:
zfit[npts]: float return fit evaluated at input points
sigma=sigma : sigma for fit
covar[10,10]: covariance matrix
sigCoef[10] : sigmas for coef.
chisq : float ..
weights[npts]: float weights
trouble : 0 converged,
-1 chisq infinite,
-2 flamdacount>30*10/flstep
-3 iteration > iterationmax def. 20
-4 alpha/c not finite.. probably 0 in partial deriv
-5 one of the coef not finite
-6 a width^2 is negative
nostop : if set then if trouble=-4 then don't stop, just return
with trouble=-4
cfplot : for curve fit: 0 noplot,1 plot,nowait, 2 plot wait at last fit, 3 wait each iter
DESCRIPTION:
Fit to the function:
z(x,y)= a0 + a1*exp[-xp^2/sigxp^2 - yp^2/sigyp^2] + a7*elOff
You need to give the coef initial values when the routine is called.
The units are:
a0 - constant
a1 : amplitude
a2 : xoffset ,az coord direction offset is az units
a3 : yoffset ,el coord direction offset in el units
a4 : FWHM_AZ ,in prime coordinate system: fwhm az units
a5 : FWHM_EL ,in prime coordinate system: fwhm el units
a6 : theta ,rotate unprimed to primed aligned along ellipsoid of beam,Deg
(only do rotation if /rot is set
a7 : zaslope ,The el slope in amplitude units per el unit
xm=(x-x0)
ym=(y-y0)
xp = xm*cos(th) + ym*sin(th)
yp =-xm*sin(th) + ym*cos(th)
angle theta rotates from az,za to axes aligned with the major axis
of the beam elipse. Only rotate if /rot is set. If not, the rot angle is forced to 0.
The x,y values are passed via a common block. The call passes in
an index to this common block.
The fit is done in the units passed in (az,el). They should be
greatcircle units with the same scale (don't mix hrs and deg)
(See /pkg/rsi/local/libao/phil/gen/gsfit2d_12m.pro)
NAME:
hansmo - hanning smooth a dataset
SYNTAX: dsmo=hansmo(d)
ARGS:
d[m,n]: float/double data to smooth
RETURNS:
dsmo[m,n]: the smoothed data.
DESCRIPTION:
hansmo will hanning smooth the data in the array d. d can be 1 or
more dimensions. It will smooth m points at a time.
The smoothing is done by convolution in the input domain.
EXAMPLE:
d=fltarr(1024,4)
...
dsmo=hansmo(d)
In this example the smoothing would be:
for i=0,3 do dsmo[*,i]=hanningSmooth(d[*,i])
(See /pkg/rsi/local/libao/phil/gen/hansmo.pro)
NAME: hardcopy - flush the postscript data to disc. SYNTAX: hardcopy ARGS: none DESCRIPTION: Flush the postscript buffers to disc. Call this routine before swithing back to x windows display. SEE ALSO: ps,pscol,psimage, x
(See /pkg/rsi/local/libao/phil/gen/hardcopy.pro)
NAME:
hdrget - input headers
SYNTAX: nhdrs=hdrget(lun,numhdrs,hdrs,scan=scan,std=std)
ARGS:
lun: int assigned to open file
numhdrs: long number of headers to read
KEYWORDS:
scan: long position to scan before listing.
std: if set then just return the standard header
RETURNS:
hdrsd[]:{hdr} return headers here (or hdrStd)
nhdrs :long number of headers found
DESCRPIPTION:
Input numhdrs headers from the current position in the file.
If the scan keyword is used, then position to start of scan before inputting.
If an integration requires more than 1 record (eg 4 correlator boards) then
each record will count as a header.
SEE ALSO:
posscan,scanlist
(See /pkg/rsi/local/libao/phil/gen/hdrget.pro)
NAME: hexpr - hex printout SYNTAX: hexpr,num ARGS: num[]: long output num as hex (1 number per line)
(See /pkg/rsi/local/libao/phil/gen/hexpr.pro)
NAME: hms1_hms3 - convert hour,min,secs 1 word to hour,min,sec separate words SYNTAX - ret=hms1_hms3(hhmmss) ARGS: hhmmss : double angle to convert RETURNS: ret[4] : double hour,min,sec, and sign
(See /pkg/rsi/local/libao/phil/gen/hms1_hms3.pro)
NAME: hms1_hr - convert hhmmss.sss as a double to hours. SYNTAX: angHr=hms1_hr(hhmmss) ARGS: hhmmss: double value to convert RETURNS: angHr: double the angle converted to hours DESCRIPTION Convert packed hours, minutes, seconds to hours. The input is a single double with hhmmss.ss with hh hours, mm minutes, ss.s seconds.
(See /pkg/rsi/local/libao/phil/gen/hms1_hr.pro)
NAME: hms1_rad - convert hhmmss.sss as a double to radians. SYNTAX: angRad=hms1_rad(hhmmss) ARGS: hhmmss: double value to convert RETURNS: angRad: double the angle converted to radians. DESCRIPTION Convert packed hours, minutes, seconds to radians. The input is a single double with hhmmss.ss with hh hours, mm minutes, ss.s seconds.
(See /pkg/rsi/local/libao/phil/gen/hms1_rad.pro)
NAME: hms1_secmid - convert hour,min,secs 1 word to seconds from midnite SYNTAX - ret=hms1_secmid(hhmmss) ARGS: hhmmss[n]: double angle to convert RETURNS: ret[n] : double secMidnite
(See /pkg/rsi/local/libao/phil/gen/hms1_secmid.pro)
NAME: hor - set horizontal scale for plotting. SYNTAX: hor,hor1,hor2 ARGS: hor1: float min horizontal value hor2: float max horizontal value. DESCRIPTION: Load the !x.range system value with the min,max xrange to plot. To reset to auto scaling call hor with no args. SEE ALSO: ver
(See /pkg/rsi/local/libao/phil/gen/hor.pro)
NAME:
ifloh10gchybrid - return true if 10 ghz hybrid in use
SYNTAX: istat=ifloh10gchybrid(iflohdr)
ARGS:
iflohdr:{hdriflo} .. iflo portion of header.
RETURNS:
istat: int 0 hybrid out, 1 hybrid in use
DESCRIPTION:
Return 1 if the hybrid in the 10 ghz upconverter is in use.
(See /pkg/rsi/local/libao/phil/gen/iflohquery.pro)
NAME:
iflohcaltype - return the type of cal used.
SYNTAX: istat=iflohcaltype(iflohdr)
ARGS:
iflohdr[]:{hdriflo} .. iflo portion of header.
RETURNS:
istat: int 0-7 setting for caltype.
DESCRIPTION:
Return the setting of the caltype when this record was written.
the values are:
0 - low correlated cal (1 diode)
1 - high correlated cal (1 diode)
2 - low crossed over cal (2 diodes)
3 - high crossed over cal (2 diodes)
4 - low uncorrelated cal (2 diodes)
5 - high uncorrelated cal (2 diodes)
6 - low correlated 90 deg phase shift cal (1 diode)
7 - high correlated 90 deg phase shift cal (1 diode)
This is the setting of the cal switch. It does not mean that
this is a cal record.
EXAMPLE:
If you have read in a correlator record:
print,corget(lun,b)
istat=iflohcaltype(b.b1.h.iflo)
will return the caltype in istat.
If iflohdr is an array then an array of ints will be returned each with a
cal type.
SEE ALSO:
chkcalonoff
(See /pkg/rsi/local/libao/phil/gen/iflohquery.pro)
NAME:
iflohlbwpol - check if hybrid used on lband wide.
SYNTAX: istat=iflohlbwpol(iflohdr)
ARGS:
iflohdr[]:{hdriflo}
RETURNS:
istat: int 1 circular pol hybrid in, 0 linear pol hybrid out.
DESCRIPTION:
The lband wide receiver has an OMT that provides linear
polarization. After the dewar there is a switchable hybrid that
converts from linear to circular. You need to know this setting for
lbw when you are using the cal values since the cals are injected as
linear and are averaged if the hybrid is inserted.
EXAMPLE:
istat=corget(lun,b)
istat=iflohlbwpol(b.b1.h.iflo)
NOTE: If iflohdr is an array, then an array of ints either 1 or 0.
(See /pkg/rsi/local/libao/phil/gen/iflohquery.pro)
NAME:
iflohrfnum - return the receiver # for this record
SYNTAX: istat=iflohrfnum(iflohdr,hdrpnt=hdrpnt)
ARGS:
iflohdr[]:{hdriflo} .. iflo portion of header.
RETURNS:
istat: int 1-16 receiver number.
DESCRIPTION:
Return the receiver number used for this record. see
helpdt feeds (in sunos) for a list of receiver numbers.
if iflohdr is an array than an array of rfnums will be returned.
(See /pkg/rsi/local/libao/phil/gen/iflohquery.pro)
NAME:
iflohstat - decode status words for iflo
SYNTAX: statInfo=iflohstat(iflohdr)
ARGS:
iflohdr[]:{hdriflo}
RETURNS:
statInfo[]{ifstat} decoded status structure
DESCRIPTION:
The iflo header contains various info that is encoded in bitmaps
(if1.st1,if1.st2,if2.st1,if2.st4). This routine decodes these
bitmaps and returns them in a structure. The input can be 1 or more
iflo headers.
Example:
print,corget(lun,b)
ifstat=iflohstat(b.b1.h.iflo)
The definition of the structure is:
rfnum: 0 ,$; rcvnum 1 to 16
if1num: 0 ,$; 1-300,2-750,3-1500 (2-12)Ghz,4-10000,5strthru
hybridIn10Ghz: 0 ,$; for 10ghz upconverter
lo1HiSid: 0 ,$; 1 yes
lbwLinPol: 0 ,$; 1 lin, 0 circular
syn1rfOn: 0 ,$; 1st lo ,1 yes
syn2rfOn: 0 ,$; sbtx synth,1 yes
lbFbA : 0 ,$; lbw filters (bit map) 1.. 9
lbFbB : 0 ,$; lbw filters (bit map) 1.. 9
useFiber: 0 ,$; 1 yes
calRcvMux: 0 ,$; rcvNUmber for upstairs cal mux
calType : 0 ,$; 0 Lcorcal,1 Hcorcal,2 Lxcal,3 Hxcal,4lcal,
; 5 Hcal,6 L90cal,7 H90cal
ac1Pwrsw : 0 ,$; ac1 strip bits on /off
ac2Pwrsw : 0 ,$; ac2 strip bits on /off
xfer1Sw : 0 ,$; 1 normal, 0 switched
sbnShClosed: 0 ,$; 1 closed
lo2Hiside : 0 ,$; 1--> high side. 4 bits
;
; from if2
;
if2inpFreq : 0 ,$; 0 spare,1 300, 2 750, 3 1500
vlbafrq2ghz: 0 ,$; 1 2000, 0 750
xfer2Sw: 0 ,$; 1 normal, 0 switched
blank430 : 0 ,$; 1 blank 0 no. was sbdoppler
noiseSrcOn: 0 ,$; 1 yes, 0 no
dualPol30If: 0 ,$; 1 2 pol, 0 bands polA
vis30MhzGr : 0 ,$; 1 greg, 0 ch
calTTlSrc : 0 ,$; 1 to 8. cal ttl pulse source
pwrMetToIF : 0 ,$; 1 yes, 0 to front panel
useAlfa : 0 ,$; 1 using alfa, 0 no
sigSrc : 0 ,$; 1 0=gr,1=ch,2=noise
if2Stat4[4]:
synDest: 0 ,$; 0-frontpanel,1-260to30conv,2-vlba/sb,3-mixers
mixerCfr: 0 ,$; 0-750,1-1250,2-1500,3-1750
ampInpSrc: 0 ,$; 0, 1-mixers,2-heliax,3=300Mhz IF
ampExtMask: 0 bit mask 7 outputs. 1->extinp, 0 from
(See /pkg/rsi/local/libao/phil/gen/iflohquery.pro)
NAME:
imgdisp - display a 2-d array as an image.
SYNTAX : imgdisp,d,ret,zx=zx,zy=zy,clip=clip,nsigclip=nsigclip,$
maxiter=maxiter,histeq=histeq,$
/rdpix,/profile,win=win,$
border=border,xrange=xval,yrange=yval,/axes,sort=sort,$
_extra=e,nomodlut=nomodlut,clip=clip,$
usecongrid=usecongrid,noscale=noscale,nodisplay=nodisplay,$
giffile=giffile ,xstyle=xstyle,ystyle=ystyle,chkinf=chkinf
double=double,imgoff=imgoff
ARGS :
d[m,n] : data to display
KEYWORDS:
zx: int x zoom factor
zy: int y zoom factor
clip[2]: float if supplied then clip the data values
before scaling to byte. Note for zx,zy<-1
this is done after rebinning (averaging).
nsigclp[2]: float if supplied then compute the clipping levels using a
robust estimator of rms (median, throw out outliers, repeat).
Data will be clipped to [-nsigClip[0],nsigClip[1]]*sig after removal
of mean. This overrides clip keyword.
if nsigclip is a single element then clip to
[-sig,sig]*nsigclip.
if nsigclip[0] is negative, force it positive.
maxiter:int If nsigclip is used, the maximum default interation
is 5. the maxiter keyword can change this.
usecongrid: if set and zx or zy, then use congrid rather than
rebin for the new image.
histeq: set for histogram equalization of image disp
sort: if set then use sort for histogram equalization.
rdpix: set to turn on read pixel function.
profile: set to turn on profiles command
win: int window number to use. default:1
poswin[2]:int x,y position for lower left corner of window
(in screen coord)
border: int number of pixels around image to allow for labels.
(ignored if ps, or !p.multi in use).
xrange[2]:float xrange for labeling the x axis (min,max)
yrange[2]:float yrange for labeling the y axis (min,max)
axes: if set then draw the axis and labels
invert: if set then invert the colors white and black
noscale: if set then do not bother to scale image.
Use this when you are redisplaying an image that
was returned by imgdisp.
nodisplay: if set then just return the scaled image, do not display
giffile: string if present then write a gif image to this file
xstyle: int if set then use for xstyle= keyword to plot. def:1
ystyle: int if set then use for ystyle= keyword to plot. def:1
chkinf: int if set then check for infinities and exclude during
nsigclip
double: int if set then use doubles when computing nsigclip
_extra=e passed to plot routine when drawing the axes.
can include, xtitle,ytitle,title, etc..
RETURNS:
ret[] : array actually displayed after scaling.(optional)
imgoff[2] : the x,y offset for img in window in normalized device coord.
use this to replot the image in the same window:
rv ret,imgoff[0],imgoff[1],/norm
DESCRIPTION:
Display an m by n array as an image. The data is passed in
via the array d. It will normally create the image in window 1. This can
be overridden with the win= keyword.
zx,zy allow you to zoom the image in the x or y dimension before display.
They must be integral values (rebin is used). gt 1 make the dimension
larger while lt -1 make the dimension smaller by this amount. If the
zoom is lt 0 then it must divide evenly into the dimension of the array.
The normal scaling of the data values is:
(d-minVal)/(maxVal-minVal) * number of ofColors in lookup table.
You can switch to histogram equalization by setting /histeq. If your data
has outliers then the /sort option will make the histogram equalization
work better (but it takes a little longer).
If the keyword axes is set, then tickmarks and labels will be drawn
around the image. xrange,yrange hold the min,max values for each
axis. If you do not specify xrange or yrange then the values 0..ndim-1
is used. You can also pass in xtitle=,ytitle=,title= to label the
axis and the plot.
Setting /rdpix will allow you to readback the position and
pixel values interactively.
You can change or load the lookup table with xloadct.
Multiple plots per page.
The routine can also be used to place multiple images per page. In
this case you must set !p.multi=[plotsLeft,ncols,nrows] outside of
this routine. You must also define the size of the window you want
to use before the first call. The border keyword is ignored.
Each image will be scaled (congrid) to fit into the plot window.
Postscript output.
You can also generate postscript files. To do this you call
psimage,args before calling imgdisp. imgdisp will check
!d.flags for scalable pixels to determine if it is writing to a postscipt
file. When postscript output is enabled the border keyword is ignored.
If !p.multi is not set then the image will fill the drawing area
(set in the call to psimage...default 7x9in) while keeping the aspect
ratio of x to y . If !p.multi is set then each image is scaled to the
plot window with congrid. When done with the postscript output you must
call hardcopy to flush the buffer.
The postscript output will look at the current lookup table. If it
is not 0-255 then it will scale the data to 0-255 and then do
an indirect lookup through the lookup table. This allows you to display
the image on the screen, change the lookup table with xloadct and then
have the changes appear on the hardcopy output. On exit the
original lookup table is restored.
Gif output
A gif file will be created if you set the giffile keyword to a file
name (include the .gif if you want it there). The routine will read the
current screen and output the file. The current lkup table is used.
For large images use border=70 so the labels fit on the page.
Examples:
assume the images are d[1024,180]
1. display image on the screen, label with a border of 50 pixels:
imgdisp,d,border=50,/axes,xra=[1400.,1500.],xtitle='freq [Mhz]',$
title='22jan01 spectra pol A'
2. above, but display in landscape mode to a postscipt file:
imgdisp,d,border=50,/axes,xra=[1400.,1500.],xtitle='freq [Mhz]',$
title='22jan01 spectra pol A'
the border keyword is ignored.
3. place 6 images per page. assume d is [1024,180,6].
window,1,xsize=1000,ysize=600 ; plot will fit in here
for i=0,5 do begin
!p.multi=[(6-i mod 6),2,3]
title=string(format='("image looking at distomat ",i2)',i+1)
imgdisp,d[*,*,i],/axes,xra=[1400.,1500.],xtitle='freq [Mhz]',$
title=title
endfor
For the same thing with postscript output:
psimage
for i=0,5 do begin
!p.multi=[(6-i mod 6),2,3]
title=string(format='("image looking at distomat ",i2)',i+1)
imgdisp,d[*,*,i],/axes,xra=[1400.,1500.],xtitle='freq [Mhz]',$
title=title
endfor
hardcopy
x ; to get back to x windows.
GOTCHAS:
1. When making images it is important to get all 256 colors for your
lookup table. The only way to guarantee this in idl is:
idl
p8 .. set to pseudo color mode
window,colors=256
Then proceed with normal processing.
2. When scaling down a plot or when putting multiple plots per page
it is easy to loose information. If you have placed horizontal or
vertical lines in the plot for reference (and they are only 1 line
wide) they may not appear on the final image.
3. landscape ps output does not always appear where you think it should. Any
page offsets are first applied and then the image is rotated by 270
degrees so that the xoffset points up and the yoffset points to the left.
psimage,/lanscape will fudge the offsets so that
x'=yoff,y'=maxLen-xoff. so the origin is at the upper left.
rotation by 270. this leaves:
y''=x, x''=(maxlen-1) so the image is inverted. It would have been nice
if idl did a rotate by 90 and then flip about the vertical axis...
For those that don't want psimage doing this magic, set xoff=0,yoff=0
with /landscape and this extra addition won't happen (but your plot
won't be visible unless you play some games..).
To get a good plot out in landscape mode try using the pstops command
outside of idl to straighten things..
psimage,/landscape
plot the image
hardcopy
x
cmd=strarr(4)
cmd[0]='/pkg/image/bin/pstops'
cmd[1]='1:0u(8.5in,11in)'
cmd[2]='/dir.../idl.ps'
cmd[3]='/dir.../idlfixed.ps'
spawn,cmd,reply,/noshell
This should put the image in the correct location on the page.
You need to replace dir... with your directory..
SEE ALSO:
imgflat,imgflaty,imghisteq,corimgdisp,corimgonoff
(See /pkg/rsi/local/libao/phil/gen/imgdisp.pro)
NAME:
imgflat - flatten an image.
SYNTAX: result=imgflat(data,code,ravg=ravg,col=[x1,x2],sig=sig,median=median,
bptouse=bptouse,nobpc=nobpc,bpZero=bpZero)
ARGS:
data[x,y] data to operate on
code : int how to flatten image:
0 bandpass correct using the average of all strips.
if /med set then use the median rather than the average.
n bandpass correct n strips at a time, averaging n strips
if /ravg is set then this is a running avg. If not,then
break image into y/n sections and do 1 avg for each section.
-n bandpass correct using bandpasses +/- n strips from the
current line.
KEYWORDS:
ravg : if set and n> 0 then use a running average of n
strips about the current line.
col[2] : after bandpass correction average columns x1 thru x2
and divide this into every column (to remove things
like continuum sources. Note that x1,x2 are zero based..
sig : if set, return map in units of sigma
median : if set then use the median rather than the average for code
bptouse[x] : if provided, use this for band pass correction. ignore code
nobpc : if set then do not do a bandpass correction.
bpZero : if set then the data has a mean close to zero. Division
by an average bandpass. In this case:
bpavg=total(data,2)
val=mean(bpavg)
bpavg=(bavg - val) + bpZero)
to correct:
data[*,i]=data[*,i]/bpavg - val
noSub : if set then divide by bandpass correction but don't subtract
1. This is usually used with mean(data) = 0 and user
supplies bptouse=
RESTRICTIONS:
if n > 0 then y must be divisible by n
if chn provided then x1 le x2 le x
DESCRIPTION:
The data array d[x,y] has x xpoints by y ypoints.
The processing steps are:
1. bandpass correct using the code provided.
2. if chn[] is specified then average columns x1 through x2 and divide this
into every column.
3. subtract 1 from the map.
4. if sig is set, compute and then divide the map by the maps sigma.
EXAMPLE:
1. assume we have correlator data of 1024 lags by 300 records. Display
sbc 1 pol A and use columns 800-900 for leveling. On display, scale
the data to 2% of Tsys:
img=imgflat(b.b1.d[*,0],0,col=[800,900])
imgdisp,(img > (-.02))<.02,zy=2
Note: corimgdisp() does all this for you in 1 routine.
2. Assume position switch correlator data of 300 recs/scan. Display
all 600 records using the 300 off records for the bandpass correction.
scale to 5% of Tsys:
print,corgetm(lun,600,b,scan=scan) ; read in 600 recs starting at on.
sbc=0
pol=0
bpc=coravg(b[300:599].(sbc).d[*,pol]) ; average the off data
img=imgflat(b.(sbc).d[*,pol],0,bptouse=bpc)
imgdisp,(img > (-.05))<.05
(See /pkg/rsi/local/libao/phil/gen/imgflat.pro)
NAME:
imgflaty - flatten an image in the y direction
SYNTAX: result=imgflaty(data,x1,x2)
ARGS:
data [m,n] data to operate on
x1 int index col average start (count from 0)
x2 int index col average end (count from 0)
DESCRIPTION:
The data array d[m,n] has m xpoints by n ypoints.
average columns located at x1 thru x2 to give a[n].
expand a to be a[m,n] by duplicating the columns
retun data/a
(See /pkg/rsi/local/libao/phil/gen/imgflaty.pro)
NAME:
imghisteq - histogram equalize an image. return byte array
SYNTAX: bytarr=imghisteq(data,stretch=stretch,invert=invert,minv=minv,
maxv=maxv,sort=sort)
ARGS:
data[m,n] : data array to equalize
KEYWORDS:
stretch[4]: after histogram equalization,
map data range s[0]-s[1] (0 to 255)
to data range s[2]-s[3] (0 to 255)
minv : float .. minv to use for histeq_no stretchclipping
maxv : float .. maxv to use for histeq_no stretchclipping
sort : if set then histeq via sort
(See /pkg/rsi/local/libao/phil/gen/imghisteq.pro)
NAME:
imghline - draw horizontal line on an image
SYNTAX: imghline,img,linind,dashlen,vlines,val,ononly=ononly
ARGS:
img[n,m] : float image to display
linind[k]: int vertical indices into img array for horizontal lines
(count from 0)
dashlen : int number of pixels for on dash.
vlines : int number of vertical lines for each dash. def:1
val : float value to use for dash.
DESCRIPTION:
Draw horizontal dashed lines into a 2d image array. By default the
dashes are 4 pixels on, 4 pixels off. You can change the length with
dashlen. The only requirement is that dashlen must divide into the
first dimension of the img array. The dashes will be 1 horizontal line
thick. You can make them wider with the vlines keyword. This is sometimes
necessary if the image is being scaled down on display).
(See /pkg/rsi/local/libao/phil/gen/imghline.pro)
NAME:
intermods - compute intermods between 2 freq.
SYNTAX: n=intermods(f1,f2,minfreq,maxfreq,maxorder,outI,nf1,nf2,
neg=neg,all=all)G
ARGS:
f1: float Mhz. First frequency to use
f2: float Mhz. 2nd freq to use.
minFreq: float Mhz. Minimum intermod freq to keep
maxFreq: float Mhz. Maximum intermod freq to keep
maxOrder: int maximum intermod to compute on each freq.
eg 5 --> up to f1^5, f2^5
KEYWORDS:
neg: If set then include intemods < 0.
all: if set then include all intermods found.
by default only unique intermods are retured even if there
different orders that generate them.
RETURNS:
n : int number of intermods we found with requested range
outI[n]: float the intermod value for each found
nf1[n] : int the order for f1 for the nth intermod
nf2[n] : int the order for f2 for the nth intermod
DESCPRIPTION:
Compute intermods between the two frequencies input (f1,f2). If only
one frequency is entered (f1) then just compute the harmonics of this
frequency.
The maxorder keyword tells how high an order to use. All of the
intermods of the various harmonics of f1,f2 (up to maxorder) are computed.
Those they lay within the range minFreq,maxFreq are kept.
The intermod differences are returned in outI. The order for each intermod
is returned in the array nf1 and nf2.
(See /pkg/rsi/local/libao/phil/gen/intermods.pro)
NAME:
intm_pdevalfa - compute alfa/pdev mixer intermods
SYNTAX:intm_pdevalfa,skycfr,rdrFrqAr,retI,code=code,lo2=lo2,$
maxorder=maxorder,minmaxfrq=minmaxfrq,all=all,$
ARGS:
skycfr: float Mhz sky cfr for alfa band. The first lo will be
250 Mhz above this value.
rdrFrqAr[n]: float Mhz freq of each radar to check.
KEYWORDS:
code: int 0 default compute intermods for the 1st mixer
1 compute intermods for the 2nd mixer
using the lower pdev band (lo1=175 1450 band)
2 compute intermods for the 2nd mixer using the
2nd pdev band (lo1=325 1300 band)
lo2: float Mhz if supplied then this is the lo2 to use.
It will override the code 1,2 selection.
maxorder: int maximum harmonic order to contemplate. Default is
10
minmaxfrq[2]: float Mhz. min,maximum frequency output from the mixer
to keep.
The default values are:
code
0 100 to 400 since:
the filters in if1 are 0..250 and 0-400.
The bbandFilters are +/= 75 Mhz.
1 -/+ 86 ... since thats the output band
2 -/+ 86 ... since thats the output band
html : if set then output html code to make a table.
RETURNS:
retI[n] : {} array of structures holding intermod info
DESCRIPTION:
Compute mixer intermods for alfa and pdev. The user specifies:
1. skycfr sky center frequency used for the alfa band. Normally this is
lo1-250.
2. rdrFrqAr[n] - array of radar frequencies to test for intermods
By default it computes the intermods in the first mixer. The keyword
code=1,2 will compute intermods in the 2nd mixers (1=175,2=325).
EXAMPLE:
skycfr=1375
rdrFrqAr=[1330,1350.]
intm_pdevalfa,skycfr,rdrFrqAr,retI
; use intm_pdevalfa_pr to print out the results
intm_pdevalfa_pr,retI
(See /pkg/rsi/local/libao/phil/gen/intm_pdevalfa.pro)
NAME:
intm_pdevalfa_pr - print out intm_pdevalfa intermod info
SYNTAX:intm_pdevalfa_pr,retI,html=html
ARGS:
retI[n]: {} structure from intm_pdevalfa to print info.
KEYWORDS:
html : if set then output html code to make a table.
RETURNS:
DESCRIPTION:
Print out the info returned from intm_pdevalfa.
(See /pkg/rsi/local/libao/phil/gen/intm_pdevalfa_pr.pro)
NAME:
inverf - compute inverse error function
SYNTAX: val=inverf(x)
ARGS:
x[n] : float/double evalute the function at x.
RETURNS:
val[n] : double the inverse error function value.
DESCRIPTION:
Compute the inverse error function.
The following approximations to the inverse of the error function are
taken from J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational
Chebyshev Approximations for the Inverse of the Error Function",
Mathematics of Computation, 30 (1976) 827-830 + microfiche appendix.
via fred schwab
(See /pkg/rsi/local/libao/phil/gen/inverf.pro)
NAME: ipsinit - initialize to use the idl ips routine for mock data SYNTAX: @ipsinit DESCRIPTION: call this routine before using any of the ips idl routines. It sets up the path for the idl ips directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/ipsinit.pro)
NAME:
isleapyear - check if year is a leap year.
SYNTAX: istat=isleapyear(year)
ARGS: year[]: int/long 4 digit year
Returns:
istat[]: int 0 if not a leap year.
1 if a leap year.
DESCRIPTION:
Determine whether a year is a leap year in the gregorian calendar.
Leap years are those years
divisible by 4 and (!(divisible by 100) or (divisible by 400)).
eg. (1900 is not a leap year, 2000 is).
The input can be a scalar or an array.
(See /pkg/rsi/local/libao/phil/gen/isleapyear.pro)
NAME:
lbgain - compute lband gain as a function of az,za
SYNTAX: gain=lbgain(az,za)
ARGS : az[n] - float azimuth in degrees.
za[n] - float zenith angle in degrees.
RETURNS:
gain[n] - float kelvins/Jy.
DESCRIPTION
Compute the gain of the lband system in Kelvins per Jansky. Data
was taken from jul00,aug00 lbn and lbwide on/off position switching
data using the correlator (see http://www.naic.edu/~phil). 1405 Mhz
was used and polA polB were averaged together.
(See /pkg/rsi/local/libao/phil/gen/lbgain.pro)
NAME:
ldcolph - load phil's colors into the lookup table
SYNTAX: ldcolph,pscol=pscol,maxlen=maxlen
ARGS:
KEYWORDS:
pscol : if set then use colors for printing (background
is white, not black
maxlen: int maximum length to load.
default: is entire lut (256 entries) with
last value extended to end of lut
DESCRIPTION:
The common block colph is loaded with 24 bit true color values.
By default 11 colors are loaded. If maxlen is < 11 then load only
the first maxlen values. The
If a function includes the colph common block, then the user
case use these values with the color= keyword in the plot routine.
The common block is defined as:
common colph,decomposedph,colph
The colors are:
screen printing (if different)
0 - black white background
1 - white black foreground
2 - red
3 - green
4 - blue
5 - yellow redish (bright yellow on white does not print well)
6 - pink
7 - aquamarine
8 - grey
9 - light brown
10 -purple
EXAMPLE:
include the common block in any function. By default @geninit will
include it in the top level.
plot,x,y,col=colph[3] this will be a green line.
(See /pkg/rsi/local/libao/phil/gen/ldcolph.pro)
NAME: lutcycle - cycle through all the idl luts.. SYNTAX: lutcycle,delay ARGS: delay: int/float. secs to wait between each step. DESCRIPTION: Cycle through all of the lookup tables supplied by idl. EXAMPLE: display an image. imgdisp,dat lutcycle,5 ; cycle through 40 luts waiting 5 seconds at each 1.
(See /pkg/rsi/local/libao/phil/gen/lutcycle.pro)
NAME: masinit - initialize to use the mock spectrometer fits routines SYNTAX: @masinit DESCRIPTION: call this routine before using any of the mas... idl routines. It sets up the path for the idl mas directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/masinit.pro)
NAME:
maskbyrms - create mask using rms of fit residuals.
SYNTAX: mask=maskbyrms(x,y,deg=deg,maxloop=maxloop,nsig=nsig,$
indxgood=indxgood,indxbad=indxbad,$
nbad=nbad,ngood=ngood,verb=verb,coef=coef)
ARGS:
x[n] : float x data
y[n] : float y data
RETURNS:
mask[n]:long holds the mask. 1 for good points, 0 for bad points
indxgood:long indices of y whose mask value is 1
indxbad:long indices of y whose mask value is 0
ngood:long the number of good (mask=1) points.
nbad:long the number of bad (mask=0) points.
coef[deg+1]: float the coef from poly_fit
KEYWORDS:
deg: int degree for polynomial fit (def:1)
maxloop: int max time to loop on fit before quitting. (def:20)
nsig: float keep points within nsig of fit (def: 3. sigma)
verb: if set then plot out the residuals and mask after the
fitting.
DESCTRIPTION:
This routine will fit a polynomial of order deg to the data x,y.
It will then remove all points whose fit residuals are greater then
nsig sigmas and iterate the fitting process. When all residuals are within
nsig, a mask is created where the remaining points have a 1 and the points
that were excluded have a mask value of 0.
NOTE: this routine uses the value of x that is provided. For large values
of deg you should scale x so that x**deg does not overflow (say -1,1).
SEE ALSO:
blmask, cormask, corblauto
(See /pkg/rsi/local/libao/phil/gen/maskbyrms.pro)
NAME: matrot - generate 1 or more rotation matrices SYNTAX: mat=matrot(order,th1,th2,th3) ARGS: order: string "xyz" order for the rotations th1 : float first angle in radians th2 : float 2nd angle in radians th3 : float 3rd angle in radians RETURNS: mat[3,3]: float rotation matrix DESCRIPTION: concatenate 3 rotation matrices. The matrices will be applied using the order string order (eg. if order='zxy' then a z rotation,followed by an x rotation followed by a rotation about the y axis. The direction of the rotation: +theta about an axis generates a CCW rotation of a vector (Note this direction was updated on 17jun20). Order of data EXAMPLE: rotate 90deg about z then x then y use vectors [1,0,0] and [0,1,0] axis='zxy' V=fltarr(3,2) V[*,0]=[1,0,0] V[*,1]=[0,1,0] thrRd=90.*!dtor M=matrot(axis,'y',thrd,thrd,thrd) to apply the matrix to a vector use: v1=transpose(m)#v of matrix_multiply(m,v,/atranspose) print,v1 1.00000 -4.37114e-08 0.00000 0.00000 1.91069e-15 1.00000
(See /pkg/rsi/local/libao/phil/gen/matrot.pro)
NAME:
mav - multiply an array by a vector
SYNTAX: val=mav(a,v,sec=sec)
ARGS:
a[n,m] : array
v[n] : vector
KEYWORDS:
sec: if set then v should match the 2nd dimension of a
returns:
val[n,m]
DESCRIPTION:
return val[i,j]= a[i,j]*v[j].. i=0,n-1,j=0,m-1
the routine will make v'[n,m] where v[i,*] is the same value
(See /pkg/rsi/local/libao/phil/gen/mav.pro)
NAME:
maxmin3d - find max, min for an array of 3d points
SYNTAX: mmAr=maxmin3d(pnts)
ARGS:
pnts[3,n]:double points to search
RETURNS:
mmAr[3,2]: double holds,max,min for each dimenstion
[*,0] are the mins, [*,1] are the max'es
DESCRIPTION
return max,min info for the 3d array of points.
EXAMPLE:
let pnts[3,n] be a double array of x,y,z points
mmAr=maxmin3d(pnts)
mmAr[*,0] = min x,y,z values
mmAr[*,1] = max x,y,z values
note that the max or min values do not come from 1 point. they are the max
min in each x,y,z dimension
(See /pkg/rsi/local/libao/phil/gen/maxmin3d.pro)
NAME:
mcalinp - input data for meascal routine.
SYNTAX: d=mcalinp(lun,type,numsteps,numloops,scan=scan)
ARGS :
lun : int of correlator file to read
type: int 1- on absorber, 2 on sky
numsteps: int number of steps to cover the frequency range
these are usually the 100 Mhz junks.
numloops: int the number of types you swept the entire frequency
range. The scans must be contiguous on disc.
KEYWORDS:
scan :long scan number to position to before reading. If
not provided then read from the current position.
not provided then read from the current position.
maskSig :float sigmas to use for rms by record masking
RETURNS:
d[numsteps*numloops*4]: {meascal} return the data here (see below).
DESCRPTION:
The meascal data acquisition routine steps through a frequency range
turning the cal off then on. It is normally run on absorber and then
on the sky. The setup is 4 sbc by 256 lags spanning 4*25 Mhz in a chunk.
An option of the routine is to loop multiple times through the frequency
range.
This routine reads the data and separates each sbc into a separate
element in the {meascal} array. The user provides the lun, type (
1 for absorber, 2 for sky) and the number of 100 Mhz steps and times
the entire frequency range was repeated.
The returned data array contains and entry for each sbc of each step:
d.frq - center freq of sbc in Mhz
d.type - 1 for absorber, 2 for sky (user supplies this value).
d.scan - scan number for the cal off scan
d.brd - board number in correlator 0..3
d.spOn[256,2] -holds the calon spectra for polA and polB
d.spOff[256,2] -holds the caloff spectra for polA and polB
d.spCal[256,2] -holds the calOn/caloff -1 spectra for polA and polB
d.tpOn[2] - total power cal on (pola,polb)
d.tpOff[2] - total power cal off (pola,polb)
(See /pkg/rsi/local/libao/phil/gen/mcalinp.pro)
NAME:
meanrob - robust mean for 1d array
SYNTAX: mean=meanrob(y,nsig=nsig,double=double,sig=sig,$
gindx=gindx,ngood=ngood,bindx=bindx,nbad=nbad,$
fpnts=fpnts,iter=iter,maxiter=maxiter,chkinf=chkinf)
ARGS:
y[n] : array to compute robust mean
KEYWORDS:
nsig : float use nsig*sigma as the threshold for the points to
keep on each iteration. The default is 3.
double: if set then force computation to be done in double
precision.
RETURNS:
mean: float/double the computed mean
sig float/double the last computed rms
fpnts : float the fraction of points used for the final
computation
gindx: long[] indices into d for the points that were used
for the computation.
ngood long number of points in gindx.
bindx: long[] indices into d for the points that were not used
nbad long number of points in bindx.
iter long number of iterations performed.
maxiter long maximum number of times to iterate. default=5
chkinf int if set then check for infinities and ignore these
points
DESCTRIPTION:
compute the robust mean for the input data array. The program loops
doing:
0. create a mask that includes all the points.
1. compute the mean, rms over the current mask
2. Find all points in the original array that are within nsig*sig of
the mean. This becomes the new mask. If the new mask has fewer
points than the old mask, go to 1.
4. Return the last mean computed. If the keywords are present, return
the sig, index for good points, index for bad points, and the fraction
of points used in the final computation.
.
(See /pkg/rsi/local/libao/phil/gen/meanrob.pro)
NAME:
meanrun - compute the running mean of a 1 or 2d array
SYNTAX: result=meanrun(data,len)
ARGS:
data[x,y] data to operate on
len : int length of running mean. If even , round up to next
odd number
DESCRIPTION:
Compute a running mean along the last dimension of the array data.
data can be a 1 or 2d array; The output will be float (unless data is
a double array in which case it will be double). For each point
data[j,i] the routine will average:
data[j,(i-len/2):(i+len/2)] points.
The edges will only average the above points which fall within
the index bounds of the array.
(eg data[j,0:len/2])
(See /pkg/rsi/local/libao/phil/gen/meanrun.pro)
NAME:
medianbychan - median 2d array by chan.
SYNTAX : result=medianbychan(d,nsections=nsections,retsection=retsection)
ARGS :
d[m,n] : input array to compute median over 2nd dimension (n).
KEYWORDS:
nsections: int if provided, then break d[*,n] up into nsections sections. Compute the
median of each of these. For each channel return the minimum
of the nsections values for each channel.
retsection: int If nsections is used then the minimum value of the nsections
measurements is returned by default. retsection lets you change
how the return value is determined.
0: return the minimum value (this is the default).
1: return the maximum value
2: return the average value
3: return the median value
RETURNS:
result[m]: result[i]= median(d[i,*])
DESCTRIPTION:
Compute the median by channel for a 2d array. If the array
is dimensioned d[m,n] then result[i]=median(d[i,*]).
The nsections keyword can be used to break the n samples up into nsections
units and compute the median of each of these separately. The minimum value
of the nsections samples is then returned for each sample. The retsection keyword
lets you changes this to the avg,max,or median.
(See /pkg/rsi/local/libao/phil/gen/medianbychan.pro)
NAME: mkallidldoc - create all html documentation. SYNTAX: @mkallidldoc DESCRIPTION: Create all of the html documentation in the directory specified by aodefdir(/doc). The routine will create a temporary file /tmp/mkallidldoc.pro and then executes it. It deletes it when done. You need write access to the aodefdir(/doc) directory and to /tmp
(See /pkg/rsi/local/libao/phil/gen/mkallidldoc.pro)
NAME:
mkazzagrid - make a grid of az,za values.
SYNTAX: mkazzagrid,az,za,azstart=azstart,azend=azend,azstep=azstep
zastart=zastart,zaend=zaend,zastep=zastep
RETURNS:
az[nptsaz,nptsza]: float return azimuth values here
za[nptsaz,nptsza]: float return za values here
KEYWORDS:
azstart : float. starting azimuth . default 0.
azend : float. ending azimuth . default 359.
azstep : float. step size for azimuth points. def:1
zastart : float. starting za. . default 0.
zaend : float. ending za. . default 20.
zastep : float. step size for za points. def:.5
DESCRIPTION:
Return the 2d arrays az,za filled in with the requested az,za.
These values can than be used to evaluate 2d functions of (az,za)
or for plotting 2d fields.
(See /pkg/rsi/local/libao/phil/gen/mkazzagrid.pro)
NAME:
mksin - make a sine wave
SYNTAX: d=mksin(len,numcycles,phase=phase,double=double)
ARGS :
len : long.. number of points
numcycles:float.. number of cycles in len
KEYWORDS:
phase : float .. starting phase in fraction of a cycle
double: if set then return double precision values
RETURNS:
d[len] :float .. the sinwave
(See /pkg/rsi/local/libao/phil/gen/mksin.pro)
NAME:
MK_HTML_HELP_PH
PURPOSE:
Given a list of IDL procedure files (.PRO), VMS text library
files (.TLB), or directories that contain such files, this procedure
generates a file in the HTML format that contains the documentation
for those routines that contain a DOC_LIBRARY style documentation
template. The output file is compatible with World Wide Web browsers.
This is a hack to mk_html_help for phil's documentation.
CATEGORY:
Help, documentation.
CALLING SEQUENCE:
MK_HTML_HELP_PH, Sources, Outfile
INPUTS:
Sources: A string or string array containing the name(s) of the
.pro or .tlb files (or the names of directories containing
such files) for which help is desired. If a source file is
a VMS text library, it must include the .TLB file extension.
If a source file is an IDL procedure, it must include the .PRO
file extension. All other source files are assumed to be
directories.
Outfile: The name of the output file which will be generated.
KEYWORDS:
TITLE: If present, a string which supplies the name that
should appear as the Document Title for the help.
VERBOSE: Normally, MK_HTML_HELP does its work silently.
Setting this keyword to a non-zero value causes the procedure
to issue informational messages that indicate what it
is currently doing. !QUIET must be 0 for these messages
to appear.
STRICT: If this keyword is set to a non-zero value, MK_HTML_HELP will
adhere strictly to the HTML format by scanning the
the document headers for characters that are reserved in
HTML (<,>,&,"). These are then converted to the appropriate
HTML syntax in the output file. By default, this keyword
is set to zero (to allow for faster processing).
BGCOLOR: background color. format is "#RRGGbb" (hex numbs) or "xlibcol"
COMMON BLOCKS:
None.
SIDE EFFECTS:
A help file with the name given by the Outfile argument is
created.
RESTRICTIONS:
The following rules must be followed in formatting the .pro
files that are to be searched.
(a) The first line of the documentation block contains
only the characters ";+", starting in column 1.
(b) There must be a line which contains the string "NAME:",
which is immediately followed by a line containing the
name of the procedure or function being described in
that documentation block. If this NAME field is not
present, the name of the source file will be used.
(c) The last line of the documentation block contains
only the characters ";-", starting in column 1.
(d) Every other line in the documentation block contains
a ";" in column 1.
Note that a single .pro file can contain multiple procedures and/or
functions, each with their own documentation blocks. If it is desired
to have "invisible" routines in a file, i.e. routines which are only
for internal use and should not appear in the help file, simply leave
out the ";+" and ";-" lines in the documentation block for those
routines.
No reformatting of the documentation is done.
MODIFICATION HISTORY:
July 5, 1995, DD, RSI. Original version.
July 13, 1995, Mark Rivers, University of Chicago. Added support for
multiple source directories and multiple documentation
headers per .pro file.
July 17, 1995, DD, RSI. Added code to alphabetize the subjects;
At the end of each description block in the HTML file,
added a reference to the source .pro file.
July 18, 1995, DD, RSI. Added STRICT keyword to handle angle brackets.
July 19, 1995, DD, RSI. Updated STRICT to handle & and ".
Changed calling sequence to accept .pro filenames, .tlb
text librarie names, and/or directory names.
Added code to set default subject to name of file if NAME
field is not present in the doc header.
sep 19, 2002, pjp . got rid of mk_html_help at the beginning,
added bgcolor keyword
(See /pkg/rsi/local/libao/phil/gen/mk_html_help_ph.pro)
NAME: mm0ninit - initialize for the new mueller 0 processing SYNTAX: @corinit DESCRIPTION: call this routine before doing the new mueller 0 processing
(See /pkg/rsi/local/libao/phil/gen/mm0init.pro)
NAME: mm0ninitwas - initialize for the new mueller 0 processing SYNTAX: @corinit DESCRIPTION: call this routine before doing the new mueller 0 processing
(See /pkg/rsi/local/libao/phil/gen/mm0initwas.pro)
NAME:
monname - return month name given month number
SYNTAX: nm=monname(monnum)
ARGS :
monnum int month number 1 to 12
DESCRIPTION:
Return the 3 character name of the month given the month number
EXAMPLE:
monnum=3
monNam=monname(monnum) ; this returns 'mar'
(See /pkg/rsi/local/libao/phil/gen/monname.pro)
NAME:
montonum - convert ascii month to number 1-12
SYNTAX: num=montonum(month)
ARGS :
month string holding 3 character month
DESCRIPTION:
Given a 3 character month abreviation return the month of year (1..12).
(See /pkg/rsi/local/libao/phil/gen/montonum.pro)
NAME:
note - write a string at the requested line on the plot.
SYNTAX: note,line,lab,xp=xp,sym=sym,lnstyle=lnstyle,dyscale=dyscale,_extra=e
ARGS:
line: float linenumber to start on. 1-33 covers the plot.
lab : string to write on plot.
xp : float xposition to start at. 0 to 1 covers the plot. default is
center each line on page.
sym : int sym number. If present then plot the symbol at the start
of the line (leave some blanks in lab at the start).
lnstyle: int if provided then draw a short line of type lnstyle at the
beginning of the line (leave blanks at start of lab).
dyscale:float if set then scale the y spacing this amount.
DESCRIPTION:
Write a line of text on a plot. The default line position runs 1 through
33. Use the xp keyword to align the text horizonally. The sym and linestyle
keywords let you put lines, symbols at the start of your text so you can
define what they are.
If !p.multi is used for multiple pages then you must recompute where
the lines go. The line number is relative to the entire page, not the
current window of !p.multi.
(See /pkg/rsi/local/libao/phil/gen/note.pro)
NAME: p50init - initialize to process p50 data SYNTAX: @p50init DESCRIPTION: call this routine before processing p50 data. It includes the pcld routines (point cloud routines).
(See /pkg/rsi/local/libao/phil/gen/p50init.pro)
NAME: p8 - set frame buffer to pseudo color. SYNTAX: p8 ARGS: none DESCRIPTION: Set the terminals frame buffer to pseudo color. This should be done before any plotting is done.
(See /pkg/rsi/local/libao/phil/gen/p8.pro)
NAME:
pagesize - set the postscript page size.
SYNTAX: pagesize,fullpage=fullpage,yoff=yoff
KEYWORDS:
fullpage: if set then set the pagesize to 7 by 10 inches.
default is 7 by 5 inches.
yoff : if supplied then move the plot from the default
position this many inches on the page.
DESCRIPTION:
Set the postscript output page size. This should only be called
when you are plotting to the postscript device (ps,psimage,pscol).
SEE ALSO:
ps, pscol, psimage
(See /pkg/rsi/local/libao/phil/gen/pagesize.pro)
NAME: pdevinit - initialize to use the pdev spectrometer SYNTAX: @pdevinit DESCRIPTION: call this routine before using any of the pdev idl routines. It sets up the path for the idl pdev directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/pdevinit.pro)
NAME: plasmaden - compute cold plasma density given plasma frequency SYNTAX: Ne=plasmaden(freqMhz) ARGS: freqMhz[]:float plasma freq in Mhz RETURNS: Ne[ ]: float in density/cc DESCRIPTION: Given the plasma freq in Mhz return the electron density. Just solves Wp^2=4*pi^2*n0/me
(See /pkg/rsi/local/libao/phil/gen/plasmaden.pro)
NAME:
pltazelusage_12m - plot the 2D az,el coverage.
SYNTAX: pltazelusage_12m,az,el,title=title,sym=sym,over=over,_extra=e,$
nodec=nodec
ARGS:
az[npts] : float azimuth positions (degrees).
el[npts] : float zenith angle positions (degrees).
KEYWORDS:
title: string label for top of plot
sym: int symbol to plot at each position.Default is *.
over: if set then overplot this data with what is there.
nodec: if set then don't plot lines of constant dec
_extra: extra keyword values to pass to plot and oplot routine.
eg (color=n).
RETURNS:
DESCRIPTION:
Plot the azimuth, el positions as a cartesian x,y plot. Elevation is the radius.
This routine can give an idea of how well a set of sources has covered the sky
(See /pkg/rsi/local/libao/phil/gen/pltazelusage_12m.pro)
NAME:
pltazzausage - plot the 2D az,za coverage.
SYNTAX: pltazzausage,az,za,title=title,sym=sym,over=over,dx=dx,_extra=e
ARGS:
az[npts] : float azimuth positions (degrees).
za[npts] : float zenith angle positions (degrees).
KEYWORDS:
title: string label for top of plot
sym: int symbol to plot at each position.Default is *.
over: if set then overplot this data with what is there.
dx: The step size in feet along the x,y axis. default is 10 feet.
_extra: extra keyword values to pass to plot and oplot routine.
eg (color=n).
RETURNS:
DESCRIPTION:
Plot the azimuth, za positions as a cartesian x,y plot. The axes are
feet from the center of the dish (projected onto z=0). This routine can
give an idea of how well a set of sources has covered the dish.
(See /pkg/rsi/local/libao/phil/gen/pltazzausage.pro)
NAME:
pltazzausage12m - plot the 2D az,za coverage. mods for 12m
SYNTAX: pltazzausage12m,az,za,title=title,sym=sym,over=over,dx=dx,$
zamax=zamax,_extra=e
ARGS:
az[npts] : float azimuth positions (degrees).
za[npts] : float zenith angle positions (degrees).
KEYWORDS:
title: string label for top of plot
sym: int symbol to plot at each position.Default is *.
over: if set then overplot this data with what is there.
dx: The step size in feet along the x,y axis. default is 10 feet.
zamax : float max za to plot def is 20
radiusza: if set then radius will be zamax. default will be 870foot radius
srcind[npts]; int hold src number for each source; 0... nsrc-1
will plot with color isrc+1
_extra: extra keyword values to pass to plot and oplot routine.
eg (color=n).
RETURNS:
DESCRIPTION:
Plot the azimuth, za positions as a cartesian x,y plot. The axes are
feet from the center of the dish (projected onto z=0). This routine can
give an idea of how well a set of sources has covered the dish.
(See /pkg/rsi/local/libao/phil/gen/pltazzausage12m.pro)
NAME:
pltbits - plot a timing diagram of the input data
SYNTAX: pltbits,x,y,bitmask,col=col,maxbits=maxbits,over=over,$
off=off,lab=lab,inc=inc,_extra=e,gaps=gaps,$
font=font,cs=cs,ncs=ncs,th=th
ARGS:
x[] : x values for the data
y[] : array holding the bits to plot
bitmask: long bits to extract and plot
KEYWORDS:
col[]: long colors for each bit to plot
maxbits : long max number of bits you want to plot. This is used to
compute the vertical positioning of the y axis. Use
this with off, if you want override the default
positioning of the traces.
over : if set then overplot from previous call
off : float . add to vertical position of each bit.Default is
.06
lab[] : string. labels for each bit.
gaps : float if provided, then an extra trace will be provided at the
top of the plot. Any points that have an x difference
>= gaps will have a transition from 1 to 0. This lets
you see where there is valid data.
font : int to use. -1 def, 1 truetype (for hardcopy)
cs : float charsize to use for axis. def=1.
ncs : float charsize to use for bit labels. def=1.
th : float line thickness. def=1. 2 doubles width
DESCRIPTION:
Suppose an int or long array holds status information that is
packed bit by bit. An example would be the vertex data that has
motor status bits encoded into a single int. This routine will
plot 1 or more of the bits versus the x axis..
EXAMPLE:
Suppose dat[100] has status info in bit0 and bit5. To plot them
versus input use:
x=finggen(100)
pltbits,x,dat,'21'x
Bit 0 will be plotted versus x with bit 5 plotted above it versus x.
(See /pkg/rsi/local/libao/phil/gen/pltbits.pro)
NAME:
pltbycol - plot values by color
SYNTAX: pltbycol,x,y,grpid,col=colar,xtitle=xtitle,ytitle=ytitle,
title=title,sym=sym,_extra=e,over=over
ARGS:
x[npts]: float xdata to plot
y[npts]: float ydata to plot
grpid[npts]: long unique number identifying a particular group. It
should run 0 through maxgrp-1. It is used to
generate the color.
KEYWORDS:
colar[ncol]: long lut values to use.
xtitle: string xlabel for plot
ytitle: string ylabel for plot
title: string title for plot
sym : int symbol to use for ploting.
over : int if set, then overplot this dataset on previous
_extra=e : will be passed to plot and oplot routine
DESCRIPTION:
Plot x,y points. Group points by color. grpind[npts] is used to identify
the points that have a common color. colar[ncol] is used for the color.
If there are more groups than colors, then the colors get reused modulo
ncolors.
EXAMPLE:
An example would be ploting the za Error for a given set of sources
by za. Suppose the array of structures src[npts] has the following elements:
src[i].name - source name
src[i].za - za for a measurement
src[i].zaErr - za error for the measurement
The unique srcnames are:
names=src[uniq(src[sort(src.name)].name)].name
You could generate the grpind[npts] array by:
nsrc=n_elements(names) ; number of unique names
grpind=lonarr(npnts) ; will hold srcid 0..nsrc-1
for i=0,nsrc-1 do begin
ind=where(src.name eq names[i],count)
if count gt 0 then grpind[ind]=i
endfor
You could then call pltbysrc with:
pltbysrc,src.za,src.zaErr,grpind,sym=1
The default color array is:
colar=findgen(10)+1
This has 10 unique colors (usually setup by ldcolph). Colors get reused
every 10 sources
The default symbol is *.
(See /pkg/rsi/local/libao/phil/gen/pltbycol.pro)
NAME:
pncodeinfo - initialize the code info for the gpsl2c med,long codes
SYNTAX: maxcodes=gpsl2ccodeinfo(gpsl2ccodeI,longcode=longcode)
KEYWORDS:
longcode: if set then return setupinfo for the gpsl2 cl codes:
767250 chips. By default return the gpsl2 cm code info:
10230 chips.
RETURNS:
maxcodes: long number of different codes we know of (satellites)
gpsl2ccodI[maxcodes]: {} codeinfo struct for each code:
DESCRIPTION:
Initialize the codeInfo structure for the gps l2c codes for the 37
assigned gps prn satellites
IDL> help,pncodeI,/st
** Structure CODEINFO, 4 tags, length=92, data length=92:
prn int 1 ; satellite number
NUM_REG LONG 1
LEN LONG 1
NUM_FDBACK LONG 1
FDBACK LONG Array[20]
galois int i ; 0==> fibonnaci, 1==> galois lfsr
startVal ulong 0 ; for galois
endVal ulong 0 ; for galois
The codeinfo struct gets passed to the routine that generates the
code.
Note that the fdback register positions are 1 based. When we use them
in the shiftregcmp routine we subtract one to get the 0 based index into
the shift register.
When calling shiftregcmp() on the long codes, the routine is slow since
idl doesn't have circular bit shifts,.
(about 20 secs..).
(See /pkg/rsi/local/libao/phil/gen/gpsl2ccodeinfo.pro)
NAME:
pncodeinfo - initialize the code info for the pncodes
SYNTAX: maxcodes=pncodeinfo(pncodeI)
RETURNS:
maxcodes: long that we know how to create
pncodeI[maxcodes]: {} codeinfo struct for each code:
DESCRIPTION:
Initialize the codeInfo structure for the pncodes that we know how to
create. These codes are the ones generated by the pncode generators
online.
Each array element contains:
IDL> help,pncodeI,/st
** Structure CODEINFO, 4 tags, length=92, data length=92:
NUM_REG LONG 1
LEN LONG 1
NUM_FDBACK LONG 1
FDBACK LONG Array[20]
galois int i ; 0==> fibonnaci, 1==> galois lfsr
startVal ulong 0 ; for galois
endVal ulong 0 ; for galois
The codeinfo struct gets passed to the routine that generates the
code.
Note that the fdback register positions are 1 based. When we use them
in the shiftregcmp routine we subtract one to get the 0 based index into
the shift register.
(See /pkg/rsi/local/libao/phil/gen/pncodeinfo.pro)
NAME:
pnthcoordsys - return coordinate system code
SYNTAX: icode=pnthcoordsys(pnthdr)
ARGS:
pnthdr:{hdrpnt} .. pnt portion of header.
RETURNS:
icode: int code for coordinate system used:
1 galactic
2 B1950
3 J2000
4 BEcliptic
5 JEcliptic
6 Ra/Dec of date (current
7 hour angle dec
8 az,za of source
9 az,za of feed (az is 180 deg from source if dome used)
10 az,za of feed with no pointing model included
EXAMPLE:
suppose we have a correlator data record:
print,corget(lun,b)
icode=pnthcoordsys(b.b1.h.pnt)
To extract the data manually (without this routine) use:
icode= ishft(b.b1.h.pnt.stat and '00078000'XL, -15)
(See /pkg/rsi/local/libao/phil/gen/pnthquery.pro)
NAME:
pnthgrmaster - return 1 if greg is master, 0 if ch master
SYNTAX: istat=pnthgrmaster(pnthdr)
ARGS:
pnthdr:{hdrpnt} .. pnt portion of header.
RETURNS:
istat: int 1 if greg master, 0 if ch master
EXAMPLE:
suppose we have a correlator data:
print,corget(lun,b)
istat=pnthgrmaster(b.b1.h.pnt)
(See /pkg/rsi/local/libao/phil/gen/pnthquery.pro)
NAME:
posscan - position to a scan/record on disc
SYNTAX: istat=posscan(lun,scan,rec,retstdhdr=retstdhdr,sl=sl)
ARGS:
lun: int .. logical unit assigned to open file
scan: long.. scan number 0--> whatever scan fits .. current or next
full scan number --> position to scan,
no rewinding allowed
rec: long grp number within scan.
0 or not included--> next record available
number --> record of current scan
keywords
retstdhdr: if valid variable, return standard header here..
(only if we positioned successfully)
skip : int .. skip this many scans forward. should use with
scan=0.
sl[]: {sl} returned from getsl routine. If provided
then routine will position directly to the scan requested.
The routine will not backup from the current position. If it finds
an increasing scan/rec number then it returns -1
returns: 1 positioned ok
0 not found
-1 found increasing scan number
-2 scan not in scanloc array
iook - 1 ok, 0 i/oerror, -1 bad headerid
(See /pkg/rsi/local/libao/phil/gen/posscan.pro)
NAME:
prfgainall- compute fractional gain do to pitch,roll,focus
SYNTAX: fracgain=prfgainall(az,za,rcvr,pitch,roll,focus,freq=freq,
foctouse=foctouse,rolltouse=rolltouse,pitchtouse=pitchtouse)
ARGS: az[npts] - float azimuth degrees
za[npts] - float za degrees
rcvr : string: lb,sb,cb,xb 1400,2380,5000,10000
pitch : optional arg. return pitch here deg.
roll : optional arg. return roll here deg
focus : optional arg. return focus here inches
KEYWORDS:
freq : Mhz. if provided, then use the specified rcvr function
and add on the relative difference from this
and the rcvr default freq.
foctouse[npts]: float use this as the focus instead of the model
rolltouse[npts]: float use this as the roll instead of the model
pitchtouse[npts]: float use this as the pitch instead of the model
DESCRIPTION:
Compute the fractional gain do to the pitch roll and focus. We use:
1. the pitch,roll,focus fits from feb00
2. aoant for pitch,roll
3. focus curve sband for the loss do to focus
4. asssume pitch roll error adds in quadrature
5. assume we multiply the focus loss by pitch,roll loss
before calling this routine do @prfinit
EXAMPLE:
make a plot of gain versus az,za
keep za above 2 degrees since za=0 not measured.
mkazzagrid,az,za,azstep=10,zastep=1,zastart=2
gain=prfgainall(az,za,'sb')
stripsxy,az,gain,0,0,/stepcol
(See /pkg/rsi/local/libao/phil/gen/prfgainall.pro)
NAME: printpath - print out the path variable SYNTAX: printpath ARGS: none DESCRIPTION: print out the path system variable one path per line. The order is the order they appear in the path variable.
(See /pkg/rsi/local/libao/phil/gen/printpath.pro)
NAME:
prwspc - compute the power spectrum of the input signal..
SYNTAX: dfrq=pwrspc(dtm)
ARGS:
dtm[npts]: real or complex input time series
RETURNS:
dfrq[npts]: real .. power spectrum squared magnitude of xform
DESCRIPTION:
Return abs(fft(dtm))^2
(See /pkg/rsi/local/libao/phil/gen/pwrspc.pro)
NAME:
ps - send plot output to postscipt file.
SYNTAX ps,filename,_extra=e,fullpage=fullpage
ARGS:
filename: string filename for outputfile. default is idl.ps
KEYWORDS:
fullpage: if set then set the pagesize 5 by 10 inches.
the default is 5 by 7.
_extra : e pass to device command.
DESCRIPTION:
Set plot output destination to a postscript file.
When done plotting use:
hardcopy
x
to return to terminal output.
SEE ALSO:
pagesize, pscol, hardcopy, x
(See /pkg/rsi/local/libao/phil/gen/ps.pro)
NAME:
pscol - send plot output to color postscipt file.
SYNTAX pscol,filename,_extra=e,fullpage=fullpage
ARGS:
filename: string filename for outputfile. default is idl.ps
KEYWORDS:
fullpage: if set then set the pagesize 5 by 10 inches.
the default is 5 by 7.
_extra : e pass to device command.
DESCRIPTION:
Set plot output destination to a color postscript file.
When done plotting use:
hardcopy
x
to return to terminal output.
SEE ALSO:
pagesize, ps, hardcopy, x
(See /pkg/rsi/local/libao/phil/gen/pscol.pro)
NAME:
psimage - prepare to send image output to a postscript file.
SYNTAX: psimage,filename,xlen=xlen,ylen=ylen,xoff=xoff,yoff=yoff,$
landscape=landscape,xroff=xroff,yroff=yroff
ARGS:
filename: string filename to write to. default is idl.ps
xlen : float number of inches for xdimension. default
is 7 inches.
ylen : float number of inches for ydimension. default
is 9 inches.
xoff : float offset in inches for the left edge of image.
the default is to center the plot.
yoff : float offset in inches for the bottom edge of image.
the default is to center the plot.
xroff : float relative offset x direction in inches. add to xoff.
yroff : float relative offset y direction in inches. add to yoff.
the default is to center the plot.
landscape: if set then plot in landscape mode.
DESCRIPTION:
Set the output plot device to a postscript file. Add keywords for
image display (8 bits per pixel). Try and center the imageon the
page. Landscape mode causes problems when offsets are used. See
imgdisp for how to get around it.
SEE ALSO:
ps,hardcopy,x,imgdisp
(See /pkg/rsi/local/libao/phil/gen/psimage.pro)
NAME: psrfinit - initialize to use the psrfits mock spectrometer routines SYNTAX: @psrfinit DESCRIPTION: call this routine before using any of the psrfits... idl routines. It sets up the path for the idl psrfits directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/psrfinit.pro)
NAME: pupfinit - initialize to use the pupfits mock spectrometer routines SYNTAX: @pupfinit DESCRIPTION: call this routine before using any of the pupfits... idl routines. It sets up the path for the idl pupfits directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/pupfinit.pro)
NAME: pupiinit - initialize to use the pupi routines SYNTAX: @pupfinit DESCRIPTION: call this routine before using any of the pupi... idl routines. It sets up the path for the idl pupi directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/pupiinit.pro)
NAME: puprinit - initialize to use the puppi raw file routines SYNTAX: @puprinit DESCRIPTION: call this routine before using any of the pupr (puppi_*.raw)... idl routines. It sets up the path for the idl pupr directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/puprinit.pro)
NAME:
pwrlawdist - generate a power law distribution.
SYNTAX: y=pwrlawdist(alpha,npts,minVal=minVal,maxVal=maxVal,
dohist=dohist,nbins=nbins,h=h,xh=xh)
ARGS:
alpha: float exponent for power law.
npts: long number of points to return
KEYWORDS:
minVal: float minimum value for the distribution. If alpha is less
than 0., then minval should be gt 0 to keep it from
blowing up. default:.01
maxVal: float maximum value for the distribution. It should be
greater tham minVal.Default:1.
dohist: if set, then compute histogram, make log,log plot
of histogram, and then do a linear fit to non zero values.
nbins: long number of bins if histogram requested. default is
npts*.005
RETURNS:
y[npts]: double the power law distributed data.
h[nbins]: long return histogram of y. (linear)
xh[nbins]: float return center of each bin.(linear)
DESCRIPTION:
Generate a random variable with a power law distribution:
y=r**(alpha). The data range will be: (minVal ge y le maxVal).
The routine will optionally compute the histogram of the data,
fit a line to a log,log version of the histogram, and then
plot the histogram and the fit. This lets you see how close to the
ideal power law you got. The fit can have trouble if the data range
does not make a good match to the binsize,number of bins.
The method blows up for alpha=-1. If alpha is within
1e-4 of -1, the routine uses -1+/-1e-4 and gives a warning..
See numerical recipes in C, section: 7.2 transformation methods, page 287.
for an description of how it's done.
EXAMPLE:
1. create a power law distribution r^1.5 ranging from 0 to 1.
Return,10000L points. Make the histogram with 100 points.
y=pwrlawdist(1.2,20000L,minval=.0,maxval=1.,/dohist,h=h,xh=xh,$
nbins=100)
..create a power law dist r^-2.2 with 10000L points.
plot out the histogram. have the values go 1 to 500.
return the histogram with it's binvalues.
kick up the number of bins to 1000 so the histogram fit works.
y=pwrlawdist(-2.2,10000L,minval=1.,maxval=500.,/dohist,h=h,xh=xh,$
nbins=1000)
(See /pkg/rsi/local/libao/phil/gen/pwrlawdist.pro)
NAME: rcvnumtonam - convert receiver number to receiver name. SYNTAX: stat=rcvnumtonam(rcvnum,rcvnam,/num) ARGS: rcvnum : int receiver number 1..16 KEYWORDS: num : if set then user inputs the rcvnam and we return the rcvnum RETURNS: rcvnam: string receiver name stat : int 1 value receiver num, 0 invalid receivernumber DESCRIPTION Convert a receiver number to receiver name. Return the receiver name in the string rcvnam. Return the status :1 ok, 0 no receiver with this num in stat.
(See /pkg/rsi/local/libao/phil/gen/rcvnumtonam.pro)
NAME:
rdcur - read cursor position multiple times
SYNTAX: n=rdcur(curI,double=double)
KEYWORDS:
double : if set then return as doubles. def: floats
RETURNS:
n: long number of points read
curI[2,n]: float/double x,y data positions for each cursor position requested.
DESCPRIPTION:
Read the cursor position multiple times. Return the number of points
read in n. Return the positions in the array curI[2,n]. The positions
are returned in data coordinates.
As you move the cursor around the cursor position will continuously
update on the terminal. When the left or center button is depressed, the
position will be stored in curI and a newline will be started on the
terminal (so you can see the value that was stored).
Mouse Button usage is:
Button Action
Depress left: record position, new terminal line
Depress center: record position, new terminal line
Depress right: exit routine.
The routine will return a maximum of 1000 points.
SEE ALSO:
cp
NOTE:
Points are only recorded on the downward press of the left or center
button.
(See /pkg/rsi/local/libao/phil/gen/rdcur.pro)
NAME:
rdcurdif - read cursor position difference multiple times
SYNTAX: n=rdcurdif(curI,/double)
RETURNS:
n: long number of points read
curI[2,3,n]: float x,y data positions for each set
curI[0,0,*] x start
curI[1,0,*] y start
curI[0,1,*] x end
curI[1,1,*] y end
curI[0,1,*] xend -xstart
curI[1,1,*] yend -ystart
DESCPRIPTION:
let the user pick two positions with the currsor and then compute the
x,y difference between the two points.
Return the positions and differences in the array curI[2,3,n]. The positions
are returned in data coordinates.
Mouse Button usage is:
Button Action
Depress left: record start position, new terminal line
Depress center: record end position, compute difference
Depress right: exit routine.
As you move the cursor around there are continual updates:
1. on the first point, the 1st position updates
2. on the 2nd opint the the 2nd cursor position and difference
will continuously update on the terminal.
The routine will return a maximum of 1000 points.
SEE ALSO:
rdcur,cp
NOTE:
Points are only recorded on the downward press of the left or center
button.
(See /pkg/rsi/local/libao/phil/gen/rdcurdif.pro)
NAME: rdevinit - initialize to use the pdev radar processor SYNTAX: @rdevinit DESCRIPTION: call this routine before using any of the rdev idl routines. It sets up the path for the idl rdev directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/rdevinit.pro)
NAME:
readasciifile - read an ascii file into strarr
SYNTAX: nlines=readasciifile(filename,inpLines,comment=comment,$
fast=fast,maxlines=maxlines)
ARGS:
filename: string the filename to read
KEYWORDS:
comment: string single character that is a comment
lines that start with this will be skipped
fast: int if set then try fast read, don't use shell,wc
maxLine:long if fast used, then this is the max number of lines
to expect in file (including comments). def:1000 lines..
if exceeded, will return an error. It can be less
RETURNS:
nlines: long number of lines read
-1 if file does not exist or read error
-2 if fast mode used and maxlines too small
inpLines[nlines]: strarr lines read
DESCRIPTION
Read an entire file into a string array. 1 line per string index.
Skip any lines that start with the comment character.
Return the string array and the number of lines read.
EXAMPLES:
filename='savfiles.dat'
nlines=readasciifile(filename,inplines,comment=';')
Timing:
For large numbers of files, /fast will run about 20 times
faster (since it doesn't make a shell call)
(See /pkg/rsi/local/libao/phil/gen/readasciifile.pro)
NAME:
recombfreq - compute recombination line freq for atoms
SYNTAX: freq=recombfreq(atom,linenum,linestep,alpha=alpha,beta=beta,gamma=gamma,
delta=delta,eps=eps,pertbl=pertbl)
ARGS:
atom : string H,he,C,N,O atom to compute. (case insensitive).
linenum[n]: long transition number (lower level)
linestep[n]: long 1=alpha,2=beta...etc
if linenum is an array and linestep is a single value
then use this value for all of linenum. If any keyword
alpha,beta.. is entered, then ignor linestep.
KEYWORDS:
alpha : if set then return alpha series (deltan=1) (ignore linestep)
beta : if set then return beta series (deltan=2) (ignore linestep)
gamma : if set then return gamma series (deltan=3) (ignore linestep)
delta : if set then return delta series (deltan=4) (ignore linestep)
eps : if set then return eps series (deltan=5) (ignore linestep)
pertbl[]:{pertbl} return period table (name,isotope%, amu for all
elements we support)
RETURNS:
freq[n] : double in Mhz.
DESCRIPTION:
compute the recombination line frequencies for an atom. Supported
atoms are: H,He4,C12,N14,O16. You can supply 1 or more linenumbers to do
at once as well as the step (1,2,3). Using any of the series keywords:
alpha,beta,gamma,... will override the linestep parameter.
Taken from "Tools of Radio Astronomy", Rohlfs&Wilson,2000. pg 334.
looks like value are good to about 1 khz for H89a at 9Ghz.
(See /pkg/rsi/local/libao/phil/gen/recombfreq.pro)
NAME:
recombsearch - search for recombination lines within a freq range.
SYNTAX: n=recombsearch(freqMin,freqMax,match,lineStepR=lineStepR,$
lineNumR=lineNumR,atoms=atoms
ARGS:
freqMin: float min frequency in Mhz to use.
freqMax: float max frequency in Mhz to use.
KEYWORDS:
lineStepR[2]: int min,max step for line transitions to use. The
Default is 1 (alfa) to 10 (??)
lineNumR[2]: int min,max line numbers to use. The default is 75 to
500.
atoms[]: string An array of atoms to search. The default is:
'H','He','C'. See recombfreq keyword pertbl for
a list of the names.
RETURNS:
n : long number of transitions found.
match[n]: {} an array of structures holding the recomb lines that
were found.
DESCRIPTION:
Look for all of the atomic transitions between frequencies freqMin and
freqmax. Return the number of transitions found and an array (match) holding
the information of each transform. The match structure contains:
** Structure <82c1a1c>, 4 tags, length=24, data length=24, refs=1:
ATOM STRING 'H' .. name of the atom
LINENUM LONG 157 .. line number
LINESTEP LONG 1 .. transition 1=alpha,2=beta..
FREQ FLOAT 1683.20 .. rest frequency of transition.
By default the routine searches the atoms H, He,C for line steps of 1..10,
and linenumbers 75 through 500.
EXAMPLES:
Suppose you have data between 1664 and 1686Mhz and you want to find
all transitions that satify: step size 1 thru 5 in atoms H,He, and C.
n=recombsearch(1664,1685,match,lineStepR=[1,5],atoms=['H','He','C'])
plot your data flagging these transitions
corplotrl,b,match
SEE ALSO: recombfreq,corplotrl
(See /pkg/rsi/local/libao/phil/gen/recombsearch.pro)
NAME:
resampleary- resample an array in the y direction
SYNTAX:nyout=resampleary(arIn,yIn,yStpOut,arOut,$
y1=y1,y2=y2,
cntyOut=cntyOut;
ARGS:
arIn[nxIn,nyIn] : float input array
yIn[nyIn] : float yvalue for arIn[*,nyIn]
yStpOut : float y step value for output array
KEYWORDS:
y1 : float first y value to use for output
def is yIn[0]
y2 : float last y value to use for output
def is yIn[nyIn-1]
RETURNS:
nyOut : long number of y entries i output:arOut[*,nout]
arOut[nxIn,nyOut] : float averaged output
yOut[nyOut] : float the y values for each arOut[*,i]
cntyOut[nyOut] : long number of arIn[*,i] averaged for
each arOut[*,i]
DESCRIPTION:
Given an array arIn[nx,ny] of unevenly sampled data in the
y direction, create an array of evenly spaced data in the
y direction arOut[nx,nyOut].
The yIn array has the y values for each arIn[*,iy].
Use the reverse indices of the histogram function on yIn
to find the arIn[*,iiy] elements that map to each arOut[*,iy]
yIn[] should be in increasing order. yStpOut defines the
fixed y spacing in arOut[*,iy]
cntYOut[nyOut] has the number of input rows arInp[*,i] averaged
for each arOut[*,j]
(See /pkg/rsi/local/libao/phil/gen/resampleary.pro)
NAME:
rfname - given the rfnumber, return the standard receiver name
SYNTAX : name=rfname(rfnum)
ARGS :
rfnum: int 1-16..
RETURNS:
name: string.. rfname.. sbn,sbw,lbn,lbw..etc..
illegal or unimplented rfnumber will be return as the string
-number.. eg rfnum 1 is not implented yet so it would return '-1'
(See /pkg/rsi/local/libao/phil/gen/rfname.pro)
NAME:
rms - compute the mean and standard deviation
SYNTAX: result=rms(x,quiet=quiet)
ARGS:
x[] : array to compute rms
KEYWORDS:
quiet: if set then don't print the rms,mean to stdout.
RETURNS:
result[2]: result[0]=mean, result[1]= std deviation
DESCTRIPTION:
compute the mean and standard deviation. Print the results to
stdio, and return in result[2]
(See /pkg/rsi/local/libao/phil/gen/rms.pro)
NAME:
rmsbychan - compute the rms/mean by chan for 2d array.
SYNTAX: result=rmsbychan(d,median=median,nodiv=nodiv,meanar=meanar)
ARGS:
d[m,n] : array to compute rms
KEYWORDS:
median: if set then use median rather than mean
nodiv : if set then don't divide by the mean
RETURNS:
result[m]: result[i]= rms(d[i,*])/mean(d[i,*])
meanar[m]: double mean(or median) for each channel
DESCTRIPTION:
compute the standard deviation/mean by channel.
(See /pkg/rsi/local/libao/phil/gen/rmsbychan.pro)
NAME:
robfit_poly - robust polyfit for 1d array
SYNTAX: coef=robfit_poly(x,y,deg,nsig=nsig,double=double,sig=sig,$
gindx=gindx,ngood=ngood,bindx=bindx,nbad=nbad,$
fpnts=fpnts,iter=iter,yfit=yfit,maxiter=maxiter,mask=mask)
ARGS:
x[n] : x array for fit
y[n] : y array for fit
deg : int deg of polynomial fit
KEYWORDS:
nsig : float use nsig*sigma as the threshold for the points to
keep on each iteration. The default is 3.
double: if set then force computation to be done in double
precision.
maxiter: maximum number of times to loop. default is 20.
mask[n] :intarr if supplied, then use this as a mask for starting
the fit. All zero elements of mask will be ignored
(on the first iteration).
RETURNS:
coef[deg+1]: float/double the fit coef.
sig float/double the last computed rms
fpnts : float the fraction of points used for the final
computation
gindx: long[] indices into d for the points that were used
for the computation.
ngood long number of points in gindx.
bindx: long[] indices into d for the points that were not used
nbad long number of points in bindx.
iter long number of iterations performed.
yfit : float/double the fit evaluated at x[n]
DESCTRIPTION:
compute a robust polynomial fit for the input data x,y. The program loops
doing:
0. create a mask that includes all the points.
1. fit the polynomial, rms residuals over the current mask
2. Find all points in the original array that are within nsig*sig of
the fit. This becomes the new mask.
3. if sig is less than the minimum sig so far, minsig=sig, and
store indices for this minimum sig
4 If the new mask does not have the same number of points than the
old mask, go to 1.
5. Check the coef of the Check if the
5. Return the last coefs computed. If the keywords are present, return
the sig, index for good points, index for bad points, and the fraction
of points used in the final computation, and yfit.
.
(See /pkg/rsi/local/libao/phil/gen/robfit_poly.pro)
NAME:
robfit_polyfft - automatic baselining of a function
SYNTAX: istat=robfit_polyfft(ydat,yfit,maskused,coef,$
edgefract=edgefract,masktouse=masktouse,$
nsig=nsig,ndel=ndel,maxloop=maxloop,$
deg=deg,sub=sub,verbose=verbose,fsin=fsin,
double=double,plver=plver
ARGS:
ydat[n] input data to baseline.
KEYWORDS :
edgefract[2]: float The fraction of each edge of the bandpass to not use
in the computations. if n=1 then use the same fraction for
both sides. If n=2 then use [0] for the left edge and [1] for
the right edge. The default value is .08.
The keyword masktouse will override this option.
masktouse:fltarr[n] if provided,then use this rather than edgefract for the
starting mask. This allows to you to force certain
areas to not be included in the fit. 0--> exclude
nsig: int The rms loop will throw out all outliers greater than
nsig*sigma. By default nsig is 3.
ndel: int The number of extra channels to remove whenever we
throw out a channel because of large rms. def=1
maxloop: int When removing outliers in the rms loop, it will continue
looping until no points are removed or maxloop iterations
is hit. By default maxloop is 15.
deg: int maximum degree of polynomial fit.
The routine starts at order 1 and iterates the
rmsloop,fit iterating the degree after each pass
until the last fitorder done is deg. By default deg is 1.
sub : if set then return the data-blfit rather than the fit.
verbose : 1 plot status after each fit, stop after last fit of plot
2 plot status each fit, stop after each plot
3 plot status each rms loop. stop each plot
-1 same as 1 but don't wait..
fsin : int the order of cos(nx),sin(nx) to also fit.
double : if set then fit using double precision
plver[2]:float vertical range [min,max ]to use for plotting bandpasses
(before subtraction). the default is autoscaling.
RETURNS :
yfit :fltarr(n) the baselined fits or ydat-yfit if /sub set.
maskused:fltarr(n) the mask used for the fits.
coefInfo:{} structure holding the coef's for the fit of each sbc
istat : int 1 if all the fits done, 0 if one or more sbc not fit.
DESCRIPTION:
robfit_polyfft will automatically create a mask and then use it to baseline
a functionn. The input is an array of data. It will
return the masks created, the coefs of the fits, and the fits. It will
optionally remove the fits from the input data and return the difference
data-fit. The verb= option will plot the data as the fitting is done.
The routine can be used for mask creation (just ignore the fits)
or mask creation and baselining.
NOTE:
You must .compile robfit_polyfft once before using the routine since
the fitting funtion is embedded in the file.
----------------------------------------------------------------------
The parameters that control the algorithm are:
edgefract: This is the fraction of points on each edge of the data that
are not used. The default is .08 . For a sbc 1024 points this
would ignore 82 points on each edge.
masktouse: This is a mask with 0's and 1s. The 0 parts will not be
used in the fits.
deg: The maximum degree for the fit. The default is 1. Values larger
than 12 or 13 may have pivot errors in the fitting.
ndel: The number of extra channels to throw out on the left and right
of each new channel we delete. This insures that skirts of
lines don't get included in the mask.
nsig : The rms computation removes all points greater than nsig*sigma.
By default nsig is set to 3.
maxLoop: The rms removal loops until there are no points within Nsigma*sigma
or maxloop iterations have occured. By default maxLoop is 15.
---------------------------------------------------------------------------
ALOGORITHM:
0. compute the starting points to use from the edge fraction or the masktouse.
The excluded points will never be included in the mask selection or fits.
let x0 by the x points, y0 be the y points after this selection.
1. set fitOrder=0, set yfit=mean(y0)
2. set yresiduals= y0 - yfit
the rms loop starts here
3. compute the rms of yresiduals
4. remove all points in y with yresiduals gt nsig*sigma .
For each point removed, also remove ndel adjacent points (determined
by the ndel keyword).
5. if the number of points removed in 4 is gt 0 and we've looped lt maxloop
times then goto 3.
fitting .. we now have a set of points within nsig of the mean.
6. fitdeg=fitdeg+1
7. fit a polynomial of fitDeg and a harmonic function of order fsin
to the points left in y.
8. evaluate the fitted function for the entire dataset y0
yfit=poly(x0,coef).
9. If fitdeg lt keyword deg goto 2.
---------------------------------------------------------------------------
The coef structure returns the fit results:
{corcoef}:
deg : int of polynomial fit
fsin : int of harmonic fit
npnts : int number pnts in ydata
coefAr[(deg+1)+fsin*2]: float coef for each fit. the polynomial
coef are followed by the cos,sin coef.
rms: float of the fitted region within the mask
maskFract: float npointsmask/pntsData
To use the coefficients to recompute the fit, the x values should go
from -pi to +pi If Nchn is the number of channels in the spectra then
x=((findgen(Nchn)/(Nchn-1.) - .5)*2.*pi
EXAMPLES:
1. Use a polynomial fit of 3, remove the fit from the data, plot the
data and stop after each sbc done:
istat=robfit_polyfft(ydat,yfit,maskused,coefinfo,deg=3,/sub,verb=1)
; plot the datafit with the mask overplotted
plot,yfit,cmask=maskused
2. Use a polynomial fit of 2 and a harmonic fit of 2. Create a mask
to use before calling this. Plot the results without stopping for keyboard
input
istat=robfit_polyfft(ydat,yfit,maskused,coefinfo,deg=2,fsin=2,/sub,$
masktouse=masktouse,verb=-1)
On return :
coefinfo.coefar[0:2, are the c0,c1,c2 of the polynomial fit and
coefinfo.coefar[3:4] are the cos(x) sin(x) amplitudes while
coefinfo.coefar[5:6] are the cos(2x) sin(2x) amplitudes
This routine was ripped off from corblauto and made to work with a
generic float array (rather than a correlator datastructure).
SEE ALSO: corblauto
(See /pkg/rsi/local/libao/phil/gen/robfit_polyfft.pro)
NAME: rotvec - rotate a vector through an angle (in deg) SYNTAX: rvec=rotvec(vec,thetaDeg,axis=axis) ARGS: vec[m,n]: float/double data before rotation. first dim is x,y,z theta : float angle in degrees to rotate vector KEYWORDS: axis: int 1=x,2=y,3=z. default axis=3 (z axis) RETURNS: rvec[m,n]: float rotated vector DESCRIPTION Rotate a vector through the angle theta (in degrees). Return the rotated vector in rvec. The first dimension of vec, rvec can be 2 (for a 2d rotations) or 3 for 3d rotations. The axis= keyword (1=x,2=y,3=z) specifies which axis to rotate about. The default is axis=3 (z). If dim1 of vec is 2 then axis= is ignored and the rotation is about the z axis. Positive rotate is counter clockwise looking down from the positive axis of rotation. The definition rotates a vector in this direction. Rotation of the coordinate system would be in the oppposite sense.
(See /pkg/rsi/local/libao/phil/gen/rotvec.pro)
NAME:
rsspecanainp - input R&S spectrum from save file
SYNTAX: npnts=rsspecanainp(fname,frq,spc,info=info,/print)
ARGS:
lun: open to xxx.txt file.. form export trace .txt mode (tabs)
print: if set the print out the ascii info header as it is read in.
RETURNS:
npnts : 0 trouble reading file
501 points found
frq[501]: freq of trace 1 in mhz
spc[501]: spectra trace 1 in dbm
info[28] : string return the 28 ascii info lines read from the file.
DESCRIPTION:
The Rhode and Schwartz spectrum analyzer can save its
trace data in a text file. This routine will input the text file and
return the spectra and frequency.
The steps in getting at trace from the ybt250 to idl is:
. in idl
@phil
openr,lun,filename,/get_lun
npnts=rsspecanainp(lun,spc1,frq1,info=info,/print)
The
(See /pkg/rsi/local/libao/phil/gen/rsspecanainp.pro)
NAME:
ruze - evaluate the ruze formula for losses from surface errors.
SYNTAX: loss=ruze(freq,surfaceErr,errlambda=errlambda,db=db)
ARGS:
freq[] : freq in Mhz
surfaceErr : rms surface error in mm
KEYWORDS:
errlambda: if set then surfaceErr is fractions of a wavelength
db : if set then return loss as db
DESCRPIPTION:
Return loss do to surface irregularities. Use the ruze formula to
evaluate it. loss=exp((-4pi*delta/2)**2 where delta is the rms surface
error. The units are milimeters unles /errlambda is set, then the
error is fractions of a wavelength (and freq is immaterial but must be
entered).
The loss is returned as a linear value unless keyword db is set. In
that case the return value is db.
(See /pkg/rsi/local/libao/phil/gen/ruze.pro)
NAME: satinit - initialize to use the satellite prediction routines. SYNTAX: @satinit DESCRIPTION: call this routine before using any of the satellite prediction routines. It sets up the path for the directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/satinit.pro)
NAME: sbinit - initialize to use the sband tx idl routines SYNTAX: @sbinit DESCRIPTION: call this routine before using any of the sb... idl routines. It sets up the path for the idl sband directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/sbinit.pro)
NAME:
scanlist - list contents of data file
SYNTAX: scanlist,lun,recperpage,scan=scan,search=search,std=std,verb=verb
ARGS:
lun: int assigned to open file
recperpage: lines per page before wait for response. def:30
KEYWORDS:
scan: long position to scan before listing. def:rewind then list
search: if set, then search for header rec if not found at current
position
std : if set then assume std data
verb : if set then print info each record read
DESCRPIPTION:
list a summary of all of the scans in a file to stdout. This works for
correlator or ri datafiles.
The output data is:
SOURCE SCAN GRPS PROCEDURE car0 lst
Where grps is the number of recs in the file, procedure is the procedure
name that was used to take the data, car0 is the first entry in the
carr[] of the header (typically it has the step in the pattern like
on or off), and lst is the local sidereal time.
SEE ALSO:
corlist
(See /pkg/rsi/local/libao/phil/gen/scanlist.pro)
NAME:
scantype - return the type of scan
SYNTAX: type=scantype(hdr)
ARGS:
hdr: {hdr}
RETURNS:
type: 0 unknown
1 calon of calonoff pattern
2 caloff of calonoff pattern
3 position on of onoff pattern
4 position off of onoff pattern
DESCRIPTION:
Return the type of scan from the header info.
EXAMPLE:
print,corget(lun,b)
type=scantype(b.b1.h)
(See /pkg/rsi/local/libao/phil/gen/scantype.pro)
NAME:
searchhdr - position to the next available hdr in the file.
SYNTAX: istat=searchhdr(lun,maxlen=maxlen)
ARGS:
lun: int lun of file to search
KEYWORDS:
maxlen: long maximum number of bytes to read. default is
1megabyte
RETURNS:
istat: int 1 --> positioned to header, 0 no header found
DESCRIPTION
Position to the next available hdr in an AO data file. Search up to
maxlen bytes before quitting. If the header is not found, return positioned
at the input position.
(See /pkg/rsi/local/libao/phil/gen/searchhdr.pro)
NAME:
secs1970tojd - convert unixsecs(1970,mjd=mjd) to jd
SYNTAX: jd[n]=secs1970tojd(secs1970)
ARGS:
secs1970[n]:double secs1970 (utc based)
KEYWORDS:
mjd: if set then return mjd rather then
jd
RETURNS:
jd[n]: double jd for above times. if /mjd this will be mjd
DESCRIPTION:
Convert from specified secsfrom 1970 to julian date.
If /mjd is set, then return mjd rather than jd.
The routine uses the 1970secs, jddate for 2000.0
to convert.
Warning: 2014:nov. this routine was accurate to about
.06 seconds.
Note: secs=systime(/sec) returns utc secs from 1970
print,systime(0) or print,systime(0,secs) take utc 1970 secs
but adjusts to the local time zone (ast) in the output string
(See /pkg/rsi/local/libao/phil/gen/secs1970tojd.pro)
NAME:
sefdget - return telescope sefd(az,za,freq) for rcvr.
SYNTAX: stat=sefdget(az,za,freq,rcvrNum,sefdVal,date=date,zaonly=zaonly)
ARGS:
az[n]: float azimuth in degrees
za[n]: float zenith angle in degrees
freq: float freq in Mhz.
rcvrNum: int receiver number 1..16 (see helpdt feeds))
KEYWORDS:
date[2]: int [year,daynumber] to use for sefd computation. The default
is to use the current day. If they sefd curves change
with time, this allows you to access the sefd curve that
was in use when the data was taken.
zaonly: if set then only return the zenith angle dependence (averages
over azimuth)
RETURNS:
sefdval[n]: float sefd in Jy
stat: int -1 --> error, not data returned
0 --> requested outside freq range of fits. Return sefd
of the closest frequency value.
1 --> frequency interpolated sefd value returned.
DESCRIPTION:
Return the sefd (Jy) for the requested receiver at the specified
frequency and az,za. The default is to use the current sefd values. The
date keyword allows you to access a sefd curve that was valid at some
other epoch.
Fits have been done for sefd(az,za) at different frequencies for various
receivers. This routine will input the fit information and compute the
sefd for the two closest frequencies and then interpolate to the
requested frequency. The input fit data is stored in a common block so
the data does not have to be input from disc a second time unless you
pick a different receiver.
NOTE:
Some receivers have no sefd fits. They will return -1 in the status.
If a requested frequency is outside the fitted values, then the value
at the closest frequency is returned (no extrapolation is done).
If you have correlator data, you can use corhsefdget() to get the
sefd value. It will figure out the receiver number and date from the
header and then call sefdget.
Some fits only have a za dependance. In these cases just enter
arbitrary values for the azimuth (but the parameter is still needed).
For a description of the sefd calibration see:
http://www.naic.edu/~phil --> sefd curves
EXAMPLES:
lbw=5
get sefd at 1400Mhz az=120,za=10 for lbw
stat=sefdget(120.,10.,1400.,lbw,sefdval)
az=fltarr(20) ; az = 0 degrees.
za=findgen(20)+1 ; za=1..20
date=[2003,200] ; for 2003, daynumber:200
stat=sefdget(az,za,1321,lbw,sefdval,date=date)
; sefdval will be an array of 20 values for za 1 to 20 degrees and azimuth
; of 0 degrees.
to convert from daynumber to day,month,year
daynum=dmtodayno(d,mon,year)
dm =daynotodm(daynum,year)
where dm=[day,month]
SEE ALSO:sefdinpdata, gainget
(See /pkg/rsi/local/libao/phil/gen/sefdget.pro)
NAME:
sefdInpData - input sefd data for rcvr.
SYNTAX:
istat=sefdInpData(rcvNum,sefdData,fname=fname,date=date)
ARGS:
rcvNum: 1 thru 16. receiver to use (same as hdr.iflo.stat1.rfnum).
KEYWORDS:
fname: to specify an alternate data file with sefd values.
The default file is aodefdir() + 'data/sefd.datR{rcvNum}
date : [year,daynum] .. if specified the data when you want the
sefd for..default is most recent.
RETURNS:
istat: 1 ok, -1 bad file/rcvnum,data probably no fit data available.
sefdData[n]:{sefdData} return fit info. 1 structure for each frequency
DESCRIPTION:
Input the sefd fit data for all the frequncies of a particular receiver
(rcvNum). The rcvNum can be extracted from the headers with iflohrfnum().
The default datafile is aodefdir() + 'data/sefd.datR{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,
- 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
freq fitType c0 c1 c2 .... cN pol calVala calvalB
c0..cN are the fit coefficients, pol is I or A or B, CalValA,B are the
cal values used for this fit.
The structure format for {sefddata} is:
sefdData.rcvNum receiver number
sefdData.numFreq number of frequencies found
sefdData.startYr for the fit
sefdData.startDaynum for the fit
sefdData.endYr for the fit
sefdData.endDaynum for the fit
sefdData.fitI[numFreq] {azzafit} structure holding the coef and other
info for each fit
sefdData.calVal[2,numFreq] cal values used when each fit was made.
See azzafit for a description of the {azzafit} structure.
This routine is called automatically by corhsefdval and sefdget().
How the different cal routines vary:
sefdinpdata() inputs the data from disc. You must specify the rcvrnum.
It defaults to the current date. It loads a table in common
holding the fit info for all of the frequencies measured.
sefdget() Pass in the frequency and the rcvrnum. It will input the
data using sefdinpdata if necessary, do the interpolation
and return the sefd.
corhsefdval() You specify the correlator sbc header (eg b.b1.h). It will
compute the frequency and then call sefdgetl(). It returns the
sefd value.
SEE ALSO: corhsefdval, sefdget ,azzafit, azzafiteval,azzafitpr
(See /pkg/rsi/local/libao/phil/gen/sefdinpdata.pro)
NAME:
select - select elements from an array
SYNTAX: anew=select(a,startInd,stepIndex)
ARGS:
a[]: array to select elements from
startInd: long starting index in a to start selecting (0 based).
stepInd: long spacing between indices to select
KEYWORDS:
RETURNS:
aNew[]: subarray extacted from a.
DESCRIPTION:
Starting at index startInd select points spaced stepInd from the array
a. Return the subarray in anew.
EXAMPLE:
a1=select(a,0,2) .. select every other sample from a starting at first
a2=select(a,1,3) .. select every third sample from a start as the 2nd
(See /pkg/rsi/local/libao/phil/gen/select.pro)
NAME: shcolsym - show the default colors and symbols SYNTAX: shcolsym ARGS: DESCRIPTION: plot the default colors and symbols.
(See /pkg/rsi/local/libao/phil/gen/shcolsym.pro)
NAME:
shiftregcmp - compute a shift register code
SYNTAX: codelen=shiftregcmp(codeI,code,lastval=lastval)
ARGS:
codeI:{} code info structure describing the shift register code
RETURNS:
codelen: long length of code computed
code[codelen]: float : computed shift register code
lastval: ulong if keyword provided then return the last value of the
shift register (only for galois lfsr codes).
DESCRIPTION:
Compute a shift register code given the shift register description in
the codeI structure. The structure contains:
help,codeI,/st
** Structure CODEINFO, 4 tags, length=92, data length=92:
NUM_REG LONG 1 ; number of registers for this code
LEN LONG 1 ; length of code
NUM_FDBACK LONG 1 ; number of taps (feedbacsk)
FDBACK LONG Array[20] ; locations of feedbacks
The indices in this array are 1 based..
the N feedback locations are stored in
fdback[0:n-1]
galois int i ; 0==> fibonnaci, 1==> galois lfsr
startVal ulong 0 ; for galois
endVal ulong 0 ; for galois
routine that will generate the code info structure:
shiftregcmp(): generates pncodes used by the sband radar.
(See /pkg/rsi/local/libao/phil/gen/shiftregcmp.pro)
NAME: sixtyunp - unpack hhmmss.s or ddmmss.s to hh mm ss.s or dd mm ss.s SYNTAX: sixtyunp(xxyyzz,sgn,result) ARGS: xxyyzz: float,long,double either hhmmss.s or ddmmss.s RETURNS: sgn : int sign either +1 or -1 result[3]: float [hh,mm,ss.s] or [dd,mm,ss.s] DESCRIPTION: split hhmmss.s or ddmmss.s to h,m,s or d,m,s. Return the positive values in the array result. The sign of the value is returned in sgn. EXAMPLES: sixtyunp,112233.3 ,sgn,result.. sgn=1., result=[11.,22.,33.3] sixtyunp,-000033.3,sgn,result.. sgn=-1., result=[0.,0.,33.3]
(See /pkg/rsi/local/libao/phil/gen/sixtyunp.pro)
NAME:
smofrqdm_1d - freq domain smoothing (1d)
SYNTAX: yNew=smofrqdm_1d(y,fracSmo=fracSmo,ftype=ftype,$
filtertouse=filtertouse,retfilt=retfilt)
ARGS:
y[n]: float data to smooth
KEYWORDS:
fracSmo:float (0 to 1.) fraction of spatial frequency to keep.
filterType: int 1 - apply fracSmo boxcar
2 - apply fracSmo boxCar*hanning window
3 - apply fracSmo boxCar*winCos4
filterToUse[n]:float if Supplied, then ignore fracSmo,filtertype. Use this
array as the multiplier in the freq domain.
(note. Dc should be at index n/2 (count from 0)
RETURNS:
ynew[n]: float the filtered data
retfilt[n]: float the filter used in the freq domain. dc is at the center
DESCRPIPTION:
Smooth a 1 d function by multiplying the fft of the function by
a smoothing function. By default this is a boxcar window of length
fracSmo*n. You can specify a tapering of the boxcar window using
fitlertype 2 or 3. You can also specify your own smoothing function
via filterToUse.
In the frequency domain, the data is rotated so that dc is at
n/2 (count from 0) before the multiply is done.
EXAMPLES:
.. let y contain a 53 cycle and 200 cycle sine wave.
.. Then try and remove the 200 cycle cosine by using
a boxcar smoothed with a cos^4 window. The nyquist rate is
511 so 200/511. = .3914. Try a filter of .38 to see how much of the
200 cycle sin wave is left (not much).
y=mksin(1024,53) + mksin(1024,200)
fract=.38
ynew=smofrqdm_1d(y,fracSmo=fract,ftype=2)
plot,abs(fft(ynew))
WARNING: The routine computes the filter size using fix(n*fract).
The smallest measurable change is fract is 1./n
SEE ALSO:
smofrqdm_2d
(See /pkg/rsi/local/libao/phil/gen/smofrqdm_1d.pro)
NAME: spcanainit - initialize to use the spectrum analyzer routines SYNTAX: @corinit DESCRIPTION: call this routine before using any of the spectrum analyzer routines (eg routines to access spectrum analyzer files.. anritsu,...)
(See /pkg/rsi/local/libao/phil/gen/spcanainit.pro)
NAME: spwrinit - initialize to use the site power idl routines. SYNTAX: @corinit DESCRIPTION: call this routine before using any of the site power idl routines. It sets up the path for the idl spwr directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/spwrinit.pro)
NAME: srcextended - compute gain correction for extended source SYNTAX: correction=srcextended(srcWidth,fwhm) ARGS: srcWidth :float source width in same units as beamwidth fwhm[n]:float beam width (fwhm) RETURNS: correction: float When computing telescope gain use correction*SrcFlux DESCRIPTION: Compute the loss from a source being extended. The user inputs the source width and beamWidth(fwhm) in the same units. The routine assumes both of these are axially symmetric. The returned correction is a fraction of the point source value. It is linear in power.
(See /pkg/rsi/local/libao/phil/gen/srcextended.pro)
NAME:
stkscalibv - stokes calibration from complex voltages
SYNTAX: istat=stkscalibv(rcvNum,cfr,smpTimeUsec,fftlen,VcalonA,VcalOnB,VcalOffA,VcalOffB,$
calibI,lcorcal=lcorcal,date=date,degPerMhz=degPerMhz,edge=edge
ARGS:
rcvNum : int receiver number to use when querying cals: 1=327,2-430,5=lbw,$
7-sbw,8-sbh,9-cb,11-xband,13-cbw,17-alfa
cfr : float center freq in Mhz. A single cal value for the whole band
will be computed at this freq.
smpTmUsec :float sample time for data in microsecs
fftlen :long fftlen to use when computing the spectra
vcalOnA[n]:complex complex voltage data for calon polA
vcalOnB[n]:complex complex voltage data for calOn polB
vcalOffA[n]:complex complex voltage data for caloff polA
vcalOffB[n]:complex complex voltage data for calOff polB
KEYWORDS:
lcorcal: if set then use low cor cal. def: high cor cal
date[2] : lonarr [year,daynumber] date to use when searching for cal value.
def: most recent cal values
degPerMhz: float initial guess for phase change polA,B per Mhz. def:8 deg/MHz
comes from the ao fiber.
edge[2] : float fraction of bandpass edge to ignore on each side. def: .06
Calibration and fitting will exclude data on the edges.
DESCRIPTION:
Using complex sampled data, compute the intensity calibration and phase correction
for polA,polB from a calon,off measurement. The user inputs the complex voltages for
calonA,calonB, caloffA,caloffB (the routine currently assumes these are all the same length).
It will compute the complex spectra, intensity calibrate, and then fix the phase difference of polA
polB using stokesU from the cal deflection.
The program steps are:
1. look up the cal value (it uses a single value for the entire band).
2. compute the complex spectra for the 4 data sets
3. compute the average spectra for the 4 buffers
4. compute the indices i0,i1 of the spectrum to use for calibrating and fitting
(this comes from the edge keyword).
5. compute the caldeflection averaging over the band i0:i1
6. compute the calScale factor (a/d counts to Kelvins)
7. compute the bandpass correction: 1/caloff normalized to mean(caloff[i0:i1])
8. compute the scaling arrays (convert to kelvins and bandpass correct)
9. compute stokesI,Q after scaling the data from the average spectra
10 compute stokesU,V : compute stokes u,v for each fft, scale and average
11. compute the cal deflection: calon - caloff for each stokes parameter
12. fit sine waves (non linear fit) to .stokes u and V : a[0]+a[1]*sin((freqMhz*a[2] + a[3])*2*pi)
13. compute the phase difference delta for polA,B using atan(fitStksV/fitStksU)*!radeg
14. unwrap any 180 deg phase jumps in delta
15. linear fit to unwrapped delta vs freqMHz.
16. using the linear fit, compute the phase correction for u,v
phCor=exp(complex(0.,-(!dtor*poly(frqAr,calibI.phFit))))
cuv=complex(calDU,calDV)
crotated=cuv*phCor
U=float(crotated)
V=imaginary(crotated)
help,calibI,/st
CFR FLOAT 2380.00 ; freq used to compute cal value
FFTLEN LONG 1024 ; length of fft when computing spectra
BWMHZ FLOAT 0.500000 ; bandwith of data in MHz
CALK FLOAT Array[2] ; cal in degK [polA,polB]
INDEDGE LONG Array[2] ; indices into fftlen to use for calibration
UI STRUCT -> Array[1] ; sine fit to stokes u (see below)
VI STRUCT -> Array[1] ; sine fit to stokes V (see below)
DELTAF FLOAT Array[1024] ; phase difference PolA-polB vs freq (in degrees)
DELTAFU FLOAT Array[1024] ; phase difference polA,B and phase unwrapping
PHFIT FLOAT Array[2] ; linear fit to deltafU (deg/Mhz). y=phfit[0] +freqMHz*a[1]
STKCALD FLOAT Array[1024, 4];stokes i,q,u,v for caldeflection after phase correction,
and intensity calibration
uI,vI sine fit info
IDL> help,calibI.uI,/st
DATA FLOAT Array[1024] ; data used for the sine fit
YFIT FLOAT Array[1024] ; sine fit to data
COEFAR FLOAT Array[4] ; coef for fit:y=a[0]+a[1](sin((freMHz*a[2]+a[3])*2*pi)
SIGMA FLOAT Array[4] ; fit sigmas for each parameter
CHISQ FLOAT 6.04890 ; chisq from fit (but weights were all set to 1.)
Notes:
1. All of the cal deflection should come out in stokesU (since the cal is injected equally into
polA and polB). For some setups i've gotten all of the cal deflection in -U (less then zero).
In this case there is a 180 deg. phase shift in delta.
2. If you need to flip the bandpass (because the original data had an odd number of high side
lo's) then remember to change the sign of stokesV.
(See /pkg/rsi/local/libao/phil/gen/stkscalibv.pro)
NAME:
stripmask - interactively make masks for strips in a map.
SYNTAX: maskAr=stripmask(d,remavg=remavg)
ARGS:
d[m,n]: float map of m samples per strips and n strips
KEYWORDS:
remavg: if set then remove the median from each strip before
displaying a strip. In this case you may not have to fiddle
with the vertical scale on each strip.
RETURNS:
maskAr[m,n]: int holds n masks of 0,1's
DESCRIPTION:
Let the user interactively define masks for a number of strips of a
map. For each of n strips, call bluser() and allow the user to define
a mask array using the cursor. When all n strips have been done,
return the maskAr.
If the user does not define a mask for a particular strip, then the
mask from the previous strip will be used.
This routine calls bluser() but it is the users responsibility to
enter the keys:
m .. then define the mask with the cursor
q .. to exit from bluser for each strip
EXAMPLE:
Suppose you've call cormapinp and you want to create a mask that
does not include continuum sources in the map. Use the total power in
polA+polB to find the continuum. Suppose there are 120 samples per
strip and 36 strips. The following code will call bluser 36 times.
istat=cormapinp(lun,scan,brdA,brdB,m,cals); input the map
polAvg=total(m.p)/2. ; average pola,polB
ver, -.001,.015 ; vertical scale for plot
maskArr=stripmasks(polAvg,/remavg)
.. maskArr will now be dimensioned maskArr[120,36]
SEE ALSO:
bluser()
(See /pkg/rsi/local/libao/phil/gen/stripmask.pro)
NAME:
strips - plot strips with offset and increment versus sample.
SYNTAX: strips,y,offset,step,smo=n,title=title,xtitle=xtitle,
ytitle=ytitle
ARGS:
y[m,n]: 2d data to plot.
offset: float offset to add to the first line plotted.
step : float value to add to plot the next line.
KEYWORDS:
smo : int smooth each line by n before plottting.
title : string title of plot
xtitle : string label for x axis
ytitle : string label for y axis
DESCRPIPTION:
Plot the 2d array y line by line. Offset the first line by off and
then space each line by step.
EXAMPLE:
dat[100,20] is the data
strips,dat,0,.02
You should setup the vertical scale with ver first.
SEE ALSO:
stripsxy,ver,hor
(See /pkg/rsi/local/libao/phil/gen/strips.pro)
NAME:
stripsxy - plot strips with offset and increment verus x.
SYNTAX: stripsxy,x,y,offset,step,over=over,smo=tosmo,dec=dec,title=title,$
xtitle=xtitle,ytitle=ytitle,stepcol=stepcol,colar=colar,_extra=e
ARGS:
x[n] : float data for x axis.
y[n,m]: float 2d data to plot. m strips of n points.
offset: float offset first strip by this amount.
step : float separate each strip by this amount.
KEYWORDS:
over : if set then continue plotting from a previous call.
smo : int number of points to smooth each line.
dec : int number of points to decimate each each line.
title : string title of plot
xtitle: string xlabel
ytitle: string ylabel
colAr[]: int lut indices to use for color
stepcol: if set then alternate colors between each line. Use color
indices 1-10. use ldcolph to load the color table.
_extra: pass this to the plot,oplot routine.
(See /pkg/rsi/local/libao/phil/gen/stripsxy.pro)
NAME:
strtovarnam - modify a string to be a valid variable name
SYNTAX: newstr=strstovarnam,str,noleading=noleading
ARGS:
str[]: str hold the string or strings to check
KEYWORDS:
noleading: if set then don't bother checking if the leading
character is a letter.
RETURNS:
newstr[]: str modified strings that are valid variable names.
DESCRIPTION:
This routine will modify the strings passed in to be valid variable
names. A variable name must start with a letter and not contain the
characters: , . ! ; : + - / =
This routine will make the following substitutions:
, -> _
. -> _
+ -> p lus
- -> m inus
/ -> _
= -> _
! -> _
) -> _
( -> _
If the string starts with a non letter then prepend the string with
the letter v.
This routine is normally used when creating a variable name
that is taken from the source name in the header.
EXAMPLE:
a='B1934+123'
aMod=strtovarnam(a)
print,aMod
.. 'B1934p123'
a='1934-123:17'
aMod=strtovarnam(a)
print,aMod
.. 'v1934m123_17'
(See /pkg/rsi/local/libao/phil/gen/strtovarnam.pro)
NAME: suninit - initialize to use idl sun mapping routines SYNTAX: @suninit DESCRIPTION: call this routine before using any of the sunxxx idl routines. It sets up the path for the idl sun directory, initilizes mas routines and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/suninit.pro)
NAME:
SVDFITpp- Perform a general least squares fit. patched version
PURPOSE:
Perform a general least squares fit with optional error estimates.
This version uses the Numerical Recipies (2nd Edition) function
SVDFIT. A user-supplied function or a built-in polynomial or
legendre polynomial is fit to the data.
CATEGORY:
Curve fitting.
CALLING SEQUENCE:
Result = SVDFIT(X, Y, [M])
INPUTS:
X: A vector representing the independent variable.
Y: Dependent variable vector. This vector must be same length
as X.
OPTIONAL INPUTS:
M: The number of coefficients in the fitting function. For
polynomials, M is equal to the degree of the polynomial + 1.
If not specified and the keyword A is set, THEN
M = N_ELEMENTS(A).
INPUT KEYWORDS:
A: The inital estimates of the desired coefficients. If M
is specified, THEN A must be a vector of M elements.
If A is specified, THEN the input M can be omitted and
M=N_ELEMENTS(A). If not specified, the initial value
of each coefficient is taken to be 1.0. If both M and A are
specified, them must agree as to the number of paramaters.
DOUBLE: Set this keyword to force double precision computations. This
is helpful in reducing roundoff errors and improves the chances
of function convergence.
WEIGHTS: A vector of weights for Y[i]. This vector must be the same
length as X and Y. If this parameter is ommitted, 1's
(No weighting) are assumed. The error for each term is
weighted by Weights[i] when computing the fit. Gaussian or
instrumental uncertianties should be weighted as
Weight = 1/Sigma where Sigma is the measurement
error or standard deviations of Y. For Poisson or statistical
weighting use Weight=1/Y, since Sigma=sqrt(Y).
FUNCTION_NAME:
A string that contains the name of an optional user-supplied
basis function with M coefficients. If omitted, polynomials
are used.
The function is called: R=SVDFUNCT(X,M)
where X and M are scalar values, and the function value is an
M element vector evaluated at X with the M basis functions.
M is the degree of the polynomial +1 if the basis functions are
polynomials. For example, see the function SVDFUNCT or SVDLEG,
in the IDL User Library:
For more examples, see Numerical Recipes in C, second Edition,
page 676-681.
The basis function for polynomials, is R[j] = x)^j.
The function must be able to return R as a FLOAT vector or
a DOUBLE vector depending on the input type of X.
LEGENDRE: Set this keyword to use the IDL function SVDLEG in the lib
directory to fit the data to an M element legendre polynomial.
This keyword overrides the FUNCTION_NAME keyword.
OUTPUTS:
SVDFIT returns a vector of the M coefficients fitted to the
supplied function.
OPTIONAL OUTPUT PARAMETERS:
CHISQ: Sum of squared errors multiplied by weights if weights
are specified.
COVAR: Covariance matrix of the coefficients.
VARIANCE: Sigma squared in estimate of each coeff(M).
That is sqrt(VARIANCE) equals the 1 sigma deviations
of the returned coefficients.
SIGMA: The 1-sigma error estimates of the returned parameters,
SIGMA=SQRT(VARIANCE).
SINGULAR: The number of singular values returned. This value should
be 0. If not, the basis functions do not accurately
characterize the data.
YFIT: Vector of calculated Y's.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
MODIFICATION HISTORY:
Adapted from SVDFIT, from the book Numerical Recipes, Press,
et. al., Page 518.
minor error corrected April, 1992 (J.Murthy)
Completely rewritten to use the actual Numerical Recipes routines
of the 2nd Edition (V.2.06). Added the DOUBLE, SIGMA, A, and
LEGENDRE keywords. Also changed Weight to Weights to match the
other fitting routines.
(See /pkg/rsi/local/libao/phil/gen/svdfitpp.pro)
NAME:
symcircle - create a circle symbol for plots
SYNTAX: symcircle,npnts=npnts,color=color,thick=thick,fill=fill
KEYWORDS:
npnts: int number of points to use in drawing the circle. The default
is 9 points. The circle will be drawn as an n-1 sided
polyhedron.
color: int number 0 to 10. Color to use when plotting symbol. The
default is the same color as the lines beging drawn. See
shcolsym for a mapping of numbers to colors.
thick:float the line thicknes to use when drawing the lines. 1. is the
default.
fill : If set then fill the circles
DESCRIPTION:
Create a cirlce symbol to use for plotting. After calling this routine,
plot,....psym=-8 will draw a circle about every point.
EXAMPLE:
symcircle,npnts=17 .. use 17 points to define circle
plot,findgen(30),psym=-8
(See /pkg/rsi/local/libao/phil/gen/symcircle.pro)
NAME:
tempplot - plot temperature in turret room for a range of days.
SYNTAX:
tempplot,yymmdd1,yymmdd2 ,title=title,cs=cs,xp=xp,off=off,lun=lun,sep=sep,$
clip=clip
ARGS:
yymmdd1 - long first day to plot
yymmdd2 - long last day to plot
KEYWORDS:
title - string title for top of each plot
cs - float scale factor for labels. For multiple windows you may
want to increase this to 1.5 or 1.8. default is 1.
xp - float xposition [0.,1.] where dates are printed. default:.02.
For hardcopy setting this to 1. puts the dates on the
right column.
off - float number of degrees to offset each day within it's
window so the plots don't lay on top of each other.
defaut:0.
lun - int if supplied, then you've already opened the input file.
Use this to access previous years:
eg: /share/obs4/temp/Y01/temp.dat
sep - if set then make 1 plot per day user should set
!p.multi before calling the routine.
clip[] -float clip the data to [mintemp,maxtemp]. allows auto scaling
with bad data points..
DESCRIPTION
plot the temperature in the turret room by hour of day for the
range of dates specified (the routine limits it to a maximum of 31 days).
The plot will have N subwindows with up to 7 consecutive days per window.
Each day will be color coded witin a window. Any extra days (>28 will appear
in the last window). The routine internally uses daynumber of year for
computations so it will not cross year boundaries gracefully. It also reads
from the temperature file that holds info by year (so the current year
is the only one that will work.
You should call ldcolph to setup the colortables for indices 0-10. You
should also set the vertical scale to temperature range you want via
ver,vmin,vmax. The horizontal scale should be set to hor,0,24.
The extra keywords that may be used are:
charsize=cs. If the multiple windows cause the letters to come out
too small.
EXAMPLES:
hor,0,24
ver,70,90
tempplot,010701,010731,title='turret room temps for jul01',$
cs=1.6
NOTES:
This routine reads from the file /share/obs4/temp/temp.dat. The routine
will only work when run at AO (since the file is not accessible at
remore sites).
The routine will not cross year boundaries very gracefully. The
data files usually contain only one years worth of data. The previous
years data is in (/share/obs4/temp/yyyy/temp.dat where yyyy is the year.
Tempread has not yet been updated to work on a pc (you need to run it
on a sun (big endian machine).
SEE ALSO:
tempread.
(See /pkg/rsi/local/libao/phil/gen/tempplot.pro)
NAME:
tempread - read receiver room temperature data
SYNTAX: dat=tempread(yymmdd,lun=lun,nrec=nrec)
ARGS:
yymmdd: long day of year to input.
KEYWORDS:
lun : int if provided then this is open to the file to read from.
The default is to read from the temperature data file of
the current year.
nrec : long . if present, read this many records starting at yymmdd.
the default is to read 1 days worth of data.
RETURNS:
dat[n]: {tempdat} temp data structure returned.
dat.time : float daynumber of sample
dat.temp : float temp in degree F.
(See /pkg/rsi/local/libao/phil/gen/tempread.pro)
NAME:
tosecs1970 - convert to unixsecs(from 1970)
SYNTAX: secs1970=tosecs1970(yymmdd,hhmmss,mjd=mjd,jd=jd,daynoyr=daynoyr)
ARGS:
yymmdd[n]: long year,mon,day for time (utc based)
hhmmss[n]: double hour,min,sec for time (utc based)
KEYWORDS:
jd : double if provided then ignore yymmdd,hhmmss and use jd
as the date to convert to secs since 1970
mjd: double if provided then ignore yymmdd,hhmmss and use mjd
as the date to convert to secs since 1970
daynoYr[2]:double [dayno,year] if provided then ignore yymmdd,hhmmss
and use dayno and year (utc based) to convert to secs1970.
DESCRIPTION:
Convert from specified date, time to secs from 1970 (unix time).
If jd,mjd, or daynoyr are provided then used them instead of yymmdd,hhmmss.
All times/dates provided are assumed to be utc based. If you want to enter an
ast time then :
secs1970=tosecs1970(080916,142022) + 4D * 3600D
Note: secs=systime(/sec) returns utc secs from 1970
print,systime(0) or print,systime(0,secs) take utc 1970 secs
but adjusts to the local time zone (ast) in the output string
(See /pkg/rsi/local/libao/phil/gen/tosecs1970.pro)
NAME: tsysinit - initialize idl to process system temperature monitoring data. SYNTAX: @tsysinit ARGS: none DESCRIPION: Initialize to process tsys data taken daily with the tsysall program online.
(See /pkg/rsi/local/libao/phil/gen/tsysinit.pro)
NAME:
tvfreq - return tv channels and frequencies
SYNTAX: tvfreq,chan1,chan2,chan,freq,sound=sound,chroma=chroma,cenfrq=cenfrq
ARGS:
chan1 - int first channel number (2-69)
chan2 - int last channel number (2-69)
RETURNS:
chan[] - long channel numbers
freq[] - float frequencies
KEYWORDS:
sound : int if set then return the sound carrier frequency
chroma: : int if set then return the color carrier
cenfrq: : int if set then return the center of the band. This should
be used for digital channels.
DESCRIPTION:
return the frequencies and channel numbers for the tv channels between
chan1 and chan2 (inclusive). Return the channel numbers in chan and the
frequencies in freq. By default the picture carrier frequency is returned
(1.25 Mhz above the lower edge). If the sound keyword is set then return
the frequency of the sound carrier (5.75 mhz above the lower edge). If the
chroma keyword is set then return the color carrier (3.579656 mhz above
the picture carrier)..
EXAMPLE:
tvfreq,2,20,chan,freq
chan[] will contain the numbers 2-20
freq[] will contain the picture carriers for each channel.
(See /pkg/rsi/local/libao/phil/gen/tvfreq.pro)
NAME: usrpinit - initialize to access usrp idlroutines. SYNTAX: @corinit DESCRIPTION: call this routine before using any of the usrp idl routines. It sets up the path for the idl correlator directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/usrpinit.pro)
NAME: usrprojinit - addpath holding routines to user projects SYNTAX: @userprojinit DESCRIPTION: This adds the path for the routines in ./usrproj . These routines are written for particular projects.
(See /pkg/rsi/local/libao/phil/gen/usrprojinit.pro)
NAME: ver - set vertical scale for plotting. SYNTAX: ver,ver1,ver2 ARGS: ver1: float min vertical value ver2: float max vertical value. DESCRIPTION: Load the !y.range system value with the min,max yrange to plot. To reset to auto scaling call ver with no args. SEE ALSO: hor
(See /pkg/rsi/local/libao/phil/gen/ver.pro)
NAME:
waitnxtgrp - wait for next group from the file to become available
SYNTAX:
istat=waitnxtgrp(lun,[maxwaitSecs],bytesingrp=bytesingrp)
ARGS:
lun: assigned to file
OPTIONAL ARGS:
maxwait: maximum number of seconds to wait before returning. Default
is 99999
RETURNS:
istat: return status. 0 ok, -1 timedout, -2 not lined up with a header.
bytesingrp: long bytes in the next group. You can use this to position
to the end of the new group
DESCRIPTION:
Wait for the next group from the file to be available. Return 0
if ok, -1 if timeout. You can then read the data in with corget. On
return, the lun is left positioned at the start of the group to read.
EXAMPLE:
Assume you are monitoring the online datafile.
if waitnxtgrp(lun) ne 0 then .. error message
istat=corget(lun,b) ;input the data
SEE ALSO:
corget
(See /pkg/rsi/local/libao/phil/gen/waitnxtgrp.pro)
NAME: wappinit - initialize to use the idl wapp pulsar routines. SYNTAX: @wappinit DESCRIPTION: call this routine before using any of the wapp idl routines. It sets up the path for the idl wapp directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/wappinit.pro)
NAME: wbffainit - initialize to use the wideband rcv and field fox SYNTAX: @wbffinit DESCRIPTION: call this routine before using any of the wbff 12m wbrcv with field fox routines (eg routines to access spectrum analyzer files.. anritsu,...)
(See /pkg/rsi/local/libao/phil/gen/wbffinit.pro)
NAME: wincos4 - generate a cos^4 window SYNTAX: win=wincos4(length) ARGS: length: long number of point in the window RETURNS: win[length]: window function. DESCRIPTION: Create a cos^4 window function of length points. The maximum value is normalized to 1. SEE ALSO: windowfunc()
(See /pkg/rsi/local/libao/phil/gen/wincos4.pro)
NAME: windinit - initialize idl to process wind monitoring data. SYNTAX: @windinit ARGS: none DESCRIPION: Initialize to process wind data.
(See /pkg/rsi/local/libao/phil/gen/windinit.pro)
NAME:
windowfunc - make a window function
SYNTAX: val=windowfunc(len,double=double,type=type)
ARGS:
len: long length of window function
KEYWORDS:
double: if set the return a double array. the default is float
type: char type of window:
'cos4' cosine^4 window
'ecb' extended cosine bell
'hsin' half sine
'han' hanning window
'ham' hamming window
The default is a hanning window
DESCRIPTION:
Create a window function that can be used in fourier transform
processing
(See /pkg/rsi/local/libao/phil/gen/windowfunc.pro)
NAME: wst - initialize to use idl weather station routines SYNTAX: @wst DESCRIPTION: call this routine before using any of the wst... orion weatherstation idl routines. This routine sets up the path for the idl wst directory and defines the necessary structures.
(See /pkg/rsi/local/libao/phil/gen/wstinit.pro)
NAME: wuse - set window for plotting SYNTAX: winuse,winnum ARGS: winnum: int window number to use KEYWORDS: Any extra keywords are passed to window,winnun DESCRIPTION: Set plotting window to winnum. If the window is not currently available then call window,winnum,_extra=e If the window is currently available then call wset,winnum
(See /pkg/rsi/local/libao/phil/gen/wuse.pro)
NAME: x - set output to xwindows device SYNTAX: x ARGS: none DESCRIPTION: Set plot device to the xwindows terminal. SEE ALSO: ps,pscol,hardcopy.
(See /pkg/rsi/local/libao/phil/gen/x.pro)
NAME:
x102combineps - combine x102 ps files into 1 file.
SYNTAX: x102combineps(m,fdscript,psdir,yymmdd1=yymmdd1,yymmdd2=yymmdd2)
ARGS: m[] : string array of filenames of a particular plot type to
combine.
fd : int .. script file to write to
psdir : string.. directory to hold the combined postscript files
KEYWORDS:
yymmdd1: long if present then first data to process
yymmdd2: long if present then last date to process
DESCRIPTION:
The mueller2.idl routine for x102 calibration creates a large number of
postscript files. For receivers with multiple frequencies it will create
1 separate file for each board or frequency. This routine will create a
csh file that will combine the N frequencies from a particular receiver into
a ps file with N pages (1 for each frequency).
To use:
1. go to the directory that has the postscript files.
2. in idl create a string array that has the plot type of the plots
to combine:
eg:
m4=spawn,'ls lbn*m4*'
2. open the script file that will hold the commands to combine the files
openw,lun,'combineps.sc',/get_lun
When done you will execute this file in the shell to do the combination.
3. decide on an output directory where you want the combined ps files
written.
x102combineps,lun,'dirpath')
free_lun,lun
4. get out of idl, look at the file to make sure it is ok.
change the mode to executable, chmod +x combineps.sc
and then execute it.. combineps.sc
The filenames look like:
lbn_1300_bd0_B0035+130_m4_22-APR-2001.ps
The routine strips off the first 3 sections to have:
B0035+130_m4_22-APR-2001.ps
It gathers together all of the files that have this base filename. For each of
these sets it sorts the full name by the second section _1300_
(the frequency).
The outputfile it chooses is:
lbn_all_B0035+130_m4_22-APR-2001.ps where B0035+130_m4_22-APR-2001.ps is the
base name. The script file will then look like:
psidlmerge -o psdir/lbn_all_B0035+130_m4_22-APR-2001.ps \
then a list of all frequencies for this source..
The script psidlmerge is in ~phil/Solaris/bin..
(See /pkg/rsi/local/libao/phil/gen/x102combineps.pro)
NAME: x111init - initialize idl to process x111 data. SYNTAX: @x111init ARGS: none DESCRIPION: Initialize to process x111 data. This is interference monitoring data taken on the telescope. The routine calls @corinit, adds the x111 path SEE ALSO: x111doc
(See /pkg/rsi/local/libao/phil/gen/x111init.pro)
NAME:
xvisualquery - query x visual information
SYNTAX: vis=xvisualquery()
ARGS:
none
RETURNS:
vis: {} structure containing x visual info found
vis.retain 0/1 ; 1--> server provides backing store
vis.pseudoCol.numvis ; number of pseudo color visuals found
vis.pseudoCol.numplanes[f]; number of planes in each of these
vis.directCol.numvis ; number of direct color visuals found
vis.directCol.numplanes[f]; number of planes in each of these
vis.trueCol.numvis ; number of true color visuals found
vis.trueCol.numplanes[f] ; number of planes in each of these
DESCRIPTION:
In idl you would like to know what kind of x visuals are available on
the x display before requesting one. The routines to query the visuals
have the unfortunate drawback that if no visual is selected, they will
select one.
This routine executes xdpyinfo and then parses the output. It keeps
track of how many pseudocolor, directcolor, and truecolor visuals were
found. For this routine, the visuals are differentiated only by the
class (pseudo, true, direct) and the number of planes (it does not
differentiate between visuals with the same class and number of planes.
The routine also returns whether the xserver can supply backing
store (vis.retain eq 1) or not (vis.retain eq 0). You could then
use this value to determine if idl should provide backing store or not.
(See /pkg/rsi/local/libao/phil/gen/xvisualquery.pro)
NAME:
ybt250inp - input ybt250 spectrum from save file
SYNTAX: ntrace=ybt250inp(lun,spc1,frq1,spc2,frq2,info=info,/print)
ARGS:
lun: open to xxx.txt file.. form export trace .txt mode (tabs)
print: if set the print out the ascii info header as it is read in.
RETURNS:
ntrace: 0 trouble reading file
1 1 trace found ..returned in spc1,frq1
2 2 trace found ..returned in spc1,frq1 and spc2,frq2
frq1[501]: freq of trace 1 in mhz
spc1[501]: spectra trace 1 in dbm
frq2[501]: freq of trace 2 in mhz
spc2[501]: spectra trace 2 in dbm
info[29] : string return the 29 ascii info lines read from the file.
DESCRIPTION:
The textronix ybt250 portable spectrum analyzer can save its
trace data in a text file. This routine will input the text file and
return the spectra and frequency.
The steps in getting at trace from the ybt250 to idl is:
1. save the trace:
On the ybt250 use:
-->file
--> save trace as
select tab separated .txt as the save option
You can change the name using the keyboard button.
2. Copy the saved file to the floppy on the ybt250:
From outside of ybt250 program (in windows)
- insert floppy
--> start programs button
-- start floppy
--> netTek icon
--> builtindisc
--> ybt250
?? may have left out a directory here..
-->appresults
copy the file you want (edit, copy)
--> back to top level of netTek
--> click on floppdisc
--> edit , copy
3. copy the file to unix/linux.. find a floppy that works ..
on solaris try volcheck (does not work on lots of machines).
4. in idl
@phil
openr,lun,filename,/get_lun
ntrace=ybt250inp(lun,spc1,frq1,spc2,frq2,info=info,/print)
(See /pkg/rsi/local/libao/phil/gen/ybt250inp.pro)
NAME: yfactor - compute Tdev for hot,cold load measurements SYNTAX: Tdev=yfactor(yhot,ycold,thot,tcold) ARGS: yhot[n]: float spectral density on hot load ycold[n]: float spectral density on cold load Thot : float temperature host load Tcold : float temperature cold load RETURNS: Tdev[n]: float device temp vs freq DESCRIPTION: Use the y factor measurement to compute the device under test temperature. The inputs are the spectral density measured on hot and cold loads. You must specify the temperature of the hot and cold load. The equations are: let c,h be hot,cold, let yc,yh be the meausured spectra, and g is the gain between loads and spectroeter/spectrum analyzer (Tc + Td)*g =Yc (Th + Td)*g =Yh yh/yc= (Th + Td)/(Tc +Td) Yf=yH/Yc Yf*(Tc + Td)= (Th + Td) Td*(Yf-1) + Yf*Tc=Th Td=(Th - Yf*Tc)/(Yf -1) or Td=(Th - Yh*Tc/Yc)/(Yh/Yc-1) Td=[(Th*Yc - Yh*Tc)/Yc] / [ (yh-yc)/yc] Td=(Th*Yc - Yh*Tc)/(yh - yc)
(See /pkg/rsi/local/libao/phil/gen/yfactor.pro)
NAME:
yymmddtodmy - convert yymmdd to ddMonyy
SYNTAX: dddMonyy=yymmddtodmy(yymmdd)
ARGS:
yymmdd: long convert yymmdd to ddMONyy (eg. 020210 to 10mar02)
RETURNS:
ddMONyy:string return '' if bad format..
DESCRIPTION
Convert from yymmdd long to ddMONyy string
EXAMPLE:
yymmdd=060428
ddMONyy=yymmddtodmy(yymmdd)
print,ddMONyy
28apr06
(See /pkg/rsi/local/libao/phil/gen/yymmddtodmy.pro)
NAME:
yymmddtojulday - convert yymmdd to julian day
SYNTAX: julday=yymmddtojulday(yymmdd)
ARGS:
yymmdd[]: long to convert
RETURNS:
julday[]: double julian day
DESCRIPTION:
Convert from yymmdd to julian day.
The input can be a scalar or an array.
(See /pkg/rsi/local/libao/phil/gen/yymmddtojulday.pro)