GBTIDL toolbox commands

GBTIDL toolbox commands#

`accumave`	Get the average from an ongoing accumulation.
`accumclear`	Clear the given accum_struct, free's the pointer, zeros' values.
`awv`	Function to find the area, width, and velocity of a galaxy profile.
`chantofreq`	Convert channel number to frequency (Hz) using
`chantovel`	Convert channel number to velcity (m/s) using
`chdoppler`	Computes the projected velocity of the telescope wrt
`compress_ints`	Take an integer array, sort it, and convert it to a comma separated
`coord_interp`	Find the value of the coordinate at the given location, watching for
`count_blanks`	Count the number of blanked values in the data array of the given
`DATA_COPY`	Copy the data container from in to out.
`DATA_FREE`	Free the data pointer in a data structure. This should only be
`DATA_NEW`	Create a new data_struct of the requested data type (spectrum or
`DATA_VALID`	Check to see if the argument is a recognized data structure.
`daynotodm`	Convert daynumber and year to day of month (1..31) and month
`dcaccum`	Add a data container into an ongoing accumulation in a given
`dcadd`	This procedure adds two data container's data arrays, returning the sum
`dcascii`	Output a data container as ascii showing the x-axis value and the
`dcbaseline`	Fit a baseline to the contents of the dc argument with various options.
`dcbias`	This procedure adds a bias factor to a data container's data
`dcboxcar`	Smooth a data container with a boxcar smoothing of a certain
`dcclip`	This procedure truncates (clips) data above and below given
`dcconvol`	A front-end to the standard IDL CONVOL function.
`dcdecimate`	This procedure decimates a spectrum (i.e. thins the spectrum by
`dcdivide`	This procedure divides the first data container's data array by the second, returning the result.
`dcdofeedoffsets`	Apply or remove the feed offsets in the supplied data container.
`dcextract`	Extract a region from a data container, producing a new data
`dcfft`	Do an FFT (forward or inverse) of a data container. For a
`dcfold`	Average two parts of in-band frequency switching (signal and
`dcfshift`	Function to calculate the shift, in channels, necessary to align in
`dchanning`	This procedure smooths a spectrum with a hanning filter.
`dcinterp`	Interpolate across blanked channels in a data container.
`dcinvert`	Flip the data end-to-end in the supplied data container.
`dcmeantsys`	Calculate the mean Tsys using the data from two spectral line data
`dcmediansub`	Subtract the median filtered values of the given width, in channels,
`dcmultiply`	This procedure multiplies two data container's data arrays, returning the product
`dcpaircheck`	This procedure checks both data containers to see if their data arrays can be used
`dcpowspec`	Compute the power spectrum for the supplied data container, which
`dcresample`	Resample a spectrum using one of 4 possible interpolation methods
`dcresize`	Change the size of an existing data container.
`dcscale`	This procedure scales a data container by a scalar value
`dcsetfframe`	Procedure to reset the frequency reference frame of the supplied data
`dcsetunits`	Set the units and convert the data to those units using the model
`dcshift`	Procedure to shift data in a data container by a given number of
`dcsmooth`	Smooth the data with a GAUSSIAN such that the spectrum in the given
`dcspurchans`	Return a vector of VEGAS ADC spur channels given the values of VSPDELT,
`dcspurinterp`	Replace data values in the given data containers at the known
`dcsubtract`	This procedure subtracts the second data container's data array from the first, returning the results.
`dcvshift`	Function to calculate the shift, in channels, necessary to align in
`dcxshift`	This function returns a value that, when given as the argument to
`DECODE_VELDEF`	Parse the SDFITS VELDEF value into its two components, the velocity
`decompress_ints`	Converts a string containing an integer list in a special syntax into
`dmtodayno`	Convert from day of month, month, and year to daynumber of year.
`doboxcar1d`	Do a boxcar smooth on a 1D array.
`docal`	This procedure gets the cal signal from a single cal-switched
`dofreqswitch`	This procedure calibrates a single integration from a frequency
`dofullsigref`	This procedure calibrates a single integration using four different
`donod`	This procedure calibrates a single integration from a Nod scan pair.
`dosigref`	Combine signal and reference data containers.
`dototalpower`	This procedure calibrates a single integration from a total power
`eqtogal`	Translate equatorial coordiantes to galactic coordinates.
`estboxres`	Estimate the new frequency_resolution after boxcar smoothing
`esthanres`	Estimate the new frequency_resolution after hanning smoothing
`extract_edges`	Given a string containing a comma separated list with optional range
`file_exists`	Check if a file name exists. Return 1 if it does, 0 if it doesn't.
`fitsdateparse`	Convert a FITS date string into its component parts.
`frame_velocity`	Toolbox front-end to chdoppler. Get the velocity of a given frame
`freqtochan`	Convert frequency (Hz) to channel number using
`freqtofreq`	Convert a frequency in one rest frame to an equivalent frequency in
`freqtovel`	Convert frequency to velocity (m/s) using the given rest
`FREQUENCY_AXIS`	Construct a frequency axis from a data structure.
`galtoeq`	Translate galactic coordinates to equatorial coordiantes.
`gauss_fits`	This code came mostly from Tom Bania's GBT_IDL work. Local
`gauss_fx`	The Gaussian function required by curvefit, and used by gauss_fits.
`gauss_plot_fn`	Generate x, y pairs to be plotted (x in channels) given a
`get_ap_eff`	Function to return a default aperture efficiency, given an
`get_chans`	Return a vector of channel numbers appropriate for the given data container.
`get_tau`	Function to return a default zenith opacity, given an
`get_usage`	Parses comment lines like this one to create a sensible usage statement.
`getazel`	Get the Azimuth and Elevation from a data container's longitude_axis
`GETDCDATA`	Convenience function for retrieving data from a data container.
`getgal`	Get galactic coordinates from a data container.
`getha`	Return the hour angle (HA) in decimal hours using the given
`gethadec`	Get the HA and DEC at the requested equinox from a data container.
`getradec`	Get the RA and DEC at the requested equinox from a data container.
`ha2ra`	Convert an hour angle (HA) in decimal hours to a Right Ascension
`isleapyear`	Determine whether a year is a leap year in the gregorian calendar.
`juldaytolmst`	Convert from Julian day to local mean sidereal time. By default the
`leftjustify`	Function to return the input string, left justified at the desired
`make_gauss_data`	Returns a gaussian, given coefficients and an x axis. Noise and offset can be added.
`makefitsdate`	Make a FITS date string from a modified Julian Date.
`ortho_fit`	Function uses general orthogonal polynomial to do least squares
`ortho_poly`	function returns the sum of the set of orthogonal polynomials
`paddedstring`	Format a number as a string, ensuring that it is padded with a
`ra2ha`	Convert a Right Ascension (RA) in degrees to an hour angle (HA) in
`scalevals`	Scale a value and determine the appropriate prefix for the scaled
`select_data`	Select data from the given io object and return the array of
`seq`	Generates a sequence of numbers in an array
`setdcdata`	Convenience function for setting data array of a data container.
`shiftfreq`	Shift a frequency using a velocity offset. Note that because of the
`shiftvel`	Shift a velocity using a velocity offset. Note that because of the
`unshiftfreq`	Undo the effects of shiftfreq using the same voffset (do not use its
`unshiftvel`	Undo the effects of shiftvel using the same voffset (do not use its
`usage`	Print out the usage for a given procedure or function.
`utcinfoinp`	Input parameters to convert utc to ut1
`utcToUt1`	utctout1 - convert utc to ut1
`VELOCITY_AXIS`	Construct a velocity axis from a data structure.
`veltochan`	Convert velocity (m/s) to channel number using
`veltofreq`	Convert velocity(m/s) to frequency using the given rest
`veltovel`	Convert a velocity in one velocity definition to an equivalent
`which`	Search for any file in the IDL !path that contains the
`which_routine`	Search for any file in the IDL !path that contains the

pro accumave(accumbuf, dc, noclear=noclear, quiet=quiet, wtarray=wtarray, count=count)#

Get the average from an ongoing accumulation.

The dc argument will contain the result. The wtarray keyword will contain the value of the weight array used at the time of the average. This can be used in a subsequent dcaccum or accum call.

The contents of the accum buffer are cleared as a consequence of calling accumave, unless the noclear keyword is set.

Note

It is a good idea to use accumclear to clear the accum buffer before using it the first time so that you can be certain it is starting from a cleared (empty) state.

Note

The frequency_resolution in dc will be the maximum resolution in all of the data containers used in the accumulation.

Params:

accumbufin, out, required, type=accum_struct: The accumulation buffer to use.
dcout, required, type=spectrum: The resulting average.

Keywords:

noclearin, optional, type=boolean, default=F: When set, the contents of the global accum buffer are not cleared. This is useful when you want to see what the current average is but also plan on continuing to add data to that average. If this is not set, you would need to restart the accumulation to average more data.
quietin, optional, type=boolean, default=F: Normally, accumave announces how many spectra were averaged. Setting this turns that announcement off. This is especially useful when multiple accum buffers are used within a procedure.
wtarrayout, optional, type=float: Optionally return the weight array (same length as the data in the dc parameter) used in constucting the average.
countout, optional, type=integer: The number of spectra that have been averaged. Returns -1 on an error. If count=0 then the other output arguments (dc and wtarray) are unchanged by this procedure.

Examples:

a = {accum_struct}
accumclear, a  ; not necessary here, but good practice
get,index=1
dcaccum, a, !g.s[0]
get,index=2
dcaccum, a, !g.s[0]
accumave, a, myavg
show, myavg
data_free, myavg ; be sure and clean up when done
; or the same step at the end, but get the average
; note that the previous dcaccum etc. must be repeated here
; since accumave cleared the average
accumave, a, myavg, wtarray=myweight
show, myavg
data_copy, myavg, mywt
setdcdata, mywt, wtarray  ; put the weight array in a copy of the data
show, mywt ; show the weight array
; clean up when done
data_free, myavg
data_free, mywt

Note

the accum_struct structure used here has internal pointers. Use accumclear to clear them, either implicitly (by use of accumave) or explicitly.

Uses:: accumclear DATA_COPY

pro accumclear(accumbuf)#

Clear the given accum_struct, free’s the pointer, zeros’ values.

Params:

accumbufin, out, required, type=accum_struct structure: the structure to clear.

function awv(data, vel, brange, erange, mode, fract, lefthorn=lefthorn, righthorn=righthorn, rms=rms, quiet=quiet)#

Function to find the area, width, and velocity of a galaxy profile.

There are four ways in which the left and right (low channel number and high channel number) edges of the galaxy profile can be determined. The mode parameter selects one of these methods:

As a fraction of the mean within the region of interest. The mean of data from brange through erange is calculated. The edges are then those locations where the data values are greater than fract*mean for 3 consecutive channels starting from the end points of the region of interest and searching towards the center.
As a fraction of the maximum value within the region of interest. The peak of data from brange to erange is found. The edges are then those locations where the data values are greater than fract*(peak-rms) for 3 consecutive channels starting from the end points of the region of interest and searching towards the center.
As a fraction of each of two peaks - identified by the user. The user uses the cursor to mark two peaks in the region of interest or those peaks are identified through the lefthorn and righthorn parameters. The maximum value within 10 channels of each user-supplied peak location is found. The left edge is where the data value falls below fract*(peak-rms) for 3 consecutive channels searched from the location of the peak. The right-channel peak is similarly used to find the right edge.
A polynomial is fit to either side of the profile between 15-85% of (peak-rms) (where the peak is from the relevent side of the profile - the two peaks are identified by the user). The velocity of the polynomical fit at fract*(peak-rms) is then found for both sides of the profile. The difference between the two is the width, the mean is the central velocity. The order of the polynomial is chosen by the user (linear fits are usually sufficient). This scheme (further described in Springob etal 2006) has the advatange of averaging out noise effects from either side of the profile.

In the first 3 modes, the final left and right edge are linear interpolations to get the fractional channel where data value crossed the threshold given by fract for that particular mode.

If an edge is not found, a warning is issued and the appropiate end-point of the region of interest is used. Only the data values within the region of interest are used here.

The returned value is a 6-element array with these values, in this order.

Area	The sum of data[i]*abs(vel[i+1]-vel[i-1])/2.0 for all i from brange to erange (modes 1 and 2) or between the channels where the data values first become negative moving out from the two peaks found in mode 3 (not including that transition channel).
Width	The absolute value of the difference between the left and right edges as determined for that mode.
Velocity	The mean of the left and right edges as determined for that mode.
Error on Area	Error estimate on the Area. The error on Area is only set for the fract=0.5 and fract=0.2 cases. Unset error estimates have a value of 0.0
Error on Width	Error estimate on the Width. The Width error is only set by mode 4. For mode 4, the Width error is always set at 2x the Velocity error. Unset error estimates have a value of 0.0
Error on Velocity	Error estimate on the Velocity. The Velocity error is only set by mode 4. Unset error estimates have a value of 0.0

Blanked data is ignored by this routine. Since the velocities used to calculate widths and centers come from the centers of the valid channels, ignoring the blanked data is equivalent to replacing the blanked data by it’s nearest non-blanked neighbor in the direction of the edge searches. If the data is all blanked in the region of interest, the returned values are all 0.

This code adapted from code in use at Arecibo. This particular version was originally from Karen O’Neil. Option 4 added by Karen Masters as part of work for GBT06A-027, GBT06B-021, GBT06C-049, GBT08B-003: “Mapping Mass in the Nearby Universe with 2MASS”.

Contributed By: Karen O’Neil, Bob Garwood, and Karen Masters.

Params:

datain, required, type=float array: The data values.
velin, required, type=float array: The velocities at each data point, in km/s.
brangein, required, type=integer: The first channel to use.
erangein, required, type=integer: The last channel to use.
modein, required, type=integer: The method to use in finding the returned values.
fractin, required, type=float: Used in locating the edges of the galaxy profile. See the documentation for more details.

Keywords:

lefthornin, optional, type=float: The location (in channels) of the left (low channel number) peak in the profile. Used only in mode 3. If this or righthorn are not provided, the user is asked to use the cursor to mark these locations.
righthornin, optional, type=float: The location (in channels) of the right (high channel number) peak in the profile. Used only in mode 3. If this or lefthorn are not provided, the user is asked to use the cursor to mark these locations.
rmsin, optional, type=float: Used in modes 2 and 3 as described above. If this is not supplied, defaults to the stddev of data within the region of interest.
quietin, optional, type=boolean: When set, the results are not printed to the terminal (they are still returned).

Returns:

the values in a 3 element array: [0] is the area, [1] is the width and [2] is the velocity. Returns 0.0 for all 3 values on error.

function chantofreq(data, chans, frame=frame, true_frame=true_frame)#

Convert channel number to frequency (Hz) using the supplied data container. Optionally return the frequency in another reference frame from that given in the data.

Params:

datain, required, type=spectrum: The spectrum data container to use to get the necessary header information.
chansin, required: The channel numbers to convert, may be an array of values.

Keywords:

framein, optional, type=string: The rest frame to convert to. Known rest frames are listed in frame_velocity. Defaults to the frame given in the data.frequency_type.
true_frameout, optional, type=string: The actual rest frame used in constructing the frequencies. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as the frame in data.frequency_type.

Returns:

frequencies (Hz)

Uses:

DATA_VALID frame_velocity

function chantovel(data, chans, frame=frame, veldef=veldef, true_frame=true_frame, true_veldef=true_veldef)#

Convert channel number to velcity (m/s) using the supplied data container. Optionally return the velocity in another reference frame.

“Params:

datain, required, type=spectrum: The spectrum data container to use to get the necessary header information.
chansin, required: The channel numbers to convert, may be an array of values.

Keywords:

framein, optional, type=string: The rest frame to convert to. Known rest frames are listed in frame_velocity. Defaults to the frame given in data.velocity_definition.
veldefin, optional, type=string: The velocity definition to use from RADIO, OPTICAL, or TRUE. Defaults to value found in data.velocity_definition.
true_frameout, optional, type=string: The actual rest frame used in constructing the velocity axis. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as the frame in data.velocity_definition.
true_veldefout, optional, type=string: The actual velocity frame used in constructing the velocity axis. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as data.velocity_definition.

Returns:

velocity (m/s)

Uses:

DATA_VALID chantofreq freqtovel

function chdoppler(ra, dec, julday, obspos=obspos, light=light)#

Computes the projected velocity of the telescope wrt six coordinate systems: geo, helio, bary, lsrk, lsrd, gal negative velocities mean approach

The standard LSR is defined as follows: the sun moves at 20.0 km/s toward ra=18.0h, dec=30.0 deg in 1900 epoch coords

Fully vectorized. All three arguments must have the same dimensions.

This code came via e-mail from Carl Heiles via Tom Bania on 11/04/2004. Updated using code found via google from same source on 11/30/2009. Local changes:

modify this documentation for use by idldoc.
removed path argument and related code, replaced by obspos argument.
default position is the GBT.
Observatory longitude was not being passed to juldaytolmst.
LSRD added
Galactocentric added
Checked against aips++ Measures. Differences were less then 20 m/s in one test case (less then 10m/s for geo, bary, and lsrk).
Double precision throughout.
Relativistic addition of velocities.

Previous revision history: carlh 29oct04

from idpppler_chl; changed calculation epoch to 2000.
19nov04: correct bad earth spin calculation
07 Jun 2005: vectorize to make faster for quantity calculations
20 Mar 2007: CH updated documentation
08 Feb 2015: CH fixed doppler additions for array inputs. See

annotated statements at end of program.

Params:

ra {in}{required} The source ra in decimal hours, equinox 2000 dec {in}{required} The source dec in decimal hours, equinox 2000 julday {in}{required} The julian day

Keywords:

obspos {in}{optional}{type=double [2]}: observatory position [East longitude, latitude] in degrees. Uses the GBT position if not specified.
lightin, optional, type=boolean: When set, returns the velocity as a fraction of c

Returns:

The velocity in km/s, or as a faction of c if the keyword /light is specified. the result is a 6-element vector whose elements are [geo, helio, bary, lsrk, lsrd, gal].

Uses:

baryvel precess juldaytolmst

function compress_ints(ints, llimit=llimit, ulimit=ulimit)#

Take an integer array, sort it, and convert it to a comma separated string where consecutive values of integers are ‘compressed’ into a range, using the syntax ‘begining:end’.

The returned string is collection of individual integers and ranges separated by commas.

Two consecutive values are left as individual integers, only sequences longer than 3 integers are compressed to a range (there’s no savings in space for 2 consecutive integers.

Optional upper and lower limits to the integer values can be specified. Values less than the lower limit are set to that limit and values greater than the upper limit are set to that limit.

Duplicate values are removed before the array is compressed to the string.

Params:

intsin, required, type=long: Integer array to convert

Keywords:

llimitin, optional, type=long: Optional lower limit to use.
ulimitin, optional, type=long: Optional upper limit to use.

Returns:

a string represent all of the unique values in the array passed in, optionally truncated to fall withing llimit to ulimit, inclusive.

function coord_interp(coords, loc, wrapsat)#

Find the value of the coordinate at the given location, watching for large discontinuities due to a coordinate wrap (e.g. time at 24/0 or degrees at 360/0).

Params:

coordsin, required, type=2-element array: The values to use in the interpolation. Only the first two values are used.
locin, required: The location to find the interpolated value at. Should be between 0 and 1.
wrapsatin, required: The value that coords wraps back to zero at.

Returns:

The interpolated coordinate value.

function count_blanks(dc, size)#

Count the number of blanked values in the data array of the given data container.

A blanked value is a NaN, this returns the number of NaNs in the given data container. It returns -1 if the dc argument is not a valid data container. The size argument will contain the total number of data elements in dc.

Params:

dcin, required, type=data container: The data container to check.
sizeout, optional, type=integer: The total number of data elements in dc (-1 if dc is an invalid data container)

Returns:

The total number of blanked values in dc. Returns -1 if dc is not a valid data container.

pro DATA_COPY(in, out)#

Copy the data container from in to out.

This procedure can be used to generate a new data container, or to copy one data container into an existing data container. However, it cannot be used to copy into one of the global data containers. To copy a data container stored as a local variable into a global data container, use the procedure set_data_container.

Params:

in ; in, required, type=data_container_struct: The data container copied. This can be identified by a local valiable or it can be a global data container, such as !g.s[0]
outout, required, type=data_container_struct: The data container to receive the copy. This can be a local variable, but NOT a global data container.

Examples:

get, mc_scan=22, cal='F', sig='T', pol='XX', if_num=1, int=1 ; sig
data_copy,!g.s[0],spec
a = getdcdata(spec)
a = a * 2.0
setdcdata,spec,a
show,spec
data_free, spec  ; clean up memory

Uses:

DATA_VALID DATA_FREE

pro DATA_FREE(data_struct)#

Free the data pointer in a data structure. This should only be used when the data structure is no longer necessary since it leaves the data pointer in an invalid state.

Params:

data_structin, out, required, type=data_container_struct: The struct to free.

function DATA_NEW(arr, spectrum=spectrum, continuum=continuum, nocheck=nocheck)#

Create a new data_struct of the requested data type (spectrum or continuum) and containing the given array, if supplied. For continuum data, when the data array is supplied, the additional pointers (date, utc, mjd, etc) are set to double precision vectors of 0s having the same number of elements as arr.

Params:

arrin, optional, type=integer, default=undefined: The data pointer that this data_struct will hold. If arr is not provided, then the data pointer will point to an undefined variable.

Keywords:

spectrumin, optional, default=set: When this is set, a spectrum structure will be returned. That is the default behavior.
continuumin, optional, default=unset: When this is set, a continuum structure will be returned. spectrum and continuum are mutually exclusive.
nocheckin, optional, default=unset: When this is set, the input parameter checking is turned off. Usefull for speed.

Returns:

requested data structure of given size or -1 on failure.

function DATA_VALID(data_struct, name=name)#

Check to see if the argument is a recognized data structure.

A data structure is invalid if its name is not ‘SPECTRUM_STRUCT’ or ‘CONTINUUM_STRUCT’.

Params:

data_structin, required, type=data structure: The structure check.

Keywords:

nameout, optional, type=string: Holds the structure name when used in the call. This will be either ‘CONTINUUM_STRUCT’ or ‘SPECTRUM_STRUCT’.

Returns:

n_elements(*data_ptr) if valid and -1 if invalid. (hence if the return value is = 0 it’s undefined).

function daynotodm(dayno, year)#

Convert daynumber and year to day of month (1..31) and month of year (1.l12).

This code came from Phil Perillat at Arecibo.

Local changes:

modify this documentation for use by idldoc.

Params:

daynoin, required, type=long integer: daynumber of year 1..365 or 366
yearin, required, type=long integer: 4 digit year

Returns:

[day,month] as a vector.

pro dcaccum(accumbuf, dc, weight=weight, quiet=quiet)#

Add a data container into an ongoing accumulation in a given accum_struct structure.

If this is the first item to be “accum”ed, it will also be used as a template, stored in accumbuf.template, to be used by accumave.

If the polarization of any items being accumed does not match that of template, the polarization of the template is changed to ‘I’.

This combines the UniPOPS functionality of ACCUM and SUM. The SUM name is already in use in IDL.

This is primarily for use inside procedures and functions where it is useful to average several data containers without disturbing the public data containers in the guide structure. Most users will find using accum preferable to using dcaccum.

The data are : sum(weight*data)
The times are : sum(duration), sum(exposure)
The weight is : sum(weight) one weight sum per channel.
The tsys is : sqrt(sum(max(weight)*Tsys^2))
The frequency resolution is the maximum of all f_res values used during the accumulation.

A warning message is shown if either the frequency_resolution or the frequency_interval do not match that in an already on-going accumulation. If the quiet flag is on, then this message is suppressed. In either case, the accumulation proceeds.

If a weight is not supplied, it will be exposure*frequency_resolution/tsys^2.

Weight can either be a scalar or it can be a vector having the same number of elements as the data in dc.

If all of data is blanked (not a number) then it is completely ignored and the accumulated weight, times, and system temperatures are unchanged. If individual regions are blanked then the weight at those channels is 0. When an average is requested (accumave) this weight array is used to rescale the data. That weight array is also available when the average is requested. If that weight array is used as input in a future average, the averaging can continue from the same point as before.

If all of the weight values (supplied or the default value as described above) are not finite (not a number) this routine behaves as if the data are blanked and the data are ignored and the accumulated values are unchanged.

If the Tsys value is not finite (not a number) but the weights and data are at least partially finite, then the Tsys value is ignored in the ongoing weighted Tsys^2 accumulation. The weights used in the Tsys^2 accumulation are kept separate from the data weights.

Params:

accumbufin, out, required, type=accum_struct: The structure containing the accumulation that you want to add to.
dcin, required, type=spectrum: The data container to accum.

Keywords:

quietin, optional, type=boolean: If set, suppress warning messages about frequency resolution and interval not matching values in accumbuf.
weightin, optional, type=float: The weight to use for this data. If this is not set, a weight of exposure*frequency_resolution/tsys^2 will be used. Weight can either be a scalar (uniform weight across at all channels) or it can be an array having the same number of elements as the data in dc.

Examples:

average some data

a = {accum_struct}
accumclear,a  ; not necessary here, but a good habit to follow
; get several records at once
s = !g.lineoutio->get_spectra(index=0)
dcaccum,a,s
data_free,s ; be sure to clean up, else leaks memory
s = !g.lineoutio->get_spectra(index=1)
dv = dcvshift, a, s ; align in velocity
dcshift, a, dv ; actually do the shift to align
dcaccum,a,s
data_free, s
accumave,a,s
show, s
data_free, s  

See ave for additional examples.

Uses:: accum DATA_VALID DATA_FREE

function dcadd(dc1, dc2)#

This procedure adds two data container’s data arrays, returning the sum

Params:

dc1in, required, type=data container: data container (spectrum or continuum)
dc2in, required, type=data container: data container (spectrum or continuum)

Returns:

float array - sum of two data containers’ data

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show,a
get, index=2
b = data_new()
data_copy,!g.s[0],b
show,b
sum = dcadd(a,b)
plot, sum

Uses:

dcpaircheck

pro dcascii(dc, brange=brange, erange=erange, file=file)#

Output a data container as ascii showing the x-axis value and the corresponding y-axis value over the given range of X values.

Params:

dcin, optional, type=data container: The data container to use. If not supplied, the plotter is queried for the currently displayed values.

Keywords:

brangein, optional, type=float, default=all: beginning of the range to plot, in current plotter units
erangein, optional, type=integer, default=all: end of the range to plot, in current plotter units
filein, optional, type=string, default=/dev/tty: The file to write to. Defaults to the current screen, using “more” to page the output.

Examples:

list the current displayed spectrum at the terminal from 1415 to 1425 MHz. The current x-axis must be in MHz.

dcascii,brange=1415,erange=1425

the output looks something like this:

Scan:     88         GalPlane 2005-11-15 +08 06 49.0
                           Ta    
       MHz-LSR                 YY
1424.9752215054413682     -0.0136188
1424.9361597211752724      0.0389956
1424.8970979369089491      0.0697150
1424.8580361526428533      0.0732256
1424.8189743683767574      0.0983463
1424.7799125841104342      0.1009806
.....
< Press Spacebar to continue, ? for help >

The first line shows the scan number, source name, date, and UT associated with this data. The next two lines describe what is in each column. The x-axis values are in the left column and the y-axis values are in the right column. The first header in the left column gives the velocity definition if velocities are show (otherwise it is empty). The second header in that column gives the units and the reference frame. The first header in the right column gives the units and the second header gives the polarization.

function dcbaseline(dc, nfit, regions, nregion, polyfit, polyrms)#

Fit a baseline to the contents of the dc argument with various options.

This fits a set of orthogonal polynomials with a maximum order given by the nfit argument to the data values in dc. The range of channels to use is given by the regions and nregion arguments. The parameters of the fit are placed in the value of the polyfit argument and the rms values of the fits are placed in the polyrms argument. If there is a problem, the return value is 0, otherwise the return value is 1.

Params:

dcin, required, type=data container: The data container to use. The values to fit come from *dc.data_ptr.
nfitin, required, type=integer: The order of polynomial to fit.
regionsin, required, type=2D array: Array describing the regions to use in the fit, along with nregion. This has the same properties as the !g.region field in the guide structure.
nregionin, required, type=integer: The number of regions in regions to actually use.
polyfitout, type=double array: Parameters describing the polynomials. Array with (nfit+1) elements. See ortho_fit for more details.
polyrmsout, type=double array: RMS values, one for each polynomial (nfit+1).

Returns:

1 on success, 0 on failure.

Examples:

ok = dcbaseline(!g.s[0], 2, !g.regions, !g.nregion, polyfit, polyrms)

Uses:

DATA_VALID get_chans ortho_fit

pro dcbias(dc, factor)#

This procedure adds a bias factor to a data container’s data

Params:

dcin, required, type=data container: data container (spectrum or continuum)
factorin, required, type=float: bias factor

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show
dcbias,a,25
show,a

pro dcboxcar(dc, width, decimate=decimate)#

Smooth a data container with a boxcar smoothing of a certain width, in channels. For odd width, this uses the built-in idl SMOOTH function. For even widths this uses doboxcar1d and the reference channel is moved left by 1/2 channel width.

Replaces the contents of the data being smoothed with the smoothed data. Use the boxcar procedure to smoothing data containers in the !g structure.

For spectrum data containers, the frequency_resolution is set using estboxres.

Params:

dcin, required, type=data container: The data container to be smoothed. The data values are modified by this procedure.
widthin, required, type=integer: Width of boxcar in channels.

Keywords:

decimatein, optional, type=boolean: If set, the data container is reduced - taking every width channels starting at channel 0.

Uses:

doboxcar1d

pro dcclip(dc, datamin, datamax, blank=blank)#

This procedure truncates (clips) data above and below given intensity limits. Data value can alternatively be blanked using the /blank flag. In that case, data values outside the limits are replaced by blanks (NaN).

Params:

dcin, required, type=data container: data container (spectrum or continuum)
dataminin, required, type=float: min data value
datamaxin, required, type=float: max data value

Keywords:

blankin, optional, type=boolean: Replace clipped values with NaN instead of the clipping limit.

Examples:

getrec,1
a = data_new()
data_copy,!g.s[0],a
show
dcclip,a,-1.0,2.0
show,a

pro dcconvol(dc, kernel, scale_factor, ok=ok, normalize=normalize, _extra=extra_keywords)#

A front-end to the standard IDL CONVOL function.

This allows users to convolve a data container with an arbitrary kernel using the IDL CONVOL function. All of the arguments except the dc parameter have the same meaning and use as in the CONVOL function. Users should consult the IDL documentation for CONVOL for a more detailed explanation than is given here. All of the CONVOL keywords are passed in through the _EXTRA keyword are so are not shown explicitly here. These keywords are /CENTER, /EDGE_WRAP, /EDGE_TRUNCATE, MISSING, /NAN and /NORMALIZE.

This procedure replaces the data values in dc with the results of the convolution.

NORMALIZE was added in IDL 6.2. It is essential when working with blanked data (NAN). It is implemented here so that it works for earlier version of IDL. For IDL 6.2 and later, /NORMALIZE is passed directly to CONVOL. For earlier versions of IDL, a second convolution is done using a vector of 1.0s having the same length as the and blanked in the same locations as the data. The convolution of the data is then divided by this convolution and the result is put back into the data container. BIAS was also added in IDL 6.2 but that functionality isn’t necessary for the rest of GBTIDL and so it has not been implemented here appart from what CONVOL provides in IDL 6.2

Params:

dcin, required, type=data container: The data container to use in the convolution.
kernelin, required, type=array: The kernel to use in the convolution. Must have fewer elements than the data container being convolved.
scale_factorin, optional, type=real, default=1: The scale factor.

Keywords:

okout, optional, type=boolean: Returns 1 if everything went ok, 0 if it did not (missing parameter, empty or invalid dc, bad kernel).
normalizein, optional, type=boolean: Set this keyword to automatically compute a scale factor and bias and apply it to the result values. If this keyword is set, the scale_factor argument and the BIAS keyword are ignored. For all input types, the scale factor is defined as the sum of the absolute values of Kernel. If blanked (NAN) values are present, the scale factor and bias are calculated without using those values so that all result values are comparable in magnitude.
_extrain, optional, type=extra keywords: Keyword arguments to CONVOL

Examples:

; hanning smoothing kernel
kernel = [0.25,0.5,0.25]
; dc is some data container that already exists and is valid
dcconvol, dc, kernel
; same kernel,ignore NAN (missing) values, truncate the 
; data at the edges and normalize.  This is how the dchanning
; procedure is implemented.
dcconvol, dc, kernel, /nan, /edge_truncate, /normalize

Uses:

setdcdata

pro dcdecimate(dc, nchan, startat=startat, ok=ok)#

This procedure decimates a spectrum (i.e. thins the spectrum by paring every nth channel). The dc argument is modified in place.

The frequency interval and reference channel are adjusted appropriately so that the x-axis labels for the decimated data are still appropriate.

This only works for spectral line data containers.

Params:

dcin, out, required, type=data container: data container (spectrum or continuum)
nchanin, required, type=integer: choose every nth channel starting at startat.

Keywords:

startatin, optional, type=integer, default=0: The starting channel.
okout, optional, type=boolean: Returns 1 if everything went ok, 0 if it did not (missing parameters, invalid or empty dc, bad startat)

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show
dchanning,a
show,a
dcdecimate,a,3
show,a
dcdecimate,a,3,startat=1
show,a

:Uses”: DATA_VALID seq

function dcdivide(dc1, dc2)#

This procedure divides the first data container’s data array by the second, returning the result.

Params:

dc1in, required, type=data container: data container (spectrum or continuum)
dc2in, required, type=data container: data container (spectrum or continuum)

Returns:

float array - division of dc1 data array by dc2 data array

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show,a
get, index=2
b = data_new()
data_copy,!g.s[0],b
show,b
div = dcdivide(a,b)
plot,div

Uses:

dcpaircheck

pro dcdofeedoffsets(dc, remove=remove, reportonly=reportonly, status=status)#

Apply or remove the feed offsets in the supplied data container.

Note: The feed offsets need only be applied for old (fitsver earlier than 1.7) data written by gbtidl. Data produced by sdfits from fitsver 1.7 on has had the feed offsets applied by sdfits. In addition, GBTIDL applies all feed offsets by default when reading from SDFITS files having an SDFITSVER keyword (files produced by sdfits). Other old SDFITS file are likely produced by GBTIDL using KEEP or SAVE and it’s impossible to tell whether any non-zero feed offsets found there have been applied to the positions and so GBTIDL does not automatically apply the feed offsets when reading that data.

This applies (default) or removes (when the /remove keyword is set) any non-zero feed offset values (feedxoff and feedeoff) found in the supplied data container. If /reportonly is used then information is printed about what the new coordinate values would have been but the data container remains unchanged.

Feed offsets can only be applied if the longitude and latitude coordinates can be converted to and from Azimuth and Elevation. Recognized coordinate systems are “RADEC”,”HADEC”,”GALACTIC”,and “AZEL”. If the offsets can not be applied (or removed) the status value will be set to 0. When successfull, the status value is 1.

This works with spectral line and continuum data containers.

Params:

dcin, required, type=data container: The data container to be used. The longitude, latitude, azimuth, and elevation values are changed by this procedure unless /reportonly is set.

Keywords:

removein, optional, type=boolean, default=0: When set, remove the feed offsets from the positions. The default is to apply the feed offsets.
reportonlyin, optional, type=boolean, default=0: When set, generate a report only. No values in dc are changed when that is set. For continuum data containers, the report is generated for a position near the middle integration.
statusout, optional, type=integer: On success, this is set to 1. On failure it is set to 0.

function dcextract(dc, startat, endat, stride)#

Extract a region from a data container, producing a new data container with fewer elements in it. The caller is responsible for eventually freeing the pointer(s) contained in the returned data container using data_free. The only header words modified are the reference_channel and frequency_interval for spectral line data. If only the dc is supplied, this function is equivalent to data_copy.

Params:

dcin, required, type=data container: The data container to extract the region from.
startatin, optional, type=integer, default=0: The first element to include in the extracted region.
endatin, optional, type=integer: The last element to include in the extracted region. If not supplied, the last element is used. endat must be greater than or equal to startat.
stridein, optional, type=integer, default=1: The increment in elements to extract, starting with startat. stride must be greater than 0.

Returns:

The extracted data container. The user is responsible for eventually freeing this using data_free. Returns -1 on error.

function dcfft(real, imag, inverse=inverse, bdrop=bdrop, edrop=edrop)#

Do an FFT (forward or inverse) of a data container. For a real-to-complex FFT, only one data container is given in the arguments. For a complex-to-complex FFT, the real part comes from the first data container and the imaginary part comes from the second data container. The returned value is always a complex array containing the results of the FFT. The input data containers are not changed by this function. It is up to the caller to determine what to do with this result (e.g. further calculations, store the real part in one data container and the imaginary part in another, etc).

This uses the builtin IDL FFT function.

For spectral-line data, when inverse is not set (the default) the builtin IDL inverse flag is used. This means that when dcfft is called with its default arguments and the input data container is a spectrum, the IDL FFT will convert properly from the frequency domain to the time domain (which is actually an inverse FFT). When inverse is set, the builtin IDL inverse flag is not set. This may be confusing to IDL users, but it will be familiar to former UniPOPS users.

For continuum data, the inverse flag here is exactly the same as the inverse flag in the builtin IDL function.

Params:

realin, required, type=data container: The real part of the data to be FFTed.
imagin, optional, type=data container: The imaginary part of the data to be FFTed. When not supplied, it is assumed that the data are pure real (imaginary part is all zero).

Keywords:

inversein, optional, type=boolean: When set, the inverse of the regular FFT is done. For spectral-line data, an inverse FFT is done when this is not set and a direct FFT is done when this is set so that the non-inverse transformation as seen by the dcfft user is frequency to time and the inverse transformation is time to frequency.
bdropin, optional, type=integer, default=0: The number of channels to exclude from the FFT at the beginning.
edropin, optional, type=integer, default=0: The number of channels to exclude from the FFT at the end.

Returns:

A complex array containing the result of the FFT. Returns -1 if no data was found in the real argument.

function dcfold(sig, ref, ftol=ftol, blankinterp=blankinterp, nomask=nomask)#

Average two parts of in-band frequency switching (signal and reference phases, with the cal-switching phases already calibrated). Typically this happens after getfs. It does not matter which data container containers which part since their relative distance in frequency is used to determine how one is shifted to align with the other. The returned result is a new data container that the user must eventually free using data_free in order to avoid any memory leaks. sig and ref must have the same number of channels and the same spacing between channels. If sig and ref are already aligned, then this is a simple average (a warning message will be printed). If the shift necessary to align ref with sig is more than the total number of channels, there is no overlap and this procedure will print an error message and return without altering sig and ref.

The “ref” data container is shifted to align in sky frequency with the “sig” data container using dcshift and the two data containers are averaged - weighting each by the inverse of square of their system temperatures. The system temperature in the result is the weighted average of the two system temperatures.

If there are any blanked channels in either “sig” or “ref” then the corresponding channels in the other spectrum, after “ref” has been shifted to align in sky frequency with “sig”, are also blanked so that the average at those channels is blanked. This is done to avoid the appearance of a spike at a blanked channel where the contribution from the other spectrum, after the shift, is typically not blanked and not equal to the local average of surrounding, non-blanked, channels. This behavior can be turned off by the nomask keyword. Alternatively, all blanked channels in “sig” and “ref” can be replaced using a linear interpolation from adjacent non-blanked channels (a simple average in the case of single, isolated, blanked channels) using the blankinterp keyword.

Params:

sigin, out, required, type=spectrum: The data container to use as the signal part in the average.
refin, out, required, type=spectrum: The data container to use as the reference part in the average. This data is shifted using dcshift to align with “sig” before averaging.

Keywords:

ftolin, optional, type=double, default=0.005: The fractional channel shift tolerance. If the fractional part of the channel shift necessary to align the two parts is less than this value, no fractional shift as described in the documentation for dcshift will be done. It might be useful to turn off the fractional shift because it can cause aliases and ringing in the case of very strong lines or other sharp features. If ftol > 0.5 no fractional shifts will be done.
blankinterpin, optional, type=boolean: When set, blanks are replaced before shifting and averaging by a linear interpolation using the finite values found in the two spectra. The dcinterp procedure is used. For single blanked channels, the replacement value is the average of the two adjacent channel values.
nomaskin, optional, type=boolean: When set, turn off the masking of blank channels from each spectrum on to the other, after the shift. This may result in spikes at the location of blanked channels. This was the original behavior of this routine. This keyword has no effect if blankinterp is set.

Returns:

data container. The user is responsible for freeing this. returns -1 on error.

Uses:

dcshift dcinterp data_free data_valid

function dcfshift(accumbuf, dc, frame=frame)#

Function to calculate the shift, in channels, necessary to align in frequency a data container with the data container template in an ongoing accumulation.

If the frame is not set, the one implied by the data header is used. Use dcxshift to align using the current settings of the plotter’s x-axis.

Params:

accumbufin, required, type=accum_struct: The ongoing accumulation buffer.
dcin, required, type=spectrum: The data container that needs to be shifted to align with the data container template in accumbuf.

Keywords:

framein, optional, type=string: The reference frame to use. If not supplied, the value implied by the last 4 characters of the velocity_definition in the ongoing accumulation will be used. See frame_velocity for a full list of supported reference frames.

Returns:

shift, in channels, to be used as argument to dcshift. Returns 0.0 on error.

Examples:

a={accum_struct}
getps,30
dcaccum, a, !g.s[0]  ; start an accum, no alignment needed yet
getps,31
fs = dcfshift(a,!g.s[0]) ; what is the shift to align 31 with 30?
; get a copy of data at !g.s[0]
data_copy,!g.s[0], d
dcshift, d, fs  ; actually shift the data
dcaccum, a, d ; and add it in
getps, 32
data_copy,!g.s[0], d
dcshift, d, dcfshift(a, d)  ; all in one line, shift 32 to align with 30
dcaccum, a, d
accumave, a, d ; result is in d now

pro dchanning(dc, decimate=decimate, ok=ok)#

This procedure smooths a spectrum with a hanning filter.

Replaces the contents of the data being smoothed with the smoothed data. Blanked data values are ignored.

For spectrum data containers, the frequency_resolution is set using esthanres.

Params:

dcin, required, type=data container: data container (spectrum or continuum)

Keywords:

decimatein, optional, type=keyword: decimate the spectrum?
okout, optional, type=boolean: Returns 1 if everything went ok, 0 if it did not (missing dc parameter, invalid or empty dc)

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show
dchanning,a,/decimate
show,a

Uses:

dcconvol dcdecimate

pro dcinterp(dc, bchan=bchan, echan=echan, linear=linear, quadratic=quadratic, lsquadratic=lsquadratic, spline=spline, quiet=quiet, ok=ok)#

Interpolate across blanked channels in a data container.

This uses the IDL INTERPOL function to replace blanked values in the data container with unblanked values according to the interpolation method selected.

You can limit the range of channels to consider using bchan and echan. When not supplied, all of the channels are used.

The default interpolation method is linear. The other interpolations may not be particularly useful across large gaps.

If all of the data in the requested range is good (unblanked) then no interpolation is done and this routine silently returns without changing anything in dc.

It is an error to request more than one interpolation method.

Params:

dcin, required, type=data container: The data container to use.

Keywords:

bchanin, optional, type=integer: The starting channel number. If not specified, bchan=0.
echanin, optional, type=integer: The last channel number. If not specified use all channels from bchan to the end.
linearin, optional, type=boolean: When set, use the linear interpolation provided by INTERPOL. This is the default interpolation when no other method is specified.
quadraticin, optional, type=boolean: When set, use the quadratic interpolation provided by INTERPOL.
lsquadraticin, optional, type=boolean: When set, use the lsquadratic (lest squares quadratic) interpolation provided by INTERPOL.
splinein, optional, type=boolean: When set, use the spline interpolation provided by INTERPOL.
quietin, optional, type=boolean: When set, do not print out informational messages. Useful primarily for supressing the message that happens when an attempt is made to interpolate across a spectrum of all blanked values. This does not supress errors due to argument errors.
okout, optional, type=boolean: This is set to 1 on success or 0 on failure (e.g. bad arguments).

Uses:

data_valid

pro dcinvert(dc)#

Flip the data end-to-end in the supplied data container.

For line data the value of frequency_increment and reference_channel are also changed appropriately so that, as displayed, there will be no change in appearance. This is useful if you need to combine (e.g. average) two data containers where the frequency increments have opposite signs.

For continuum data (where the need to invert is less obvious), all of the time-dependent arrays are also flipped (utc, mjd, etc).

The invert is done in place.

This procedure can not be used on the global guide data containers because they are passed by value.

Params:

dcin, required, type=data container: The data container to be smoothed. The data values are modified by this procedure.

function dcmeantsys(dc_nocal, dc_withcal, tcal=tcal, used_tcal=used_tcal)#

Calculate the mean Tsys using the data from two spectral line data containters, one with the CAL on and one with CAL off.

\[mean_tsys = tcal * mean(nocal) / (mean(withcal-nocal)) + tcal/2.0\]

where nocal and withcal are the data values from dc_nocal and dc_withcal and tcal is as described below.

The outer 10% of all channels in both data containers are ignored.
Blanked data values are ignored.
The tcal value used here comes from the dc_nocal data container unless the user supplies a value in the tcal keyword.
The tcal value actually used is returned in used_tcal.

This is used by the GUIDE calibration routines and is encapsulated here to ensure consistency.

Params:

dc_nocalin, required, type=spectrum data container: The data with no cal signal.
dc_withcalin, required, type=spectrum data container: The data with a cal signal.

Keywords:

tcalin, optional, type=float: A scalar value for the cal temperature (K). If not supplied. dc_nocal.mean_tcal will be used.
used_tcalout, optional, type=float: The tcal value actually used.

pro dcmediansub(dc, width, ok=ok)#

Subtract the median filtered values of the given width, in channels, from the data. The result replaces the original data values.

Uses the IDL MEDIAN function to get the median filtered array.

Params:

dcin, required, type=data container: The data container to smooth.
widthin, required, type=integer: The desired number of channels to use in performing the median filter.

Keywords:

okout, optional, type=boolean: Returns 1 if everything went ok, 0 if it did not (missing parameter, empty or invalid dc, bad width).

Examples:

; dc already exists and is a valid data container
; subtract a median filter of width 200
dcmediansub,dc,200

:Uses”: data_valid

function dcmultiply(dc1, dc2)#

This procedure multiplies two data container’s data arrays, returning the product

Params:

dc1in, required, type=data container: data container (spectrum or continuum)
dc2in, required, type=data container: data container (spectrum or continuum)

Returns:

float array - product of dc1 and dc2’s data arrays

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show,a
get, index=2
b = data_new()
data_copy,!g.s[0],b
show,b
product = dcmultiply(a,b)
plot, product

Uses:: dcpaircheck

function dcpaircheck(dc1, dc2, msg)#

This procedure checks both data containers to see if their data arrays can be used for mathematical functions, such as dcadd, dcsubtract, etc.

Params:

dc1in, required, type=data container: data container (spectrum or continuum)
dc2in, required, type=data container: data container (spectrum or continuum)
msgout, required, type=string: string containing error if found

Returns:

0 - pair is incompatible; 1 - pair is compatible

Examples:

get, index=1
a = data_new()
data_copy,!g.s[0],a
show,a
get, index=2
b = data_new()
data_copy,!g.s[0],b
show,b
status = dcpaircheck(a,b,msg)
if (status ne 1) then print, msg 
"data containers must contain data of equal length" 

pro dcpowspec(dc)#

Compute the power spectrum for the supplied data container, which is overwritten with the result. Uses dcfft to do the fft. This is simply the sum of squares of the real and imaginary parts of the fft on the data.

Params:

dcin, out, required, type=data container: The data container to use for both input and output.

Uses:

dcfft

pro dcresample(dc, newinterval, keychan, nearest=nearest, linear=linear, lsquadratic=lsquadratic, quadratic=quadratic, spline=spline, ok=ok)#

Resample a spectrum using one of 4 possible interpolation methods onto a new frequency axis.

Takes the data values in the given data container and interpolates them onto a new frequency axis derived from the newinterval and keychan arguments. The reference_frequency is unchanged by this operation. The new frequency_interval is given by the newinterval argument. The new frequency axis is determined by choosing one channel in the original frequency axis and requiring that the frequency at that channel be centered on one of the channels in the interpolated data (the keychan argument). This sets the new reference frequency. When not supplied, keychan defaults to the channel nearest to the original reference channel.

The dc argument is altered by this operation. The data values will contain the interpolated values and the header will describe the new frequency axis. The number of channels in the resulting interpolated data is just enough to use as much of the original data as possible without trying to interpolate beyond the ends of that data.

Interpolation is done using INTERPOL or internally (nearest). The channel number are used as the abscissa. Blanked values (NaN) are ignored. See the documentation for INTERPOL for details about the various interpolation methods.

This only works with spectrum data containers.

Params:

dcin, out, required, type=data container: The data container to work on.
newintervalin, required, type=real: The new frequency_interval to use in the interpolation. This must be non-zero. If it has opposite sign from the input frequency interval, dcinvert will be used to first reverse the sense of the frequency axis.
keychanin, optional, type=integer: The new frequency axis will have one channel where the frequency value equals the original frequency value at the keychan channel. When not supplied, this defaults to the channel nearest to the reference_channel.

Keywords:

nearestin, optional, type=boolean: When set do nearest-neighbor interpolation.
linearin, optional, type=boolean: When set (the default) do a linear interpolation.
lsquadraticin, optional, type=boolean: When set do least squares quadratic fit interpolation.
quadraticin, optional, type=boolean: When set do a quadratic fit interpolation.
splinein, optional, type=boolean: When set do a spline interpolation.
okout, optional, type=boolean: This is set to 1 on success, otherwise it is 0 (bad arguments, bad data).

Examples:

; get a data container from the GUIDE structure
data_copy,!g.s[2], mydc ; this is a copy
; do a linear interpolation to 1.5 x the channel spacing
dcresample,dc,dc.frequency_interval*1.5
show,dc
; this is exactly equivalent to "dcinvert"
dcresample,dc,-dc.frequency_interval
; get a fresh copy, use channel 0 as the keychan and
; do a spline interpolation, make sure things are ok
data_copy,!g.s[2],mydc
dcresample,dc,dc.frequency_interval*1.5,0,/spline,ok=ok
if ok then show,dc

Uses:

DATA_VALID dcinvert

pro dcresize(dc, newsize, zero=zero, beginning=beginning)#

Change the size of an existing data container.

If new channels are added to the end of an existing data container, they are filled with the blanked value (not a number, NAN) unless the /zero flag is set (in which case they are replaced with zeros). If /beginning is set then the extra channels are added to the beginning of the data and reference_channel value is adjust accordingly.

This only works for spectrum data containers where the x-axis is linear in frequency and so the x-axis values are always well-determined when extending a spectrum data container.

Params:

dcin, out, required, type=data container: The data container to resize.
newsizein, required, type=integer: The new number of channels. This must be > 0.

Keywords:

zeroin, optional, type=boolean: If set, any new channels are filled with zero instead of NAN.

Examples:

; double the size of dc, new channels are filled with NAN
nels = data_valid(dc)
dcresize,dc,nels*2
; back to its original size
dcresize,dc,nels
; add 100 channels, fill with 0.0
dcresize,dc,(nels+100),/zero
; to use this, or any toolbox function, on a DC in the !g
; structure, do this
dc = !g.s[0]
nels = data_valid(dc)
dcresizse,dc,nels+200

pro dcscale(dc, factor)#

This procedure scales a data container by a scalar value

Params:

dcin, required, type=data container: data container (spectrum or continuum)
factorin, required, type=float: scale factor

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show
dcscale,a,100
show,a

pro dcsetfframe(dc, toframe=toframe)#

Procedure to reset the frequency reference frame of the supplied data container(s) to the indicated frame.

The frequency axis reference frame is given by the frequency_type header value. Normally this remains “TOPO”, the topocentric frame in which the data were originally taken. The default GBTIDL display routines use the reference frame found in the velocity_definition header value. This procedure can be used to set the frequency axis values to that frame (or any other valid reference frame). This can be useful when exporting the data to an environment that doesn’t know about anything except the frequency information.

This routine will not change the displayed frequency axis.

If the frame is not set, the one implied by the data header velocity_definition is used.

This routine adjusts the following header values: reference_frequency, frequency_interval, frequency_type and center_frequency. When written to an SDFITS table, the first three of these four values correspond to CRVAL1, CDELT1, and CTYPE1, in that order.

Nothing is changed if the frame is already set to the desired frame.

Nothing is changed if the frame is unrecognized.

Params:

dcin, required, type=spectrum: The data container(s) to be adjusted.

Keywords:

toframein, optional, type=string: The reference frame to use. If not supplied, the value implied by the last 4 characters of the velocity_definition the data container. See frame_velocity for a full list of supported reference frames.

pro dcsetunits(dc, units, tau=tau, ap_eff=ap_eff, ret_tau=ret_tau, ret_ap_eff=ret_ap_eff, ok=ok)#

Set the units and convert the data to those units using the model described here.

This assumes that the units in the data container are in Ta, which is what they are immediately after the initial standard calibration (e.g. dofreqswitch). There is no attempt to look at the units field in the data container to see if this is appropriate or to undo any previous unit conversion prior to applying this conversion.

The recognized units and their associated scale factors are:

Ta - data are unchanged
Ta* - scale = exp(tau/sin(elevation))/0.99
Jy - scale = exp(tau/sin(elevation))/0.99 / (2.85*ap_eff)

The elevation is taken directly from the data container. Tau, the zenith opacity, may be supplied by the user. If not supplied, then the get_tau function is used, using the observed_frequency value from the data container. The aperture efficiency (ap_eff) can also be supplied by the user. If not supplied then the get_ap_eff function is used, using the observed_frequency from the data container.

Users are strongly encouraged to supply values for these keywords since the defaults are not very accurate.

Params:

dcin, required, type=data container: The data container to use.
unitsin, optional, type=string, default=’Ta’: The units to set, chosen from ‘Ta’,’Ta*’,and ‘Jy’.

Keywords:

tauin, optional, type=float: tau at zenith
ap_effin, optional, type=float: aperture efficiency
ret_tauout, optional, type=float: The tau actually used here.
ret_ap_effout, optional, type=float: The ap_eff actually used here.
okout, optional, type=boolean: This is 0 (false) if the units were unrecognized or tau and ap_eff were not between 0 and 1. In that case, the data are unchanged. Otherwise this is 1 (true).

pro dcshift(dc, offset, wrap=wrap, ftol=ftol, linear=linear, quadratic=quadratic, lsquadratic=lsquadratic, spline=spline, cubic=cubic, nowelsh=nowelsh, nopad=nopad, blanks=blanks, ok=ok)#

Procedure to shift data in a data container by a given number of channels.

Shift the data in the data container by the given number of channels (which may be a floating point number containing fractional channels) such that data in channel i before the shift is found in channel i+offset after the shift. The data’s reference channel is also shifted so that features should remain stationary when displayed except when the x-axis is channels. Shifted data is replaced with zeros unless the /wrap keyword is set in which case the data wraps around as necessary during the shift. Data are first shifted by the nearest integer number of channels and then the result is shifted by the remaining fractional channel if that fraction channel is more than ftol.

Currently, the default fractional shifting is done using an FFT windowed using a Welch function to reduce ringing as a result of the FFT. Set ftol to a number >= 1 to turn off all fractional shifting. See below for alternative methods. FFT is the default because that is the method that UniPOPS used in its SHIFT verb.

Use one of xshift, vshift, or fshift to calculate the offset to align this data with data in an ongoing accumulation.

Windowing using the Welsh function reduces ringing at the expense of the high-order terms in the FFT. This very slightly broadens features. This can be turned off using the /nowelsh keyword.

The data are padded with 0s to the next higher power of two as an intermediate step (the final result will have the original size of the data). This is done to reduce ringing due to any discontinuities between the two ends of the data. This can be turned off using the /nopad keyword. If there are any bad data points, they are interpolated using a linear interpolation prior to the FFT and then reblanked after the FFT.

Linear, quadratic, least squares quadratic, spline, and cubic interpolation are all available as alternatives to the default FFT based interpolation for the fractional part of the shift. The first 4 of these methods use the INTERPOL function and the cubic method uses the INTERPOLATE functions. Please see the IDL help on those functions for more information about them. The cubic function is a close approximation to a sinc function.

In some limited testing with GBT data, all of these interpolation methods agree reasonably well with each other. The linear interpolation is, not surprisingly, the fastest since it only uses 2 data points for each interpolated point. The quadratic and cubic methods both use 3 points for each interpolated point and both of those methods are not quite 2 times slower than the linear method. The FFT method (the default) is about another factor of 2 slower. Spline and least squares quadratic both use 4 data points for each interpolated point and they are substantially slower. Spline takes about 20 times longer and lsquadratic takes about 100 times as long as the linear interpolation.

If none of /linear, /spline, /quadratic, /lsquadratic, or /cubic are specified then an FFT is used for the fractional shift. It is an error to use more than one of these flags in the same call.

Params:

dcin, out, required, type=spectrum: The data container to shift. The shift is done in place.
offsetin, required, type=floating point: The number of channels to shift the data (positive shifts things towards higher channels, negative shifts things towards lower channels).

Keywords:

wrapin, optional, type=boolean: Data shifted off one end of the array appears on the other end of the array (it wraps around as a result of the shift) when this is set. Otherwise, as data is shifted it is blanked (replaced by NaNs) and data shifted off the end is lost.
ftolin, optional, type=floating point, default=0.01: Fractional shifts (the non-integer portion of offset) are only done when they are larger than ftol. Set this value to >= 1.0 to turn off all fractional shifts.
linearin, optional, type=boolean: When set, use the linear interpolation provided by INTERPOL for any fractional shift larger than ftol.
quadraticin, optional, type=boolean: When set, use the quadratic interpolation provided by INTERPOL for any fractional shift larger than ftol.
lsquadraticin, optional, type=boolean: When set, use the lsquadratic (lest squares quadratic) interpolation provided by INTERPOL for any fractional shift larger than ftol.
splinein, optional, type=boolean: When set, use the spline interpolation provided by INTERPOL for any fractional shift larger than ftol.
cubicin, optional, type=boolean: When set, use the cubic interpolation provided by INTERPOLATE for any fractional shift larger than ftol. The value of the CUBIC keyword in the INTERPOLATE call is set to -0.5.
nowelshin, optional, type=boolean: When set, the shifted data is NOT windowed using the Welsh function. This is ignored when a non-FFT-based fraction shift is done
nopadin, optional, type=boolean: When set, the data is NOT padded with 0s to the next higher power of 2 prior to the FFT and shift. The data are never padded for the non-FFT-based fractional shifts.
blanksout, optional, type=integer array: The shifted channel locations of all non-finite (blanks) channel values in the original, unshift data container. If there were no blank channels, this value is -1. If a blank channel is shifted out of the spectrum it is not included in this list (that may also result in a return value of -1). If there was a fractional shift and the resulting blank channel is now an extra channel wide then that is also reflected in this keyword. This array does not include the new blank values added during the shift to replace the shifted values. If the returned value is not -1 then all of the values in blanks are valid channel numbers and are blanked (not a number) in the shifted spectrum. This is primarily used by the dcfold procedure.
okout, optional, type=boolean: This is set to 1 on success or 0 on failure (e.g. bad arguments).

Uses:

DATA_VALID

pro dcsmooth(dc, newres, decimate=decimate, ok=ok)#

Smooth the data with a GAUSSIAN such that the spectrum in the given data container with an original resolution of frequency_resolution now has a resolution of newres, where newres is expressed in channels.

Optionally also decimate the spectrum by keeping only every newres channels.

The frequency_resolution field is set to newres * abs(frequency_interval) after this procedure is used.

The width of the smoothing Gaussian is sqrt(newres2-oldres^2).

Params:

dcin, required, type=data container: The data container to smooth.
newresin, required, type=real: The desired new resolution in units of channels. This must be >= the frequency_resoltuion also expressed in channels. If it is equal to the oldres value this procedure does not change the data.

Keywords:

decimatein, optional, type=boolean: When set, only every NEWRES channels are kept, starting from the original 0 channel. If NEWRES is not an integer, this may not be a wise thing to do (the decimation rounds to the nearest integer).
okout, optional, type=boolean: Returns 1 if everything went ok, 0 if it did not (missing parameter, empty or invalid dc, bad kernel).

Examples:

; dc already exists and is a valid data container
; smooth to 2 channels
dcsmooth,dc,2
; do it again, now to 4 channels, decimate this time
dcsmooth,dc,4,/decimate

Uses:

DATA_VALID dcconvol dcdecimate make_gauss_data

function dcspurchans(vsprval, vsprpix, vspdelt, nchan, docenterspur=docenterspur, count=count)#

Return a vector of VEGAS ADC spur channels given the values of VSPDELT, VSPRPIX, VSPRVAL, and NCHAN.

The VSP* values should be taken from a GBTIDL data container, where the VSPRPIX value has already had 1 subtracted from the value found in the FITS file (FITS counts channels from 1, GBTIDL counts channels from 0).

Spurs are located at the following channels: \(spur_j = (j - VSPRVAL)*VSPDELT + VSPRPIX\), where j is an integer that goes from 0 through 33. The VSP* values are double precision floats but the spur_j is assumed to be an integer. Conversion to an integer is done using the IDL round function.

NCHAN is used to determine which spur_j values should be returned. Valid spur_j values are equal to or greater than 0 and less than NCHAN. The count keyword is set to the number of valid spur_j values returned. If there are none, then the returned spur_j value is -1 and count is 0.

The value at j=16, the center ADC spur, is not included in the returned values unless the docenterspur keyword is set. Normal sdfits use will replace that center ADC spur with the average of the two adjacent channels and so it’s usually not necessary in GBTIDL to know where that center ADC spur is since it does not need to be reflagged or interpolated across.

These spurs do not include the center spur. That spur always occurs at the center channel (NCHAN/2 when counting from 0, as is done in IDL) and is replaced with the average of the two adjacent channels by sdfits unless the “-nointerp” option to sdfits is used. That spur does not arise in the ADC and so does not move as the spectral window is tuned across the ADC bandpass.

This routine does not check the validity of the input values.

Params:

vsprvalin, required, type=double: The spur reference value.
vsprpixin, required, type=double: The channel number (pixel) corresponding to the spur reference value (vsprval).
vspdeltin, required, type=double: The channel increment between adjacent spur numbers (j to j+1).
nchanin, required, type=integer: The number of channels to use in determining which spur_j values to return.

Keywords:

docenterspurin, optional, type=boolean: When set, the center ADC spur will be included in the array of returned channel numbers. The center ADC spur is normally interpolated across by sdfits and is typically not needed when reflagging or interpolating across the spurs in GBTIDL. The default behavior is to not return that channel in the list of spur channel numbers.
countout, optional, type=integer: The number of valid VEGAS spur channels returned.

Returns:

an array listing the channel numbers associated with all of the VEGAS spurs associated with the input parameters.

pro dcspurinterp(dc)#

Replace data values in the given data containers at the known locations of VEGAS spurs with the average of the two adjacent channel values.

VEGAS produces spurs (spikes) at channels corresponding to integer multiples of the ADC sampler frequency divided by 64. The normal behavior of sdfits is to flag these channels when the data are filled. When flagged data is retrieved by GBTIDL, the flagged channels are replaced with NaNs and will not contribute to most operations. If the skipflag option is used to retrieve the data then flagging is not applied and the original data values will appear at all channels.

It is sometimes useful to replace the spur data values (either NaNs if the flags have been used or the original data values if the flags are ignored) with interpolated values. This routine replaces those values with the average of the two channels on either side of the spur (or the adjacent channel if the spur is at the end of the spectrum).

This spur locations are determined using the VSPDELT, VSPRPIX, and VSPRVAL header values. For data filled using older versions of sdfits, these values are not present in the SDFITS table and will appear as NaN in GBTIDL. In that case, no spur interpolation is possible and a warning message will be printed. The data should be refilled by the most recent version of sdfits to make use of this procedure.

Params:

dcin, out, required, type=spectrum data container(s): The data container(s) to alter. May be an array of data containers.

Uses:

dcspurchans

function dcsubtract(dc1, dc2)#

This procedure subtracts the second data container’s data array from the first, returning the results.

Params:

dc1in, required, type=data container: data container (spectrum or continuum)
dc2in, required, type=data container: data container (spectrum or continuum)

Returns:

float array: dc2 data subtracted from dc1 data

Examples:

get,index=1
a = data_new()
data_copy,!g.s[0],a
show,a
get, index=2
b = data_new()
data_copy,!g.s[0],b
show,b
diff = dcsubtract(a,b)
plot,a

Uses:

dcpaircheck

function dcvshift(accumbuf, dc, frame=frame, veldef=veldef, voffset=voffset)#

Function to calculate the shift, in channels, necessary to align in velocity a data container with the data container template in an ongoing accumulation.

If the frame is not set, the one implied by the data header is used. Use dcxshift to align using the current settings of the plotter’s x-axis.

Params:

accumbufin, required, type=accum_struct: The ongoing accumulation buffer.
dcin, required, type=spectrum: The data container that needs to be shifted to align with the data container template in accumbuf.

Keywords:

framein, optional, type=string: The reference frame to use. If not supplied, the value implied by the last 4 characters of the velocity_definition in the ongoing accumulation will be used. See frame_velocity for a full list of supported reference frames.
veldefin, optional, type=string: The velocity definition to use. If not supplied, the value implied by the first 4 characters of the velocity_definition in the ongoing accumulation will be used. See frame_velocity for a full list of supported reference frames.
voffsetin, optional, type=double: A velocity offset in km/s to apply before aligning. If not set, this defaults to a value of 0.0. A typical use would be to set this equal to the source velocity (at the velocity header value). Not that the units expected here are km/s, in agreement with those expected for setvoffset. The units for the velocity header field are m/s.

Returns:

shift, in channels, to be used as argument to dcshift. Returns 0.0 on error.

Examples:

a={accum_struct}
getps,30
dcaccum, a, !g.s[0]  ; start an accum, no alignment needed yet
getps,31
vs = dcvshift(a,!g.s[0]) ; what is the shift to align 31 with 30?
; get a copy of data at !g.s[0]
data_copy,!g.s[0], d
dcshift, d, vs  ; actually shift the data
dcaccum, a, d ; and add it in
getps, 32
data_copy,!g.s[0], d
dcshift, d, dcvshift(a, d)  ; all in one line, shift 32 to align with 30
dcaccum, a, d
accumave, a, d ; result is in d now

function dcxshift(accumbuf, dc)#

This function returns a value that, when given as the argument to shift, will shift the given data container such that it is aligned in the current displayed x-axis with the data in the ongoing accumulation buffer. If there is no ongoing accumulation then this function returns 0. The units of the returned value are channels. The primary data container must be shifted that many channels in order to align in the current x-axis with the data in the accumulation buffer.

Params:

accumbufin, required, type=accum_struct: The ongoing accumulation buffer.
dcin, required, type=spectrum: The data container that needs to be shifted to align with the data container template in accumbuf.

Returns:

shift, in channels, to be used as argument to dcshift. Returns 0.0 on failure.

Examples:

Accumulate several PS scans

a={accum_struct}
getps,30
dcaccum, a, !g.s[0]  ; start an accum, no alignment needed yet
getps,31
xs = dcxshift(a,!g.s[0]) ; what is the shift to align 31 with 30?
; get a copy of data at !g.s[0]
data_copy,!g.s[0], d
dcshift, d, xs  ; actually shift the data
dcaccum, a, d ; and add it in
getps, 32
data_copy,!g.s[0], d
dcshift, d, dcxshift(a, d)  ; all in one line, shift 32 to align with 30
dcaccum, a, d
accumave, a, d ; result is in d now

function DECODE_VELDEF(veldef, velocity_definition, reference_frame)#

Parse the SDFITS VELDEF value into its two components, the velocity definition and velocity reference frame.

This value must contain 8 characters where the first 4 characters describe the velocity definition and the last 4 characters describe the reference frame. If the first 4 characters are recognized, they are expanded to be used in other toolbox functions (e.g. velocity_axis). Recognized velocity definitions and their expanded form are: RADI (RADIO), OPTI (OPTICAL), and RELA (TRUE). Any leading dash is removed from the last 4 characters before that value is set. In addition, OBS is translated to TOPO for topocentric. The return value is 0 if there were problems (non-standard velocity definition, wrong number of characters, etc) and 1 if everything appears to be okay.

Params:

veldefin, required, type=string: The value to decode.
velocity_definitionout, required, type=string: This will be one of the standard types, RADIO, OPTICAL, or TRUE, or a copy of the first 4 characters (if the latter, the return value will be 0).
reference_frameout, required, type=string: This will be the last 4 characters minus any leading dash.

Returns:

1 on success and 0 on failure.

function decompress_ints(strints, unique=unique, llimit=llimit, ulimit=ulimit)#

Converts a string containing an integer list in a special syntax into an integer array. The syntax is comma separated values, with : used for ranges. Example: 3,5,8:12 -> [3,5,8,9,10,11,12]

Optional upper and lower limits to the integer values can be specified. Values less than the lower limit are set to that limit and values greater than the upper limit are set to that limit.

If the unique keyword is set, then the returning array is sorted and contains only unique values.

Params:

strintsin, required, type=string: string containing integers

Returns:

integer array

function dmtodayno(day, mon, year)#

Convert from day of month, month, and year to daynumber of year. It also works with arrays.

This code came from Phil Perillat at Arecibo.

Local changes:

modify this documentation for use by idldoc.

Params:

dayin, required, type=long integer: day of month
monin, required, type=long integer: month of year 1..12
yearin, required, type=long integer: 4 digit year

Returns:

daynum, int/long daynumber of year. First day of year is 1.

function doboxcar1d(array, width, edge_truncate=edge_truncate, nan=nan, missing=missing)#

Do a boxcar smooth on a 1D array.

This uses the IDL smooth for odd smoothing widths but even widths are handled here since smooth rounds even widths to the next higher integer. Only 1-d arrays are handled here.

Params:

arrayin, required, type=1d array: The array to be smoothed.
widthin, required, type=integer: The width of the boxcar.

Keywords:

edge_truncatein, optional, type=boolean: Same meaning as for smooth.
nanin, optional, type=boolean: Same meaning as for smooth.
missingin, optional, type=same type as array, default=NaN: Same meaning as for smooth.

Returns:

the smoothed result.

pro docal(result, sig_off, sig_on, tcal=tcal)#

This procedure gets the cal signal from a single cal-switched integration.

The result is the difference of the data in the two data containers:

(*result.data_ptr) = (*sig_on.data_ptr - *sig_off.data_ptr)/2.0

The tsys in the result is meanTsys as calculated by dcmeantsys. The integration and exposure times in the result are the sum of those two times from each data container. All other header parameters in the result are copies of their values in the sig_off spectrum. dcmeantsys uses the mean_tcal value found in the sig_off data container unless the user supplies a tcal value using the tcal keyword. The mean_tcal value in result will reflect the actual tcal value used (as resulted by dcmeantsys).

This simple routine is designed to be called from a more complicated routine like getcal. This does not check the arguments for consistency or type.

It is the responsibility of the caller to ensure that result is freed using data_free when it is no longer needed (i.e. at the end of all anticipated calls to this function before returning to the calling level). Failure to do that will result in memory leaks. It is not necessary to free these data containers between consecutive calls to this function at the same IDL level (e.g. inside the same procedure).

Params:

resultout, required, type=spectrum: The result as described above.
sig_offin, out, required, type=spectrum: An uncalibrated spectrum with no cal signal.
sig_onin, required, type=spectrum: An uncalibrated spectrum with a cal signal.

Keyword:

tcalin, optional, type=float: A scalar value for the cal temperature (K). If not supplied. sig_off.mean_tcal will be used.

pro dofreqswitch(sigwcal, sig, refwcal, ref, smoothref, tsys=tsys, tau=tau, tcal=tcal, sigResult=sigResult, refResult=refResult)#

This procedure calibrates a single integration from a frequency switched scan.

The expected 4 spectra that are used here are the signal with no cal, the signal with cal, the reference with no cal and the reference with cal.

dototalpower is used to get the total power for the two signal spectra and for the two reference spectra.
These are then combined using dosigref to get the calibrated results.
- sigresult is done using the signal total power as the “signal” and the reference total power as the “reference”
- refresult is done using the reference reference total power as the “signal” and the signal total power as the “reference”.
See dototalpower and dosigref for details about each step in the process of producing the two results.

dcfold can then be used to combine sigresult and refresult to produce a folded result. That step does not happen here.

The user can optionally over-ride the reference system temperature calculated in dototalpower and used in dosigref by supplying a value for the tsys and tau keywords here. tsys is the system temperate at tau=0. If the user supplies this keyword, tsys is first adjusted to the elevation of the reference spectrum : \(tsys_used = tsys*exp(tau/sin(elevation)\).

If tau is not supplied, then the get_tau function is used, using the reference observed_frequency to arrive at a very rough guess as to the zenith opacity, tau. Users are encouraged to supply tau when they also supply tsys to improve the accuracy of this calculation. The adjusted tsys then becomes the reference spectrum’s tsys value for use in dosigref.

The units of sigresult and refresult are “Ta”. Use dcsetunits to change these units to something else.

This is used primarily by getfs and this code does almost no argument checks or sanity checks. The calling routine is expected to check that the 4 input spectra are compatible (all are valid data containers and all have the same number of data points).

It is the responsibility of the caller to ensure that sigResult and refResult are freed using data_free when their use is finished (i.e. at the end of all anticipated calls to this function before returning to the calling level). Failure to do that will result in memory leaks. It is not necessary to free these data containers between consecutive calls to this function at the same IDL level (e.g. inside the same procedure).

Params:

sigwcalin, required, type=spectrum: An uncalibrated spectrum from the signal phase with the cal on.
sigin, required, type=spectrum: An uncalibrated spectrum from the signal phase with the cal off.
refwcalin, required, type=spectrum: An uncalibrated spectrum from the reference phase with the cal on.
smoothrefin, optional, type=integer: Boxcar smooth width for reference spectrum. No smoothing if not supplied or if value is less than or equal to 1.

Keywords:

tsysin, optional, type=float: tsys at zenith, this is converted to a tsys at the observed elevation. If not suppled, the tsys for each integration is calculated as described elsewhere.
tauin, optional, type=float: tau at zenith, if not supplied, it is estimated using get_tau tau is only used when the requested units are other than the default of Ta and when a user-supplied tsys value at zenith is to be used.
sigResultout, required, type=spectrum: The result when using the signal phases as “sig” in dosigref.
refResultout, optional, type=spectrum: This result when using the reference phases as “sig” in dosigref.

pro dofullsigref(result, sigwcal, sig, refwcal, ref, smoothref, avgref, signocal=signocal, tsys=tsys, tau=tau, tcal=tcal, retsigtsys=retsigtsys, retreftsys=retreftsys)#

This procedure calibrates a single integration using four different uncalibrated spectrum data containers.

One pair of data containers are from the signal (source) and the other pair are the reference set of data containers. Within each pair, one spectrum has signal from the cal and the other does not.

dototalpower is used to get the total power in each pair (signal and reference).
These are then combined in dosigref to get the calibrated result.

See dototalpower and dosigref for details about each step.

The user can optionally over-ride the reference system temperature calculated in dototalpower and used in dosigref by supplying a value for the tsys and tau keywords here. tsys is the system temperate at tau=0. If the user supplies this keyword, tsys is first adjusted to the elevation of the reference spectrum : \(tsys_used = tsys*exp(tau/sin(elevation)\).

If tau is not supplied, then the get_tau function is used, using the reference observed_frequency to arrive at a very rough guess as to the zenith opacity, tau. Users are encouraged to supply tau when they also supply tsys to improve the accuracy of this calculation. The adjusted tsys then becomes the reference spectrum’s tsys value for use in dosigref.

The units of result is “Ta”. Use dcsetunits to change these units to something else.

This is used primarily by getsigref and this code does almost no argument checks or sanity checks. The calling routine is expected to check that the 4 input spectra are compatible (all are valid data containers and all have the same number of data points).

It is the responsibility of the caller to ensure that result is freed using data_free when it is no longer needed (i.e. at the end of all anticipated calls to this function before returning to the calling level). Failure to do that will result in memory leaks. It is not necessary to free these data containers between consecutive calls to this function at the same IDL level (e.g. inside the same procedure).

Params:

resultout, required, type=spectrum: The resulting spectrum data container.
sigwcalin, required, type=spectrum: An uncalibrated spectrum from the signal scan with the cal on.
sigin, required, type=spectrum: An uncalibrated spectrum from signal scan with the cal off.
refwcalin, required, type=spectrum: An uncalibrated spectrum from reference scan with the cal on.
refin, required, type=spectrum: An uncalibrated spectrum from reference scan with the cal off.
smoothrefin, optional, type=integer: Boxcar smooth width for reference spectrum. No smoothing if not supplied or if value is less than or equal to 1.
avgrefin, optional, type=boolean: If set, then the refwcal; is assumed to already be the average ref scan and no total power processing of the ref data is done here. The ref argument is unused in this case.

Keywords:

signocalin, optional, type=boolean: If set, then there is no CAL data for the sig state. sigwcal argument is ignored, only the sig argument is used as is.
tsysin, optional, type=float: tsys at zenith, this is converted to a tsys at the observed elevation. If not suppled, the tsys for each integration is calculated as described elsewhere.
tauin, optional, type=float: tau at zenith, if not supplied, it is estimated using get_tau tau is only used when the requested units are other than the default of Ta and when a user-supplied tsys value at zenith is to be used.
tcalin, optional, type=float: Cal temperature (K) to use in the Tsys calculation. If not supplied, the mean_tcal value from the header of the cal_off switching phase data in each integration is used. This must be a scalar, vector tcal is not yet supported.
retsigtsysout, optional, type=float: The reference Tsys used here.
retsigtsysout, optional, type=float: The signal Tsys calculated here. If tsys is supplied, then retsigtsys is equal to retreftsys. The signal Tsys is not actually used, mearly calculated and reported back through this keyword.

pro donod(result, s1_b1_off, s1_b1_on, s2_b1_off, s2_b1_on, s1_b2_off, s1_b2_on, s2_b2_off, s2_b2_on, smoothref, ret_tsys, avgref=avgref, tsys=tsys, tau=tau, tcal=tcal, eqweight=eqweight, single_beam=single_beam)#

This procedure calibrates a single integration from a Nod scan pair.

An integration in each scan in a standard Nod scan pair contains data from two beams (b1 and b2 here). And the data from each beam contains one spectrum with the cal on and one with no cal. There are then 8 spectra to be combined here. This assumes that data identified with scan 1 (s1) here has the signal data in beam 1 and the reference data in beam 2. Then in scan 2 (s2) the reference data is beam 1 and the signal data is in beam 2.

Each cal on, cal off pair for a specific beam and scan are combined using dototalpower to produce 4 spectra. Each spectra also has a tsys associated with it.
The signal from scan 2 (the beam 2 data) is combined with the reference from scan 1 (the beam 1 data) using dosigref
The signal from scan 1(the beam 1 data) is combined with the reference from scan 2 (the beam 2 data) also using dosigref.
These two calibrated spectra are combined using dcaccum to produce the final calibrated result.

See dototalpower and dosigref for details about each step.

The returned tsys (ret_tsys) are the 4 system temperatures calculed in the first step (from dototalpower). They are, in order, the system temperatures from scan 1 and beam 1 (signal in scan 1), scan 1 and beam 2 (reference in scan 1), scan 2 and beam 1 (reference in scan 2) and scan 2 and beam 2 (signal in scan 2). The system temperature in the result is a weighted combination of the two reference system temperatures as described in dcaccum.

The user can optionally over-ride the system temperatures calculated in dototalpower and used in dosigref by supplying a value for the tsys and tau keywords here. tsys is the system temperate at tau=0. If the user supplies this keyword, tsys is first adjusted to the elevation of the reference spectrum \(tsys_used = tsys*exp(tau/sin(elevation)\).

If tau is not supplied, then the get_tau function is used, using the reference observed_frequency to arrive at a very rough guess as to the zenith opacity, tau. Users are encouraged to supply tau when they also supply tsys to improve the accuracy of this calculation. The adjusted tsys then becomes the reference spectrum’s tsys value for use in dosigref. In that case, the returned tsys values are the corrected, user-supplied tsys values using the elevation from the appropriate scan.

The units of result is “Ta”. Use dcsetunits to change these units to something else.

This is used primarily by getnod and this code does almost no argument checks or sanity checks. The calling routine is expected to check that the 8 input spectra are compatible (all are valid data containers and all have the same number of data points).

It is the responsibility of the caller to ensure that result is freed using DATA_FREE when it is no longer needed (i.e. at the end of all anticipated calls to this function before returning to the calling level). Failure to do that will result in memory leaks. It is not necessary to free these data containers between consecutive calls to this function at the same IDL level (e.g. inside the same procedure).

Params:

resultout, required, type=data container: This is the resulting data container.
s1_b1_off: in, required, type=spectrum: An uncalibrated spectrum from feed 0 (source/tracking feed) in scan 1 with the cal off.
s1_b1_on: in, required, type=spectrum: An uncalibrated spectrum from feed 0 (source/tracking feed) in scan 1 with the cal on.
s2_b1_off: in, required, type=spectrum: An uncalibrated spectrum from feed 0 (reference/non-tracking feed) in scan 2 with the cal off.
s2_b1_on: in, required, type=spectrum: An uncalibrated spectrum from feed 0 (reference/non-tracking feed) in scan 2 with the cal on.
s1_b2_offin, required, type=spectrum: An uncalibrated spectrum from feed 1 (reference/non-tracking feed) in scan 1 with the cal off.
s1_b2_onin, required, type=spectrum: An uncalibrated spectrum from feed 1 (reference/non-tracking feed) in scan 1 with the cal on.
s2_b2_offin, required, type=spectrum: An uncalibrated spectrum from feed 1 (source/tracking feed) in scan 2 with the cal off.
s2_b2_onin, required, type=spectrum: An uncalibrated spectrum from feed 1 (source/tracking feed) in scan 2 with the cal on.
smoothrefin, optional, type=integer: Boxcar smooth width for reference spectrum. No smoothing if not supplied or if value is less than or equal to 1.
ret_tsysout, optional, type=float: internal use

Keywords:

avgref: in, optional, type=boolean: If set, then the off-signal terms (reference) are assumed to be already averaged and no total power processing of that data is done here. These terms are s1_b2_off, s1_b2_on, s2_b1_off, and s2_b1_on. Of these 4 terms, only s1_b2_off and s2_b1_off are actually used when this keyword is set.
tsysin, optional, type=float: tsys at zenith, this is converted to a tsys at the observed elevation. If not suppled, the tsys for each integration is calculated as described elsewhere.
tauin, optional, type=float: tau at zenith, if not supplied, it is estimated using get_tau tau is only used when the requested units are other than the default of Ta and when a user-supplied tsys value at zenith is to be used.
tcalin, optional, type=float: Cal temperature (K) to use in the Tsys calculation. If not supplied, the mean_tcal value from the header of the cal_off switching phase data in each integration is used. This must be a scalar, vector tcal is not yet supported.
eqweightin, optional, type=boolean: When set, average the two beams with equal weight, else use default weighting.
single_beamout, optional, type=boolean: When set, one of the two beams was blanked. Useful for reporting back to the user.

pro dosigref(dcresult, dcsig, dcref, smoothref)#

Combine signal and reference data containers.

The result is \(((dcsig-dcref)/dcref) * dcref.tsys\), where \(dcsig\) and \(dcref\) are the data values in the signal and reference data containers, respectively. The tsys in the result, dcresult.tsys, is equal to dcref.tsys.

The exposure time of the result is \(t_{res} = t_{sig} * t_{ref} * smoothref / (t_{sig}+t_{ref} * smoothref)\), where \(t_{sig}\) and \(t_{ref}\) are the exposures of the signal and reference data container and where smoothref=1 in when smoothref has not been supplied. That is the effective integration time appropriate for use in the radiometer equation. Note that in the two limits of smoothref (1=no smoothing and large smoothing), if \(t_{sig} = t_{ref}\) then \(exposure = 1/2 * t_{sig}\) or \(exposure = t_{sig}\). With no smothing, the noise increases due to the subtraction of an equally noisy reference spectrum. With smoothing, the noise is as it would be for just the signal spectrum. Smoothing should be done with caution, though, because it can emphasize systematic glitches present in both signal and reference spectra that are not fully subtracted as a consequence of the smoothing. All of the other header values are copied from dcsig.

If smoothref is set, the dcref values are smoothed using the IDL smooth function. This does a simple boxcar smooth where the smoothref parameter is the width of the boxcar. smoothref is only used if it is larger than 1. In certain cases this can improve the signal to noise ratio, but it may degrade baseline shapes and artificially emphasize spectrometer glitches. Use with care. A value of smoothref=16 is often a good choice

The dcresult data container is created as necessary. If it already exists, the internal pointer will be reused. It is the responsibility of the calling procedure or functiosdn to free that pointer using data_free. Failure to do that will result in a memory leak.

Params:

dcresultout, required, type=data_container: The result.
dcsigin, required, type=data_container: The signal data container.
dcrefin, required, type=data_container: The reference data container.
smoothrefin, optional, type=integer: Boxcar smooth width for reference spectrum. No smoothing if not supplied or if value is less than or equal to 1.

pro dototalpower(result, sig_off, sig_on, tcal=tcal)#

This procedure calibrates a single integration from a total power scan.

The result is the average of the data in the two data containers:

\[(*result.data_{ptr}) = (*sig_{off}.data_{ptr} + *sig_{on}.data_{ptr})/2.0\]

The tsys in the result is meanTsys as calculated by dcmeantsys. The integration and exposure times in the result are the sum of those two times from each data container. All other header parameters in the result are copies of their values in the sig_off spectrum. dcmeantsys uses the mean_tcal value found in the sig_off data container unless the user supplies a tcal value using the tcal keyword. The mean_tcal value in result will reflect the actual tcal value used (as resuted by dcmeantsys).

This simple routine is designed to be called from a more complicated routine like gettp. This does not check the arguments for consistency or type.

It is the responsibility of the caller to ensure that result is freed using DATA_FREE when it is no longer needed (i.e. at the end of all anticipated calls to this function before returning to the calling level). Failure to do that will result in memory leaks. It is not necessary to free these data containers between consecutive calls to this function at the same IDL level (e.g. inside the same procedure).

Params:

resultout, required, type=spectrum: The result as described above.
sig_offin, out, required, type=spectrum: An uncalibrated spectrum with no cal signal.
sig_onin, required, type=spectrum: An uncalibrated spectrum with a cal signal.

Keywords:

tcalin, optional, type=float: A scalar value for the cal temperature (K). If not supplied. sig_off.mean_tcal will be used.

function eqtogal(ra, dec, equinox)#

Translate equatorial coordiantes to galactic coordinates.

Params:

rain, required, type=double: RA, in degrees.
decin, required, type=double: Dec, in degrees.
equinoxin, required, type=double: The coordinate system EQUINOX in years.

Returns:

double precision [glon, glat] in degrees.

Uses:

glactc precess

function estboxres(width, oldres)#

Estimate the new frequency_resolution after boxcar smoothing something with a known old frequency_resolution.

This estimates the new frequency_resolution by assuming that the frequeny_resolution is the FWHM of a gaussian response function and that the ratio of the widths is equal to the inverse ratio of the heights of the response functions before and after the smoothing. Since a boxcar function isn’t a Gaussian, this will be inaccurate to some degree. The initial response function may also not be Gaussian. Tests against a full convolution as well as tests involving noise from real GBT data indicate that this estimate is accurate to within a few percent.

The boxcar procedure uses this function to adjust the frequency_resolution of the data container.

Params:

widthin, required, type=integer: The width of the boxcar.
oldresin, required, type=float: The frequency_resolution in channels before boxcar smoothing.

Returns:

the frequency resolution in channels after boxcar smoothing to the given width.

Examples:

; what would the new resolution be of the PDC after boxcar
; smoothing with a width of 5 channels.
print,

esthanres(5,!g.s[0].frequency_resolution/abs(!g.s[0].frequency_interval))

function esthanres(oldres)#

Estimate the new frequency_resolution after hanning smoothing something with a known old frequency_resolution.

This estimates the new frequency_resolution by assuming that the frequeny_resolution is the FWHM of a gaussian response function and that the ratio of the widths is equal to the inverse ratio of the heights of the response functions before and after the smoothing. Since a hanning function isn’t a Gaussian, this will be inaccurate to some degree. The initial response function may also not be Gaussian. Tests against a full convolution as well as tests involving noise from real GBT data indicate that this estimate is accurate to within a few percent.

The hanning procedure uses this function to adjust the frequency_resolution of the data container.

Params:

oldresin, required, type=float: The frequency_resolution in channels before hanning smoothing.

Returns:

the frequency resolution in channels after hanning smoothing.

Examples:

; what would the new resolution be of the PDC after hanning
; smoothing
print,

esthanres(!g.s[0].frequency_resolution/abs(!g.s[0].frequency_interval))

function extract_edges(strlist)#

Given a string containing a comma separated list with optional range strings, turn that into two strings giving the start and end points that would be used to turn that string into integers.

Params:

strlistin, required, type=string: Comma separated list

Returns:

2-element string array

Examples:

extract_edges(["1,2,4:8"]) returns ["1,2,4","1,2,8"]

function file_exists(filename, fullname, dir=dir, size=size)#

Check if a file name 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.

NOTE: this routine will only find regular files, a directory name will return a non-existant file.

This code came from Phil Perillat at Arecibo. Local changes:

modify this documentation for use by idldoc.

Params:

filenamein, required, type=string: filename to search for
fullnameout, optional, type=string: full directory/filename where file was found.

Keywords:

dirin, optional, type=string: If supplied then search through these directories. May be a vector.

sizein, out, optional: If supplied then set this to the file size in bytes, if it exists.

Returns:

1 file found, 0 not found

Examples:

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)

function fitsdateparse(dateobs)#

Convert a FITS date string into its component parts.

Params:

dateobsin, required, type=string: A FITS date string following the pattern: YYYY-MM-DD[Thh:mm:ss[.sss…]] or DD/MM/YY.

Returns:

double precision vector containing [year,month,day,hour,minute,second]. If the time fields are not present, 0s are returned. If the string is invalid, -1 is returned.

function frame_velocity(data, toframe, fromframe, bootstrap=bootstrap, status=status)#

Toolbox front-end to chdoppler. Get the velocity of a given frame relative to another frame.

Recognized frames:

TOPO	topocentric. The observed (sky) frame.
GEO	geocentric.
HEL	heliocentric.
BAR	barycentric.
LSR	Local standard of rest (kinematic).
LSD	Dynamic LSR.
GAL	galactocentric.

Params:

datain, required, type=spectrum: Data container to get pointing direction, time, and telescope location from.
toframein, required, type=string: The desired frame.
fromframein, optional, type=string: The originating frame. Defaults to “TOPO”.

Keywords:

bootstrapin, optional, type=boolean: If set, the data.frame_velocity will be used to get to the frame given in data.velocity_definition. The other data header information will only be used to get to other frames. Bootstraping is desirable because the sdfits tool may not be getting the correct times and positions for each integration - early versions didn’t. However, the frame_velocity is correct and can be relied on.
statusout, optional, type=integer: This is 1 if there were no problems and is 0 if anything unexpected happened (invalid data container type, velocity definition, or unrecognized data frames).

Returns:

Velocity difference between fromframe and toframe along the line of sight implied by the data header.

Uses:

chdoppler decode_veldef

function freqtochan(data, freqs, frame=frame, true_frame=true_frame)#

Convert frequency (Hz) to channel number using the supplied data container. The frequency may be expressed in another reference frame from that given in the data.

Params:

datain, required, type=spectrum: The spectrum data container to use to get the necessary header information.
freqsin, required;: The frequencies to convert (Hz), may be an array of values.

Keywords:

framein, optional, type=string: The rest frame that the frequencies are in. Known rest frames are listed in frame_velocity. Defaults to the frame given in the data container.
true_frameout, optional, type=string: The actual rest frame used in converting the frequencies. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as the frame in data.frequency_type.

Returns:

channel number.

Uses:

data_valid frame_velocity

function freqtofreq(dc, freq, toframe, fromframe)#

Convert a frequency in one rest frame to an equivalent frequency in a new rest frame.

Params:

dcin, required, type=spectrum: The spectrum data container to use to get the necessary header information for the conversion between the two frames.
freqin, required, type=double: The frequencies, in Hz, to convert.
toframein, required, type=string: The desired rest frame. Known rest frames are listed in frame_velocity.
fromframein, required, type=string: The rest frame appropriate for freq. Known rest frames are listed in frame_velocity.

Returns:

the converted frequency

function freqtovel(freq, restfreq, veldef=veldef)#

Convert frequency to velocity (m/s) using the given rest frequency and velocity definition. The units (Hz, MHz, GHz, etc) of the frequencies to convert must match that of the rest frequency argument.

Params:

freqin, required: Frequency. Units must be the same as the units of restfreq.
restfreqin, required: Rest frequency. Units must be the same as those of freq.

Keywords:

veldefin, optional, type=string: The velocity definition which must be one of OPTICAL, RADIO, or TRUE. Defaults to RADIO.

Returns:

velocity in m/s

function FREQUENCY_AXIS(data, frame=frame, true_frame=true_frame)#

Construct a frequency axis from a data structure.

Params:

datain, required, type=spectrum: The data structure to use in constructing a frequency axis.

Keywords:

framein, optional, type=string: The reference frame to use. If not supplied, the data.frequency_type value will be used. See frame_velocity for a full list of supported reference frames.
true_frameout, optional, type=string: The actual rest frame used in constructing the frequency axis. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as the frame in data.frequency_type.

Returns:

A vector of frequencies (Hz). May also set the keywords. Returns -1 on a severe error. If the ref_frame is invalid this will use the values found in the data structure and the values of the true_frame keyword will be set.

Uses:

chantofreq

function galtoeq(glon, glat, equinox)#

Translate galactic coordinates to equatorial coordiantes.

Params:

glonin, required, type=double: Galactic longitude, in degrees.
glatin, required, type=double: Galactic latitude, in degrees.
equinoxin, required, type=double: The coordinate system EQUINOX in years.

Returns:

double precision [ra,dec] in degrees.

Uses:

glactc precess

function gauss_fits(xx, yy, nregions, regions, inits, ngauss, max_iters, coefficients, errors, parinfo=parinfo, _EXTRA=ex)#

This code came mostly from Tom Bania’s GBT_IDL work. Local modifications include:

Allows multiple gauss fits in multiple regions
Results not printed here, nor analyzed: fits are just returned
All additional keywords used mpcurvefit can be used here as well.
Array syntax changed to [] and compile_opt idl2 used.
Indententation used to improve readability.
Some argument checks added.

Extra parameters are passed to mpcurvefit.

Params:

xxin, required, type=array: The x-values to use in the fit.
yyin, required, type=array: The data to be fit at xx.
nregionsin, required, type=long: The number of regions in which to fit gaussians
regionsin, required, type=2-D array: 2-D array marking ends of each region.
initsin, required, type=2-D array: 2-D array of the form [[h,c,w],[h,c,w],[h,c,w],…], where h = height, c = center, w = full width half maximum. Each [h,c,w] corresponds to a gaussian to fit.
ngaussin, required, type=integer: The total number of gaussians to fit.
max_itersin, required, type=long: max interations
coefficientsout, required, type=2-D array: 2-D array of the form [[h,c,w],[h,c,w],[h,c,w],…], where h, c, w are the results for the fit of the height, center, and width. If one of these was specified as fixed, it will return identical to its value in the inits param.
errorsout, required, type=2-D array: 2-D array of the form [[h,c,w],[h,c,w],[h,c,w],…], where h, c, w are the 1-sigma errors for the results returned in the coefficients parameter

Keywords:

parinfoin, optional, type=array

Array of structures for placing different constraints on each parameter: Each parameter is associated with one element of the array, in numerical order. The structure can have the following entries (none are required):

.VALUE

the starting parameter value (but see the START_PARAMS parameter for more information).

.FIXED

a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation.

.LIMITED

a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together.

.LIMITS

a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together.

.PARNAME

a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available.

.STEP

the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value.

.RELSTEP

the relative step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen.

.MPSIDE

the sidedness of the finite difference when computing numerical derivatives. This field can take four values:

0 - one-sided derivative computed automatically

1 - one-sided derivative (f(x+h) - f(x) )/h

-1 - one-sided derivative (f(x) - f(x-h))/h

2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)

Where h is the STEP parameter described above. The “automatic” one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0.

.MPMAXSTEP

the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value in one iteration. A value of 0 indicates no maximum. Default: 0.

.TIED

a string expression which “ties” the parameter to other free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo(2).tied = ‘2 * P(1)’. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [NOTE: the PARNAME can’t be used in expressions.]

.MPPRINT

if set to 1, then the default ITERPROC will print the parameter value. If set to 0, the parameter value will not be printed. This tag can be used to selectively print only a few parameter values out of many. Default: 1 (all parameters printed)

_EXTRAin, optional, type=record: Extra keywords are passed to mpcurvefit.

Returns:

fits array: the fit for each region, back to back. Use nregions and regions parameters to unwrap this result.

Examples:

simple example: one gaussian

; make a simple gauss
x = lindgen(150)
h = 400000.
c = 75.
w = 15.
noise = 10000
data=h*exp(-4.0*alog(2)*(x-c)^2/w^2)+(randomn(seed,n_elements(x))*noise)
; make an initial guess to this guassian
h = 400000.
c = 75.
w = 15.
inits = [h,c,w]
nregions = 1
regions = [[20,120]]
ngauss = 1
max_iters = 500
yfit = gauss_fits(x,data,nregions,regions,inits,ngauss,max_iters,coefficients,errors,quiet=1)
; view the results
plot, data
gbtoplot, x[regions[0]:regions[1]], yfit, color=!red, /chan

complex examle: multiple gaussians in multiple regions

; create 5 gaussians in the same plot
a1 = [400000.,35.,15.]
a2 = [100000.,15,7.5]
a3 = [200000.,110,8.0]
a4 = [100000.,150,5.5]
a5 = [100000.,170,5.5]
a = [a1,a2,a3,a4,a5]
x = lindgen(200)
data = make_gauss(x,a,10000.)
plot, data

; specify 3 regions
nregions = 3
regions = [[5,75],[90,130],[135,190]]
inits = [[a1],[a2],[a3],[a4],[a5]]
ngauss = 5
max_iters = 500
p = replicate({value:0.D, fixed:0, limited:[0,0], ;                      limits:[0.D,0]}, 15) 
                   ; 15 = 5 gauss * 3 parameter per guass
p[*].value = a
; hold the first gaussians height fixed
p[0].fixed = 1

; find all the fits at once
yfit = gauss_fits(x,data,nregions,regions,inits,ngauss,max_iters,coefficients,errors,parinfo=p,quiet=1)

; unwrap the results and plot them
ystart = 0
for i=0,(nregions-1) do begin
    b = regions[0,i]
    e = regions[1,i]
    yend = ystart + (e-b)
    y = yfit[ystart:yend]
    ystart = yend + 1
    gbtoplot, x[b:e], y, color=!red, /chan
endfor

pro gauss_fx(x, a, f, pder)#

The Gaussian function required by curvefit, and used by gauss_fits. Used to fit the sum of a number of Gaussians. Taken from ASTROLIB.

Params:

xin, required, type=float array: The x-location to evaluate the gaussians at.
ain, required, type=float array: The gaussians are described by a.
fout, required, type=float array: The evaluated value.
pderout, required, type=float array: The derivatives with respect to each parameter.

function gauss_plot_fn(params, minchan, maxchan, chanperpix, count=count)#

Generate x, y pairs to be plotted (x in channels) given a structure describing the number of gaussians and the current limits and resolution of the plotter.

This is used by gshow to call oplotfn so that the function is used to regenerate the gaussian values to be plotted as the plotter changes (e.g. zooming) which might otherwise lead to poorly sampled gaussians on unzoomed plots being unnecessarily jagged when zoomed in.

Most users will never use this directly although it could be used as a template when generating functions to use with oplotfn.

Params:

paramsin, required, type=structure: The fields in this structure are the names of the last 3 parameters is make_gauss_data. See the documention there for more information (fields are a, noise and offset).
minchanin, required, type=double: The minimum channel number of the current plot window’s x-axis.
maxchanin, required, type=double: The maximum channel number of the current plot window’s x-axis.
chanperpixin, required, type=double: The channels per pixel of the current plot window’s x-axis. This value is always positive.

Keywords:

countout, required, type=integer: The total number of points in arr to plot. Should be equal to 0 (do not plot any points) or the number of elements along the second axis of arr.

Returns:

arr - 2D array, arr[0,*] is the x values and arr[1,*] the y.

function get_ap_eff(freq, elev)#

Function to return a default aperture efficiency, given an observing frequency. Used in the calibration routines.

The formula comes from the Ruze equation, using a surface accuracy of 390 microns. Currently the aperture efficiency returned is independent of elevation. For observations at low elevations it is especially important for the observer to bypass the default value returned by this routine.

Params:

freqin, required, type=float: Observing frequency in GHz
elevin, optional, type=float: Observing elevation in degrees; currently not used.

Returns:

the aperture efficiency

Examples:

apeff = get_ap_eff(18.5)
print,apeff

function get_chans(dc, nregion, regions)#

Return a vector of channel numbers appropriate for the given data container.

If nregion is not set or is 0, all channels are returned. If nregion is set and is not zero then regions must be provided. Ignores any regions or parts of region that extend beyond the end of the data container.

Params:

dcin, required, type=data container: The data container to use.
nregionin, optional, type=integer: The number of regions to use. If not set, then all channels are returned.
regionsin, optional, type=2D array: Limit the channel numbers returned to only those channels that fall in these regions. This assumes the regions follows the same rules as used in !g.nregions. That the regions don’t overlap and are properly sorted with the first channel of each region at regions[0,*] and the last channel of each region at regions[1,*]. The nregion argument determines how many regions are used here. If nregion is supplied and is > 0 then regions is required.

Examples:

; equivalent to dindgen(n_elements(*dc.data_ptr))
c = get_chans(dc)
; get the chans corresponding to recently set global regions
c = get_chans(dc, !g.nregions, !g.regions)

Returns:

Array of channel numbers. Returns -1 if there was a problem.

Uses:

data_valid

function get_tau(freq)#

Function to return a default zenith opacity, given an observing frequency. Used in the calibration routines.

These values are simply best guesses. The user should replace this routine as required for the program at hand.

Params:

freqin, required, type=float: Observing frequency in GHz

Returns:

the zenith opacity

Examples:

apeff = get_ap_eff(18.5)
tau = get_tau(18.5)
print, apeff, tau

function get_usage(filename, method=method, verbose=verbose)#

Parses comment lines like this one to create a sensible usage statement. Remember that more then one method can reside in a single .pro file, so the method keyword may be needed.

Params:

filenamein, required, type=string: full path name to the .pro file whose comments are to be parsed

Keywords:

methodin, optional, type=string: the name of the method whose comments are to be parsed.
verbosein, optional, type=boolean: if not set, just the one line usage statement is returned, if set, full comments, parameter descriptions and other methods used list are also returned.

Returns:

the usage as a vector of strings. See usage for an example where this function is used.

function getazel(dc, quiet=quiet)#

Get the Azimuth and Elevation from a data container’s longitude_axis and latitude_axis header words.

This returns the Azimuth and Elevation as a 2-element vector (for spectral line data) or array (with dimension of [2, number of integrations] for continuum data) appropriate for the given data container. The values in the longitude_axis and latitude_axis are used as is the coordinate_mode. If the coordinate_mode is AZEL then those values are returned as is.

If the data container’s coordinate_mode field is OTHER then the values are returned as is and a warning message is emited. That warning message can be suppressed if /quiet is set.

Params:

dcin, required, type=data container: The data container to get the coordinate values from.

Keywords:

quietin, optional, type=boolean: When set, warning messages are suppressed.

Returns:

vector (spectral line) or array (continuum) holding Azimuth and Elevation in degrees.

function GETDCDATA(dc, elements)#

Convenience function for retrieving data from a data container. Can retrieve the entire data array, just a range, or a single element of the data.

Params:

dcin, required, type=struct: data container (spectrum or continuum)
elementsin, optional, type=long array: elements of data to return, one integer, or two element array specifiying range

Returns:

either the entire data array of data container, or part of it

Examples:

; dc already exists
; all of it
thedata = getdcdata(dc)
; some of it
somedata = getdcdata(dc,[100:400])

Uses:

DATA_VALID

function getgal(dc, quiet=quiet)#

Get galactic coordinates from a data container.

This returns the galactic longitude and latitude as a 2-element vector (for spectral line data) or array (with dimension of [2, the number of integrations] for continuum data) appropriate for the given data container. The values in the longitude_axis and latitude_axis are used as is the coordinate_mode. If the coordinate_mode is GALACTIC then those values are returned as is.

If the data container’s coordinate_mode field is OTHER then the values are returned as is and a warning message is emited. That warning message can be suppressed if /quiet is set.

Params:

dcin, required, type=data container: The data container to get the coordinate values from.

Keywords:

quietin, optional, type=boolean: When set, warning messages are suppressed.

Returns:

vector (spectral line) or array (continuum) holding galactic longitude and latitude, in degrees.

function getha(dc)#

Return the hour angle (HA) in decimal hours using the given data container. Note that this returns the value in hours, not decimal degrees as getradec does.

For continuum data, a vector of hour angles (one at each integration) is returned.

Params:

dcin, required, type=data_container_struct: The data container to use.

Returns:

HA in hours.

function gethadec(dc, equinox, quiet=quiet)#

Get the HA and DEC at the requested equinox from a data container.

This returns the Hour Angle and Declination as a 2-element vector (for spectral line data) or array (with dimension of [2, the number of integrations] for continuum data) appropriate for the given data container given the requested equinox. If equinox is not supplied, the value in the header is used if present and set, otherwise 2000.0 is used (e.g. if the data container holds GALACTIC coordinates). The values in the longitude_axis and latitude_axis are used as is the coordinate_mode. If the coordinate_mode is HADEC and the equinox argument is not supplied or it matches the equinox in the data container, then those values are returned as is.

If the data container’s coordinate_mode field is OTHER then the values are returned as is and a warning message is emited. That warning message can be suppressed if /quiet is set.

Params:

dcin, required, type=data container: The data container to get the coordinate values from.
equinoxin, optional, type=double: The equinox to use (decimal years). Defaults to the value of dc.equinox.

Keywords:

quietin, optional, type=boolean: When set, warning messages are suppressed.

Returns:

vector (spectral line) or array (continuum) holding HA and DEC in degrees

function getradec(dc, equinox, quiet=quiet)#

Get the RA and DEC at the requested equinox from a data container.

This returns the Right Ascension and Declination as a 2-element vector (for spectral line data) or array (with dimension of [2,the number of integrations] for continuum data) appropriate for the given data container given the requested equinox. If equinox is not supplied, the value in the header is used if present and set, otherwise 2000.0 is used (e.g. if the data container holds GALACTIC coordinates). The values in the longitude_axis and latitude_axis are used as is the coordinate_mode. If the coordinate_mode is RADEC and the equinox argument is not supplied or it matches the equinox in the data container, then those values are returned as is.

If the data container’s coordinate_mode field is OTHER then the values are returned as is and a warning message is emited. That warning message can be suppressed if /quiet is set.

Params:

dcin, required, type=data container: The data container to get the coordinate values from.
equinoxin, optional, type=double: The equinox to use (decimal years). Defaults to the value of dc.equinox.

Keywords:

quietin, optional, type=boolean: When set, warning messages are suppressed.

Returns:

vector (spectral line) or array (continuum) holding RA and DEC in degrees.

function ha2ra(ha, lst)#

Convert an hour angle (HA) in decimal hours to a Right Ascension (RA) in degrees using the given LST in seconds.

The returned value (RA) is always a postive number between 0 and 360. This works for vectors (lst must either be a single scalar or a vector of the same length as the HA argument).

Params:

hain, required, type=floating point: The hour angle(s) to convert (in degrees).
lstin, required, type=floating point: The LST (in seconds) to use in the conversion. If this is a vector, it must have the same number of elements as ha.

Returns:

RA (right ascension) in degrees. The returned value is always between 0.0 and 360.0. Returns NaN on error (missing arguments or bad number of elements in lst).

function isleapyear(year)#

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.

This code came from Phil Perillat at Arecibo.

Local changes:

modify this documentation for use by idldoc.

Params:

yearin, required, type=long integer: 4 digit year

Returns:

istat: int 0 if not a leap year, 1 if a leap year.

function juldaytolmst(juldat, obslong=obslong)#

Convert from Julian day to local mean sidereal time. By default the longitude of the GBT is used.

If you need local apparent sidereal time, then add the equation of the equinox to these values (see nutation_m()).

This code came from Phil Perillat at Arecibo.

Local changes:

use the position of the GBT as the default position.
modify this documentation for use by idldoc.
modify to make use of the obslong keyword, which is now east longitude in degrees to match contents of gbtidl data container.

Params:

juldatin, required, type=double: Array of julian days to convert.

Keywords:

obslongin, optional, type=double, default=GBT: East longitude of observatory in degrees.

Returns:

lmst[n]: double local mean sidereal time in radians

Uses:

utcinfoinp utcToUt1

function leftjustify(in, width, flag)#

Function to return the input string, left justified at the desired width. Whitespace is first trimmed from the input string according to flag and then enough spaces are appended to the input string to pad it out to the desired width. If the input string is wider than the desired width, after trimming, it is truncated so that the returned string always is of the desired width.

Params:

inin, required, type=string: The input string to pad.
widthin, required, type=integer: The desired width. Must be a positive integer.
flagin, optional, type=integer, default=0: This is passed directly to strtrim to control on in is trimmed prior to padding. If flag is 0, trailing whitespace is trimmed, if flag is 1, leading whitepsace is trimed, and if flag is 2, both are trimmed.

Returns:

the padded string

function make_gauss_data(x, a, noise, offset)#

Returns a gaussian, given coefficients and an x axis. Noise and offset can be added. Used in guide layer and in tests.

Params:

xin, required, type=array: x-axis
an, required, type=array: coefficents for gaussian: [height, center, width]
noisein, required, type=double: noise to inject
offsetin, optional, type=double: offset to add to gaussian

Returns:

the gaussian evaluated at x with noise added

Examples:

; create two gaussians
a = [400000.,75.,15.]
a = [a,[200000,35,5.]]
x = lindgen(150)
y = make_gauss_data, x, a, 100.
plot, x, y

function makefitsdate(mjd, precision=precision)#

Make a FITS date string from a modified Julian Date.

Params:

mjdin, required, type=double: Modified Julian Date, in days.

Keywords:

precisionin, optional, type=integer: Number of digits after the decimal point in the seconds field of the output FITS date. If precision is 0 (the default) no decimal point is used.

Returns:

a FITS date of the form YYYY-MM-DDThh:mm:ss[.sss].

Uses:

paddedstring

function ortho_fit(xx, yy, nfit, cfit, rms)#

Function uses general orthogonal polynomial to do least squares fitting.

This code came from Tom Bania’s GBT_IDL work. Local modifications include:

Documentation modified for use by idldoc
Array syntax changed to [] and compile_opt idl2 used.
Indententation used to improve readability.
Unnecessary code removed.
Some argument checks added.

Params:

xxin, required: The x-values to use in the fit.
yyin, required: The data to be fit at xx.
nfitin, required, type=integer: The order of the polynomial to fit.
cfitout, required: On return, cfit has coefficients of the polynomial \(fit = sum(m) a_m x^m\) because of round-off using this form gives unreliable results above about order 15 even for double precision. See ortho_poly for a better discussion on the contents of cfit.
rmsout, required: The rms error for each polynomial up to order nfit.

Returns:

polyfit array: contains information for fitting to arbitrary points using the recursion relations

Examples:

fit a 3rd order polynomial to the data in !g.s[0]

yy = *(!g.s[0].data_ptr)
xx = dindgen(n_elements(yy))
f = ortho_fit(xx, yy, 3, cfit, rms)

function ortho_poly(x, polyfit)#

function returns the sum of the set of orthogonal polynomials described by the polyfit argument evaluated at x.;

This code came from Tom Bania’s GBT_IDL work. Local modifications include:

Documentation modified for use by idldoc
Array syntax changed to [] and compile_opt idl2 used.
Indententation used to improve readability.
Unnecessary code removed.
Some argument checks added.
Changed name of cfit argument to polyfit to avoid confusion with cfit argument of ortho_fit, which is NOT the same thing.

Params:

xin, required

The x-values to use in evaluating the polynomials.

polyfitin, required, type=2D array

The array describing the polynomials. Typically this will be the return value from ortho_fit. The dimensionality is [4,(nfit+1)] * polyfit(3,m) gives the weighting coefficient for the m-th orthogonal polynomial * polyfit(0:2,m) for m=2 are the recursion coefficients for

\[p^(m) = x*p^(m-1)*polyfit(2,m) + p^(m-1)*polyfit(1,m) + p^(m-2)*polyfit(0,m) \]

polyfit(0,0) is the value of the constant term

\[p^(1) = polyfit(1,1)*x + polyfit(0,1)\]

Returns:

Array giving the evaluated polyfit at x.

Examples:

Fit a 7-th order polynomial to the data in !g.s[0] and subtract the result of the fit from that data, replacing the data in !g.s[0].

yy = *(!g.s[0].data_ptr)
xx = dindgen(n_elements(yy))
poly = ortho_fit(xx, yy, 7, cfit, rms)
thefit = ortho_poly(xx,poly)
*(!g.s[0].data_ptr) = *(!g.s[0].data_ptr) - thefit

function paddedstring(number, precision=precision)#

Format a number as a string, ensuring that it is padded with a leading zero so that there are always 2 digits before the optional decimal point.

This could be generalized for use with numbers larger than 2 digits before the decimal point, but that is all I needed for now.

Params:

numberin, required, type=numeric: The number to convert.

Keywords:

precisionin, optional, type=integer, default=0: The number of characters after the decimal point. If precision is 0, no decimal point appears in the returned value.

Returns:

string representation of number at given precision.

function ra2ha(ra, lst)#

Convert a Right Ascension (RA) in degrees to an hour angle (HA) in degrees using the given LST in seconds.

The returned value (HA) is always between -180.0 and +180.0. This works for vectors (lst must either be a single scalar or a vector of the same length as the RA argument).

Params:

rain, required, type=floating point: The right ascensions(s) to convert (in degrees).
lstin, required, type=floating point: The LST (in seconds) to use in the conversion. If this is a vector, it must have the same number of elements as ra.

Returns:

HA (hour angle) in degrees. The returned value is always between -180. and 180.0. Returns NaN on error (missing arguments or bad number of elements in lst).

pro scalevals(thevalues, scaledvalues, prefix)#

Scale a value and determine the appropriate prefix for the scaled result.

This could be more general, but it is only intended to work out scaling for Giga, Mega, and kilo (G,M,k). If thevalues is a vector, only the first value is examined and all values are scaled by the same factor.

Params:

thevaluesin, required, type=float: The values to be scaled.
scaledvaluesout, required, type=float: The scaled values.
prefixout, required, type=string: One of “G”, “M” or “k” for scalings by 1e9, 1e6, or 1e3 respectively.

function select_data(io_object, count=count, _EXTRA=ex)#

Select data from the given io object and return the array of matching indices.

Data can be selected based on entries in the index file, such as source name, polarization type, IF number, etc. For a complete list of eligible parameters use the procedure listcols

See the discussion on “Select” in the GBTIDL User’s Guide here for a summary of selection syntax.

The selection criteria are ultimately passed to the io class’s search_index via the _EXTRA parameter.

Params:

io_objectin, required, type=sdfits io object: The io object on which the selection is performed

Keywords:

countout, optional, type=integer: The number of matches found.
_EXTRAin, optional, type=extra keywords: These are the selection parameters.

Returns:

an array of indicies. Returns a value of -1 if no match was found.

function seq(s_beg, s_end, s_inc)#

Generates a sequence of numbers in an array

Params:

s_begin, required, type=int: first value
s_endin, required, type=int: last value
s_incin, optional, type=int: increment

Returns:

long integer array

pro setdcdata(dc, value, elements)#

Convenience function for setting data array of a data container.

Params:

dcin, required, type=struct: data container (spectrum or continuum)
valuein, required, type=float: single float or array float
elementsin, optional, type=long: elements of data to set, one integer, or two element array specifiying range

function shiftfreq(freq, voffset, veldef=veldef)#

Shift a frequency using a velocity offset. Note that because of the nature of the non-true velocity definitions, this is not a simply reversable shift for anything by a TRUE velocity. Namely, f=shiftfreq(shiftfreq(f,voffset),-voffset) will not result in the same value of f as you started out with for anything except veldef='TRUE'. If you want to undo a frequency shift using the same voffset, use unshiftfreq.

Params:

freqin, required, type=double: The frequencies, in Hz, to shift.
voffsetin, required, type=double: The velocity offset (m/s) to add to the frequencies (resulting frequencies are then appropriate for a frame moving at voffset relative to the original frame of freq).

Keywords:

veldefin, optional, type=string, default=’RADIO’: The velocity definition from one of RADIO, OPTICAL or TRUE.

Returns:

the shifted frequencies in Hz.

function shiftvel(vel, voffset, veldef=veldef)#

Shift a velocity using a velocity offset. Note that because of the nature of the non-true velocity definitions, this is not a simply reversable shift for anything by a TRUE velocity. Namely, v=shiftvel(shiftvel(v,voffset),-voffset) will not result in the same value of v as you started out with for anything except veldef='TRUE'. If you want to undo a velocity shift using the same voffset, use unshiftvel.

Params:

velin, required, type=double: The velocities, in m/s, to shift.
voffsetin, required, type=double: The velocity offset (m/s) to add to the velocities (resulting velocities are then appropriate for a frame moving at voffset relative to the original frame of vel).

Keywords:

veldefin, optional, type=string, default=’RADIO’: The velocity definition from one of RADIO, OPTICAL or TRUE.

Returns:

the shifted velocities in m/s.

function unshiftfreq(freq, voffset, veldef=veldef)#

Undo the effects of shiftfreq using the same voffset (do not use its negative here). See the comments in shiftfreq. This is intended to only be used on a quantity that is the result from shiftfreq.

Params:

freqin, required, type=double: The frequencies, in Hz, to unshift.
voffsetin, required, type=double: The velocity offset (m/s) to remove. This should be the same value used in a previous call to shiftfreq.

Keywords:

veldefin, optional, type=string, default=’RADIO’: The velocity definition from one of RADIO, OPTICAL or TRUE.

Returns:

the unshifted frequencies in Hz.

function unshiftvel(vel, voffset, veldef=veldef)#

Undo the effects of shiftvel using the same voffset (do not use its negative here). See the comments in shiftvel. This is intended to only be used on a quantity that is the result from shiftvel.

Params:

velin, required, type=double: The velocities, in m/s, to unshift.
voffsetin, required, type=double: The velocity offset (m/s) to remove. This should be the same value used in a previous call to shiftvel.

Keywords:

veldefin, optional, type=string, default=’RADIO’: The velocity definition from one of RADIO, OPTICAL or TRUE.

Returns:

the unshifted velocities in m/s.

pro usage(proname, verbose=verbose, source=source)#

Print out the usage for a given procedure or function.

If the procedure or function is found in the gbtidl pro tree then get_usage is used to produce the output, otherwise the argument is passed to doc_library.

GBTIDL uses idldoc to create our reference manual. The format used there does not lead to very useful output from doc_library (which simply shows the comment block in the file). get_usage knows what idldoc tags to look for and it formats things accordingly. It can also generate just a single “Usage:” line to remind the user what the syntax of a particular procedure or function is.

Params:

pronamein, required, type=string: The name of the procedure or function.

Keywords:

verbosein, optional, type=boolean, default=0: When set, a verbose usage is printed.
sourcein, optional, type=boolean, default=0: When set, the entire contents of the relevent .pro file is sent to the terminal through “more”. This keyword trumps verbose and no usage information is printed.

Uses:

get_usage which_routine

function utcinfoinp(juliandate, utcInfo)#

Input parameters to convert utc to ut1

Read the utc to ut1 information from the utcToUt1.dat file. The UTC_INFO structure will return the information needed to go from utc to UT1. the routine utcToUt1 converts from utc to ut1 using the information read in here into the UTC_INFO structure. The conversion algorithm is:

\[utcToUt1= ( offset + ((julDay - julDayAtOffset))*rate\]

The offset, rate, data are input from the utcToUt1.dat file.

The user passes in the julian date and the utcToUt1.dat file will be searched for the greatest julian date that is less than or equal to the date passed in. If all of the values are after the requested juliandate, the earliest value in the file will be used and and error will be returned.

Note

@utc_info must be done before this function is called.

This code came from Phil Perillat at Arecibo. Local changes:

modify to find the local copy of utcToUt1.dat file.
modify this documentation for use by idldoc.

Note

The file is updated whenever a leap second occurs or whenever the drift rate changes (usually every 6 months or a year). If you have downloaded this file from ao, then you need to redownload the newer versions occasionally. Check the file aodefdir()/data/pnt/lastUpdateTmStamp for when your file was last updated.

Params:

juliandatein, required, type=double: Julian date
utcInfoout, required, type=utcInfo structure;: utcInfo structure

Returns:

status; 0=problem, 1=ok

Uses:

file_exists dmtodayno daynotodm

function utcToUt1(juldat, utcInfo)#

utctout1 - convert utc to ut1

Return the offset from utc to ut1 as a fraction of a day. The returned value (dut1Frac) is defined as ut1Frac=utcFrac + dut1Frac; The fraction of a day can be less then 0.

The utc to ut1 conversion info is passed in via the structure UTC_INFO.

This code came from Phil Perillat at Arecibo.

Local changes:

modify this documentation for use by idldoc.

Params:

juldatin, required, type=double: julian date, may be a vector
utcInfoin, required, type=utc_info structure: The utc to ut1 conversion information.

Returns:

ut1FracOffset[n]: double add this to utc based times to get ut1

function VELOCITY_AXIS(data, veldef=veldef, frame=frame, true_frame=true_frame, true_veldef=true_veldef)#

Construct a velocity axis from a data structure.

Params:

datain, required, type=spectrum: The data structure to use in constructing a velocity axis.

Keywords:

veldefin, optional, type=string: The velocity definition to use, chosen from ‘RADIO’, ‘OPTICAL’, and ‘TRUE’. If not supplied, the data.velocity_definition value will be used.
framein, optional, type=string: The reference frame to use. If not supplied, the data.velocity_definition value will be used. See frame_velocity for a full list of supported reference frames.
true_frameout, optional, type=string: The actual rest frame used in constructing the velocity axis. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as the frame in data.velocity_definition.
true_veldefout, optional, type=string: The actual velocity frame used in constructing the velocity axis. The only way this will not equal the frame argument is if that argument was invalid. In that case, this keyword will be the same as data.velocity_definition.

Returns:

A vector of velocities (m/s). May also set the keywords. Returns -1 on a severe error. If the veldef or ref_frame are invalid this will use the values found in the data structure and the values of the true_veldef and true_frame keywords will be set.

Uses:

chantovel

function veltochan(data, vels, frame=frame, veldef=veldef)#

Convert velocity (m/s) to channel number using the supplied data container. The velocity’s frame and definition can also be supplied.

Params:

datain, required, type=spectrum: The spectrum data container to use to get the necessary header information.
velsin, required: The velocities (m/s) to convert, may be an array of values.

Keywords:

framein, optional, type=string: The rest frame that the velocities are in. Known rest frames are listed in frame_velocity. Defaults to the frame given in data.velocity_definition.
veldefin, optional, type=string: The velocity definition in use from RADIO, OPTICAL, or TRUE. Defaults to the value found in data.velocity_definition.

Returns:

channel number.

Uses:

DATA_VALID DECODE_VELDEF chantofreq freqtovel veltofreq freqtochan

function veltofreq(vel, restfreq, veldef=veldef)#

Convert velocity(m/s) to frequency using the given rest frequency and velocity definition. The units of the returned result (Hz, MHz, GHz, etc) are the same as that of the rest frequency argument. The velocity must be in m/s.

Params:

velin, required: Velocity in m/s
restfreqin, required: Rest frequency. The units of restfreq the units of the returned result.
veldefin, optional, type=string: The velocity definition which must be one of OPTICAL, RADIO, or TRUE. Defaults to RADIO.

Returns:

frequency in same units as restfreq argument.

function veltovel(vel, toveldef, fromveldef)#

Convert a velocity in one velocity definition to an equivalent velocity in another definition.

Params:

velin, required, type=double: The velocity, in m/s, to convert.
toveldefin, required, type=string: The desired velocity definition. Must be one of ‘TRUE’, ‘OPTICAL’ and ‘RADIO’.
fromveldefin, required, type=string: The input velocity definition. Must be one of ‘TRUE’, ‘OPTICAL’ and ‘RADIO’

Returns:

the converted velocity in m/s.

pro which(name)#

Search for any file in the IDL !path that contains the user-supplied IDL routine (procedure or function) name. Also indicates compilation status of each routine (in IDL lingo, whether or not the routine is “resolved”.)

Restrictions: The IDL !path is searched for file names that are simply the module (in IDL documentation, “module” and “routine” are used interchangeably) name with a “.pro” suffix appended to them. A module stored inside a file whose name is different than the module name (followed by a “.pro”) will not be found UNLESS that module happens to be the currently-resolved module! E.g., if the module “pro test_proc” lives in a file named “dumb_name.pro”, then it will not be found:

IDL> which, 'test_proc'
Module TEST_PROC Not Compiled.
% WHICH: test_proc.pro not found on IDL !path.

unless it happens to be resolved:

IDL> .run dumb_name
% Compiled module: TEST_PROC.
IDL> which, 'test_proc'
Currently-Compiled Module TEST_PROC in File:
/home/robishaw/dumb_name.pro

However, this is terrible programming style and sooner or later, if you hide generically-named modules in inappropriately-named files, bad things will (deservedly) happen to you.

The routine further assumes that a file named “dumb_name.pro” actually contains a module named “dumb_name”! If it doesn’t, then you are a bad programmer and should seek professional counseling.

Finally, if the user has somehow compiled a module as a procedure and then compiled a module of the same name as a function, they will both be available to the user, therefore both are listed. This situation should probably be avoided.

Notes: First, all currently-compiled procedures and functions are searched. Then the remainder of the IDL !path is searched. The current directory is searched before the IDL !path, whether or not the current directory is in the IDL !path, because this is the behavior of .run, .compile, .rnew, DOC_LIBRARY, etc.

Modification History:

2003, May 30	Written by Tim Robishaw, Berkeley
2004, Feb 17	Fixed oddity where user tries to call a function as if it were a procedure, thus listing the module in both the Compiled Functions and Compiled Procedures list.
2005, Jun 14	Split code into which_routine function for use elsewhere in GBTIDL and which. Reformatted comments for use with idldoc.
2006, May 11	Incorporate T. Robishaw changes found in original version of which into this version. Fixed scenario where two modules with the same name are compiled, one as a procedure, the other as a function. Add current directory to the top of the path and make sure the path is a unique list of directories. Added capability to cope with symbolic links. Fixed up warning and documentation for the strange case of having both a procedure and function of the same name compiled.

Params:

namein, required, type=string: The procedure or function name to search for.

Examples:

You haven’t yet resolved (compiled) the routine (module) DEFROI. Let’s look for it anyway:

IDL> which, 'defroi
Module DEFROI Not Compiled.

Other Files Containing Module DEFROI in IDL !path:
/usr/local/rsi/idl/lib/defroi.pro

For some reason you have two modules with the same name. (This can occur in libraries of IDL routines such as the Goddard IDL Astronomy User’s Library; an updated version of a routine is stored in a special directory while the old version is stored in its original directory.) Let’s see which version of the module ADSTRING we are currently using:

IDL> which, 'adstring.pro'
Currently-Compiled Module ADSTRING in File:
/home/robishaw/idl/goddard/pro/v5.4+/adstring.pro

Other Files Containing Module ADSTRING in IDL !path:
/home/robishaw/idl/goddard/pro/astro/adstring.pro

Uses:

which_routine

function which_routine(name, unresolved=unresolved)#

Search for any file in the IDL !path that contains the user-supplied IDL routine (procedure or function) name.

If the returned string has non-zero length, the routine has been compiled (resolved in IDL lingo) otherwise, it has not. Any unresolved files that contain this routine are also returned in the unresolved keyword value. If both the returned string and the resolved keyword values have zero length, then the routine could not be found in the !path.

This is the code behind <a href=”which.html”>which</a>. It was originally all in which but was split so that other code could make use of this functionality (e.g. driving a browser to the appropriate reference manual location, finding and summarizing the help information for a routine).

Restrictions: The IDL !path is searched for file names that are simply the module (in IDL documentation, “module” and “routine” are used interchangeably) name with a “.pro” suffix appended to them. A module stored inside a file whose name is different than the module name (followed by a “.pro”) will not be found UNLESS that module happens to be the currently-resolved module! E.g., if the module “pro test_proc” lives in a file named “dumb_name.pro”, then it will not be found:

IDL> a=which_routine('test_proc',unresolved=unresolved)
IDL> print, strlen(a), strlen(unresolved)
          0           0

unless it happens to be resolved:

IDL> .run dumb_name
% Compiled module: TEST_PROC.
IDL> print,which_routine('test_proc')
/hvc/robishaw/dumb_name.pro

However, this is terrible programming style and sooner or later, if you hide generically-named modules in inappropriately-named files, bad things will (deservedly) happen to you.

The routine further assumes that a file named “dumb_name.pro” actually contains a module named “dumb_name”! If it doesn’t, then you are a bad programmer and should seek professional counseling.

Note

First, all currently-compiled procedures and functions are searched. Then the remainder of the IDL !path is searched.

Modification History:

2003, May 30	Written by Tim Robishaw, Berkeley
2004, Feb 17	Fixed oddity where user tries to call a function as if it were a procedure, thus listing the module in both the Compiled Functions and Compiled Procedures list.
2005, Jun 14	Split code into which_routine function for use elsewhere in GBTIDL and which. Reformatted comments for use with idldoc.

Params:

namein, required, type=string: The procedure or function name to search for.

Keywords:

unresolvedout, type=string: The paths to files that likely contain this name but are not the currently compiled version containing name.

Examples:

You haven’t yet resolved (compiled) the routine (module) DEFROI. Let’s look for it anyway:

IDL> a=which_routine('defroi',unresolved=unresolved)
IDL> print,strlen(a)
      0
IDL> print, unresolved
/usr/local/rsi/idl/lib/defroi.pro

For some reason you have two modules with the same name. (This can occur in libraries of IDL routines such as the Goddard IDL Astronomy User’s Library; an updated version of a routine is stored in a special directory while the old version is stored in its original directory.) Let’s see which version of the module ADSTRING we are currently using:

IDL> a=which_routine('adstring.pro',unresolved=unresolved)
IDL> print,a
/hvc/robishaw/idl/goddard/pro/v5.4+/adstring.pro
IDL> print,unresolved
/hvc/robishaw/idl/goddard/pro/astro/adstring.pro

Returns:

String containing the file name from which the resolved (compiled) version of name was found. This is a zero-length string if name has not been resolved (compiled) yet.

.VALUE	the starting parameter value (but see the START_PARAMS parameter for more information).
.FIXED	a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation.
.LIMITED	a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together.
.LIMITS	a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together.
.PARNAME	a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available.
.STEP	the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value.
.RELSTEP	the relative step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen.
.MPSIDE	the sidedness of the finite difference when computing numerical derivatives. This field can take four values: 0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h -1 - one-sided derivative (f(x) - f(x-h))/h 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) Where h is the STEP parameter described above. The “automatic” one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0.
.MPMAXSTEP	the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value in one iteration. A value of 0 indicates no maximum. Default: 0.
.TIED	a string expression which “ties” the parameter to other free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo(2).tied = ‘2 * P(1)’. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [NOTE: the PARNAME can’t be used in expressions.]
.MPPRINT	if set to 1, then the default ITERPROC will print the parameter value. If set to 0, the parameter value will not be printed. This tag can be used to selectively print only a few parameter values out of many. Default: 1 (all parameters printed)