; NAME:
; EQ2HOR
;
; PURPOSE:
; Convert celestial (ra-dec) coords to local horizon coords (alt-az).
;
; CALLING SEQUENCE:
;
; eq2hor, ra, dec, jd, alt, az, [ha, LAT= , LON= , /WS, OBSNAME= , $
; /B1950 , PRECESS_= 0, NUTATE_= 0, REFRACT_= 0, $
; ABERRATION_= 0, ALTITUDE= , /VERBOSE, _EXTRA= ]
;
; DESCRIPTION:
; This is a nice code to calculate horizon (alt,az) coordinates from equatorial
; (ra,dec) coords. It is typically accurate to about 1 arcsecond or better (I
; have checked the output against the publicly available XEPHEM software). It
; performs precession, nutation, aberration, and refraction corrections. The
; perhaps best thing about it is that it can take arrays as inputs, in all
; variables and keywords EXCEPT Lat, lon, and Altitude (the code assumes these
; aren't changing), and uses vector arithmetic in every calculation except
; when calculating the precession matrices.
;
; INPUT VARIABLES:
; RA : Right Ascension of object (J2000) in degrees (FK5); scalar or
; vector.
; Dec : Declination of object (J2000) in degrees (FK5), scalar or vector.
; JD : Julian Date [scalar or vector]
;
; Note: if RA and DEC are arrays, then alt and az will also be arrays.
; If RA and DEC are arrays, JD may be a scalar OR an array of the
; same dimensionality.
;
; OPTIONAL INPUT KEYWORDS:
; lat : north geodetic latitude of location in degrees
; lon : EAST longitude of location in degrees (Specify west longitude
; with a negative sign.)
; /WS : Set this to get the azimuth measured westward from south (not
; East of North).
; obsname: Set this to a valid observatory name to be used by the
; astrolib OBSERVATORY procedure, which will return the latitude
; and longitude to be used by this program.
; /B1950 : Set this if your ra and dec are specified in B1950, FK4
; coordinates (instead of J2000, FK5)
; precess_ : Set this to 1 to force precession [default], 0 for no
; precession correction
; nutate_ : Set this to 1 to force nutation [default], 0 for no nutation.
; aberration_ : Set this to 1 to force aberration correction [default],
; 0 for no correction.
; refract_ : Set to 1 to force refraction correction [default], 0 for no
; correction.
; altitude: The altitude of the observing location, in meters. [default=0].
; verbose: Set this for verbose output. The default is verbose=0.
; _extra: This is for setting TEMPERATURE or PRESSURE explicity, which are
; used by CO_REFRACT to calculate the refraction effect of the
; atmosphere. If you don't set these, the program will make an
; intelligent guess as to what they are (taking into account your
; altitude). See CO_REFRACT for more details.
;
; OUTPUT VARIABLES: (all double precision)
; alt : altitude (in degrees)
; az : azimuth angle (in degrees, measured EAST from NORTH, but see
; keyword WS above.)
; ha : hour angle (in degrees) (optional)
;
; DEPENDENCIES:
; NUTATE, PRECESS, OBSERVATORY, SUNPOS, ADSTRING() (from the astrolib)
; CO_NUTATE, CO_ABERRATION, CO_REFRACT, ALTAZ2HADEC
;
; BASIC STEPS
; Apply refraction correction to find apparent Alt.
; Calculate Local Mean Sidereal Time
; Calculate Local Apparent Sidereal Time
; Do Spherical Trig to find apparent hour angle, declination.
; Calculate Right Ascension from hour angle and local sidereal time.
; Nutation Correction to Ra-Dec
; Aberration correction to Ra-Dec
; Precess Ra-Dec to current equinox.
;
;
;CORRECTIONS I DO NOT MAKE:
; * Deflection of Light by the sun due to GR. (typically milliarcseconds,
; can be arseconds within one degree of the sun)
; * The Effect of Annual Parallax (typically < 1 arcsecond)
; * and more (see below)
;
; TO DO
; * Better Refraction Correction. Need to put in wavelength dependence,
; and integrate through the atmosphere.
; * Topocentric Parallax Correction (will take into account elevation of
; the observatory)
; * Proper Motion (but this will require crazy lookup tables or something).
; * Difference between UTC and UT1 in determining LAST -- is this
; important?
; * Effect of Annual Parallax (is this the same as topocentric Parallax?)
; * Polar Motion
; * Better connection to Julian Date Calculator.
;
; EXAMPLE
;
; Find the position of the open cluster NGC 2264 at the Effelsburg Radio
; Telescope in Germany, on June 11, 2023, at local time 22:00 (METDST).
; The inputs will then be:
;
; Julian Date = 2460107.250
; Latitude = 50d 31m 36s
; Longitude = 06h 51m 18s
; Altitude = 369 meters
; RA (J2000) = 06h 40m 58.2s
; Dec(J2000) = 09d 53m 44.0s
;
; IDL> eq2hor, ten(6,40,58.2)*15., ten(9,53,44), 2460107.250d, alt, az, $
; lat=ten(50,31,36), lon=ten(6,51,18), altitude=369.0, /verb, $
; pres=980.0, temp=283.0
;
; The program produces this output (because the VERBOSE keyword was set)
;
; Latitude = +50 31 36.0 Longitude = +06 51 18.0
; Julian Date = 2460107.250000
; Ra, Dec: 06 40 58.2 +09 53 44.0 (J2000)
; Ra, Dec: 06 42 15.7 +09 52 19.2 (J2023.4422)
; Ra, Dec: 06 42 13.8 +09 52 26.9 (fully corrected)
; LMST = +11 46 42.0
; LAST = +11 46 41.4
; Hour Angle = +05 04 27.6 (hh:mm:ss)
; Az, El = 17 42 25.6 +16 25 10.3 (Apparent Coords)
; Az, El = 17 42 25.6 +16 28 22.8 (Observer Coords)
;
; Compare this with the result from XEPHEM:
; Az, El = 17h 42m 25.6s +16d 28m 21s
;
; This 1.8 arcsecond discrepancy in elevation arises primarily from slight
; differences in the way I calculate the refraction correction from XEPHEM, and
; is pretty typical.
;
; AUTHOR:
; Chris O'Dell
; Univ. of Wisconsin-Madison
; Observational Cosmology Laboratory
; Email: odell@cmb.physics.wisc.edu
; NAME:
; FILTER_IMAGE
;
; PURPOSE:
; Identical to MEDIAN or SMOOTH but handle edges and allow iterations.
; EXPLANATION:
; Computes the average and/or median of pixels in moving box,
; replacing center pixel with the computed average and/or median,
; (using the IDL SMOOTH() or MEDIAN() functions).
; The main reason for using this function is the options to
; also process the pixels at edges and corners of image, and,
; to apply iterative smoothing simulating convolution with Gaussian,
; and/or to convolve image with a Gaussian kernel.
;
; CALLING SEQUENCE:
; Result = filter_image( image, SMOOTH=width, MEDIAN = width, /ALL_PIXELS
; /ITERATE, FWHM =, /NO_FT_CONVOL)
;
; INPUT:
; image = 2-D array (matrix)
;
; OPTIONAL INPUT KEYWORDS:
; SMOOTH = scalar (odd) integer specifying the width of a square box
; for moving average, in # pixels. /SMOOTH means use box
; width = 3 pixels for smoothing.
;
; MEDIAN = scalar (usually odd) integer specifying the width of square
; moving box for median filter, in # pixels. /MEDIAN means use
; box width = 3 pixels for median filter.
;
; /ALL_PIXELS causes the edges of image to be filtered as well. This
; is accomplished by reflecting pixels adjacent to edges outward
; (similar to the /EDGE_WRAP keyword in CONVOL).
; Note that this is a different algorithm from the /EDGE_TRUCATE
; keyword to SMOOTH or CONVOL, which duplicates the nearest pixel.
;
; /ITERATE means apply smooth(image,3) iteratively for a count of
; (box_width-1)/2 times (=radius), when box_width >= 5.
; This is equivalent to convolution with a Gaussian PSF
; of FWHM = 2 * sqrt( radius ) as radius gets large.
; Note that /ALL_PIXELS is automatically applied,
; giving better results in the iteration limit.
; (also, MEDIAN keyword is ignored when /ITER is specified).
;
; FWHM_GAUSSIAN = Full-width half-max of Gaussian to convolve with image.
; FWHM can be a single number (circular beam),
; or 2 numbers giving axes of elliptical beam.
;
; /NO_FT_CONVOL causes the convolution to be computed directly,
; with intrinsic IDL CONVOL function. The default is to use
; FFT when factors of size are all LE 13. Note that
; external function convolve.pro handles both cases)
;
; OPTIONAL INPUT/OUTPUT KEYWORD:
; PSF = Array containing the PSF used during the convolution. This
; keyword is only active if the FWHM_GAUSSIAN keyword is also
; specified. If PSF is undefined on input, then upon output it
; contains the Gaussian convolution specified by the FWHM_GAUSSIAN
; keyword. If the PSF array is defined on input then it is used
; as the convolution kernel, the value of the FWHM_GAUSSIAN keyword
; is ignored. Typically, on a first call set PSF to an undefined
; variable, which can be reused for subsequent calls to prevent
; recalculation of the Gaussian PSF.
; RESULT:
; Function returns the smoothed, median filtered, or convolved image.
; If both SMOOTH and MEDIAN are specified, median filter is applied first.
;
; EXAMPLES:
; To apply 3x3 moving median filter and
; then 3x3 moving average, both applied to all pixels:
;
; Result = filter_image( image, /SMOOTH, /MEDIAN, /ALL )
;
; To iteratively apply 3x3 moving average filter for 4 = (9-1)/2 times,
; thus approximating convolution with Gaussian of FWHM = 2*sqrt(4) = 4 :
;
; Result = filter_image( image, SMOOTH=9, /ITER )
;
; To convolve all pixels with Gaussian of FWHM = 3.7 x 5.2 pixels:
;
; Result = filter_image( image, FWHM=[3.7,5.2], /ALL )
;
; EXTERNAL CALLS:
; function psf_gaussian
; function convolve
; pro factor
; function prime ;all these called only if FWHM is specified
;
; PROCEDURE:
; If both /ALL_PIXELS (or /ITERATE) keywords are set then
; create a larger image by reflecting the edges outward, then call the
; IDL MEDIAN() or SMOOTH() function on the larger image, and just return
; the central part (the original size image).
;
; NAN values are recognized during calls to MEDIAN() or SMOOTH(), but
; not for convolution with a Gaussian (FWHM keyword supplied).
; HISTORY:
; Written, 1991, Frank Varosi, NASA/GSFC.
; FV, 1992, added /ITERATE option.
; FV, 1993, added FWHM_GAUSSIAN= option.
; Converted to IDL V5.0 W. Landsman September 1997
; Use /EVEN call to median, recognize NAN values in SMOOTH
; W. Landsman June 2001
; Added PSF keyword, Bjorn Heijligers/WL, September 2001
; NAME:
; FIND_WITH_DEF()
; PURPOSE:
; Searches for files with a default path and extension.
; EXPLANATION:
; Finds files using default paths and extensions, similar to using the
; DEFAULT keyword with the OPEN statement in VMS. Using this routine
; together with environment variables allows an OS-independent approach
; to finding files.
; CALLING SEQUENCE:
; Result = FIND_WITH_DEF( FILENAME, PATHS [, EXTENSIONS ] )
;
; INPUTS:
; FILENAME = Name of file to be searched for. It may either be a
; complete filename, or the path or extension could be left
; off, in which case the routine will attempt to find the
; file using the default paths and extensions.
;
; PATHS = One or more default paths to use in the search in case
; FILENAME does not contain a path itself. The individual
; paths are separated by commas, although in UNIX, colons
; can also be used. In other words, PATHS has the same
; format as !PATH, except that commas can be used as a
; separator regardless of operating system. The current
; directory is always searched first, unless the keyword
; NOCURRENT is set.
;
; A leading $ can be used in any path to signal that what
; follows is an environmental variable, but the $ is not
; necessary. (In VMS the $ can either be part of the path,
; or can signal logical names for compatibility with Unix.)
; Environmental variables can themselves contain multiple
; paths.
;
; OPTIONAL INPUTS:
; EXTENSIONS = One or more extensions to append to end of filename if the
; filename does not contain one (e.g. ".dat"). The period
; is optional. Multiple extensions can be separated by
; commas or colons.
; OUTPUTS:
; The result of the function is the name of the file if successful, or
; the null string if unsuccessful.
; OPTIONAL INPUT KEYWORDS:
; NOCURRENT = If set, then the current directory is not searched.
;
; RESET = The FIND_WITH_DEF routine supports paths which are
; preceeded with the plus sign to signal that all
; subdirectories should also be searched. Often this is
; used with logical names. It can be rather slow to search
; through these subdirectories. The /RESET keyword can be
; used to redefine an environment variable so that
; subsequent calls don't need to look for the
; subdirectories.
;
; To use /RESET, the PATHS parameter must contain the name
; of a *single* environment variable. For example
;
; setenv,'FITS_DATA=+/datadisk/fits'
; file = find_with_def('test.fits','FITS_DATA',/reset)
;
; EXAMPLE:
;
; FILENAME = ''
; READ, 'File to open: ', FILENAME
; FILE = FIND_WITH_DEF( FILENAME, 'SERTS_DATA', '.fix' )
; IF FILE NE '' THEN ...
;
;
; PROCEDURE CALLS:
; BREAK_PATH(), FIND_ALL_DIR()
; REVISION HISTORY:
; Version 1, William Thompson, GSFC, 3 May 1993.
; Removed trailing / and : characters.
; Fixed bugs
; Allow for commas within values of logical names.
; Added keyword NOCURRENT.
; Changed to call BREAK_PATH
; Version 2, William Thompson, GSFC, 3 November 1994
; Made EXTENSIONS optional.
; Version 3, William Thompson, GSFC, 30 April 1996
; Call FIND_ALL_DIR to resolve any plus signs.
; Version 4, S.V. Haugan, UiO, 5 June 1996
; Using OPENR,..,ERROR=ERROR to avoid an IDL 3.6
; internal nesting error.
; Version 5, R.A. Schwartz, GSFC, 11 July 1996
; Use SPEC_DIR to interpret PATH under VMS
; Version 6, William Thompson, GSFC, 5 August 1996
; Took out call to SPEC_DIR (i.e., reverted to version 4). The
; use of SPEC_DIR was required to support logical names defined
; via SETLOG,/CONFINE. However, it conflicted with the ability
; to use logical names with multiple values. Removing the
; /CONFINE made it unnecessary to call SPEC_DIR in this routine.
; Version 7, William Thompson, GSFC, 6 August 1996
; Added keyword RESET
; Converted to IDL V5.0 W. Landsman October 1997
; Use STRTRIM instead of TRIM, W. Landsman November 1998
; Use STRSPLIT instead of STR_SEP(), V5.3 or later W.L. July 2002
; NAME:
; FITSDIR
; PURPOSE:
; Display selected FITS keywords from the headers of FITS files.
; EXPLANATION:
; The values of either user-specified or default FITS keywords are
; displayed in either the primary header and/or the first extension header.
; Unless the /NOSIZE keyword is set, the data size is also displayed.
; The default keywords are as follows (with keywords in 2nd row used if
; those in the first row not found, and the 3rd row if neither the keywords
; in the first or second rows found:)
;
; DATE-OBS TELESCOP OBJECT EXPTIME
; TDATEOBS TELNAME TARGNAME INTEG ;First Alternative
; DATE OBSERVAT EXPOSURE ;Second Alternative
; INSTRUME EXPTIM ;Third Alternative
;
; FITSDIR will also recognize gzip compressed files (must have a .gz
; extension).
; CALLING SEQUENCE:
; FITSDIR , [ directory, TEXTOUT =, /FLAT, KEYWORDS=, /NOSIZE, /NoTELESCOPE
; ALT1_KEYWORDS= ,ALT2_KEYWORDS = ,ALT3_KEYWORDS =
;
; OPTIONAL INPUT PARAMETERS:
; DIRECTORY - Scalar string giving file name, disk or directory to be
; searched. Wildcard file names are allowed. Examples of
; valid names include 'iraf/*.fits' (Unix), d:\myfiles\f*.fits',
; (Windows) or 'Macintosh HD:Files:*c0f.fits' (Macintosh).
;
; OPTIONAL KEYWORD INPUT PARAMETER
; KEYWORDS - FITS keywords to display, as either a vector of strings or as
; a comma delimited scalar string, e.g.'testname,dewar,filter'
; If not supplied, then the default keywords are 'DATE-OBS',
; 'TELESCOP','OBJECT','EXPTIME'
; ALT1_KEYWORDS - A list (either a vector of strings or a comma delimited
; strings of alternative keywords to use if the default
; KEYWORDS cannot be found. By default, 'TDATEOBS', is the
; alternative to DATE-OBS, 'TELNAME' for 'TELESCOP','TARGNAME'
; for 'OBJECT', and 'INTEG' for EXPTIME
; ALT2_KEYWORDS - A list (either a vector of strings or a comma delimited
; strings of alternative keywords to use if neither KEYWORDS
; nor ALT1_KEYWORDS can be found.
; ALT3_KEYWORDS - A list (either a vector of strings or a comma delimited
; strings of alternative keywords to use if neither KEYWORDS
; nor ALT1_KEYWORDS nor ALT2_KEYWORDS can be found.
; /NOSIZE - if set then information about the image size is not displayed
; TEXTOUT - Controls output device as described in TEXTOPEN procedure
; textout=1 TERMINAL using /more option
; textout=2 TERMINAL without /more option
; textout=3 .prt
; textout=4 laser.tmp
; textout=5 user must open file
; textout=7 Append to existing .prt file
; textout = filename (default extension of .prt)
; /EXTEN - If set, then the first FITS extension is also checked for the
; desired keywords.
; /NOTELESCOPE - If set, then if the default keywords are used, then the
; TELESCOPE (or TELNAME, OBSERVAT, INSTRUME) keywords are omitted
; to give more room for display other keywords. The /NOTELESCOP
; keyword has no effect if the default keywords are not used.
; OUTPUT PARAMETERS:
; None.
;
; EXAMPLES:
; (1) Print info on all'*.fits' files in the current directory using default
; keywords. Include information from the extension header
; IDL> fitsdir,/exten
;
; (2) Write a driver program to display selected keywords in HST/ACS drizzled
; (*drz) images
; pro acsdir
; keywords = 'date-obs,targname,detector,filter1,filter2,exptime'
; fitsdir,'*drz.fits',key=keywords,/exten
; return & end
;
; (3) Write info on all *.fits files in the Unix directory /usr2/smith, to a
; file 'smith.txt' using the default keywords, but do not display the
; value of the TELESCOPE keyword
;
; IDL> fitsdir ,'/usr2/smith/*.fits',t='smith.txt', /NoTel
;
; PROCEDURE:
; FINDFILE (or FILE_SEARCH if since V5.5) is used to find the specified
; FITS files. The header of each file is read, and the selected
; keywords are extracted. The formatting is adjusted so that no value
; is truncated on display.
;
; SYSTEM VARIABLES:
; TEXTOPEN (called by FITSDIR) will automatically define the following
; non-standard system variables if they are not previously defined:
;
; DEFSYSV,'!TEXTOUT',1
; DEFSYSV,'!TEXTUNIT',0
;
; PROCEDURES USED:
; FDECOMP, FXMOVE, MRD_HREAD, REMCHAR, SPEC_DIR(),
; TEXTOPEN, TEXTCLOSE
; MODIFICATION HISTORY:
; Written, W. Landsman, HSTX February, 1993
; Converted to IDL V5.0 W. Landsman September 1997
; Search alternate keyword names W.Landsman October 1998
; Avoid integer truncation for NAXISi >32767 W. Landsman July 2000
; Don't leave open unit W. Landsman July 2000
; Added EXTEN keyword, work with compressed files, additional alternate
; keywords W. Landsman December 2000
; Don't assume floating pt. exposure time W. Landsman September 2001
; Major rewrite, KEYWORD & ALT*_KEYWORDS keywords, no truncation,
; /NOSIZE keyword W. Landsman, SSAI August 2002
; NAXIS* values must be integers W. Landsman SSAI June 2003
; NAME:
; FITS_READ
;*PURPOSE:
; To read a FITS file.
;
;*CATEGORY:
; INPUT/OUTPUT
;
;*CALLING SEQUENCE:
; FITS_READ, filename_or_fcb, data [,header, group_par]
;
;*INPUTS:
; FILENAME_OR_FCB - this parameter can be the FITS Control Block (FCB)
; returned by FITS_OPEN or the file name of the FITS file. If
; a file name is supplied, FITS_READ will open the file with
; FITS_OPEN and close the file with FITS_CLOSE before exiting.
; When multiple extensions are to be read from the file, it is
; more efficient for the user to call FITS_OPEN and leave the
; file open until all extensions are read.
;
;*OUTPUTS:
; DATA - data array. If /NOSCALE is specified, BSCALE and BZERO
; (if present in the header) will not be used to scale the data.
; If Keywords FIRST and LAST are used to read a portion of the
; data or the heap portion of an extension, no scaling is done
; and data is returned as a 1-D vector. The user can use the IDL
; function REFORM to convert the data to the correct dimensions
; if desired. If /DATA_ONLY is specified, no scaling is done.
; HEADER - FITS Header. If an extension is read, and the /NO_PDU keyword
; parameter is not supplied, the primary data unit header
; and the extension header will be combined. The header will have
; the form:
;
;
; BEGIN MAIN HEADER --------------------------------
;
; BEGIN EXTENSION HEADER ---------------------------
; 1. (Default=0, the first group)
;
;*OUTPUT KEYWORD PARAMETERS:
; ENUM - Output extension number that was read.
; MESSAGE = value: Output error message
;
;*NOTES:
; Determination or which extension to read.
; case 1: EXTEN_NO specified. EXTEN_NO will give the number of the
; extension to read. The primary data unit is refered
; to as extension 0. If EXTEN_NO is specified, XTENSION,
; EXTNAME, EXTVER, and EXTLEVEL parameters are ignored.
; case 2: if EXTEN_NO is not specified, the first extension
; with the specified XTENSION, EXTNAME, EXTVER, and
; EXTLEVEL will be read. If any of the 4 parameters
; are not specified, they will not be used in the search.
; Setting EXTLEVEL=0, EXTVER=0, EXTNAME='', or
; XTENSION='' is the same as not supplying them.
; case 3: if none of the keyword parameters, EXTEN_NO, XTENSION,
; EXTNAME, EXTVER, or EXTLEVEL are supplied. FITS_READ
; will read the next extension in the file. If the
; primary data unit (PDU), extension 0, is null, the
; first call to FITS_READ will read the first extension
; of the file.
;
; The only way to read a null PDU is to use EXTEN_NO = 0.
;
; If FIRST and LAST are specified, the data is returned without applying
; any scale factors (BSCALE and BZERO) and the data is returned in a
; 1-D vector. This will allow you to read any portion of a multiple
; dimension data set. Once returned, the IDL function REFORM can be
; used to place the correct dimensions on the data.
;
; IMPLICIT IMAGES: FITS_READ will construct an implicit image
; for cases where NAXIS=0 and the NPIX1, NPIX2, and PIXVALUE
; keywords are present. The output image will be:
; image = replicate(PIXVALUE,NPIX1,NPIX2)
;
;*EXAMPLES:
; Read the primary data unit of a FITS file, if it is null read the
; first extension:
; FITS_READ, 'myfile.fits', data, header
;
; Read the first two extensions of a FITS file and the extension with
; EXTNAME = 'FLUX' and EXTVER = 4
; FITS_OPEN, 'myfile.fits', fcb
; FITS_READ, fcb,data1, header2, exten_no = 1
; FITS_READ, fcb,data1, header2, exten_no = 2
; FITS_READ, fcb,data3, header3, extname='flux', extver=4
; FITS_CLOSE, fcb
;
; Read the sixth image in a data cube for the fourth extension.
;
; FITS_OPEN, 'myfile.fits', fcb
; image_number = 6
; ns = fcb.axis(0,4)
; nl = fcb.axis(1,4)
; i1 = (ns*nl)*(image_number-1)
; i2 = i2 + ns*nl-1
; FITS_READ,fcb,image,header,first=i1,last=i2
; image = reform(image,ns,nl,/overwrite)
; FITS_CLOSE
;
;*PROCEDURES USED:
; FITS_CLOSE, FITS_OPEN, IEEE_TO_HOST, IS_IEEE_BIG()
; SXADDPAR, SXDELPAR, SXPAR()
;*HISTORY:
; Written by: D. Lindler, August 1995
; Converted to IDL V5.0 W. Landsman September 1997
; Avoid use of !ERR W. Landsman August 1999
; Read unsigned datatypes, added /no_unsigned W. Landsman December 1999
; Don't call FITS_CLOSE unless fcb is defined W. Landsman January 2000
; Set BZERO = 0 for unsigned integer data W. Landsman January 2000
; Only call IEEE_TO_HOST if needed W. Landsman February 2000
; Ensure EXTEND keyword in primary header W. Landsman April 2001
; Don't erase ERROR message when closing file W. Landsman April 2002
; Assume at least V5.1 remove NANValue keyword W. Landsman November 2002
; Work with compress files (read file size from fcb),
; requires updated (Jan 2003) version of FITS_OPEN W. Landsman Jan 2003
; NAME:
; FITS_WRITE
;
;*PURPOSE:
; To write a FITS primary data unit or extension.
;
;*CATEGORY:
; INPUT/OUTPUT
;
;*CALLING SEQUENCE:
; FITS_WRITE, filename_or_fcb, data, [header_in]
;
;*INPUTS:
; FILENAME_OR_FCB: name of the output data file or the FITS control
; block returned by FITS_OPEN (called with the /WRITE or
; /APPEND) parameters.
;
;*OPTIONAL INPUTS:
; DATA: data array to write. If not supplied or set to a scalar, a
; null image is written.
; HEADER_IN: FITS header keyword. If not supplied, a minimal basic
; header will be created. Required FITS keywords, SIMPLE,
; BITPIX, XTENSION, NAXIS, ... are added by FITS_WRITE and
; do not need to be supplied with the header. If supplied,
; their values will be updated as necessary to reflect DATA.
;
;*INPUT KEYWORD PARAMETERS:
;
; XTENSION: type of extension to write (Default="IMAGE"). If not
; supplied, it will be taken from HEADER_IN. If not in either
; place, the default is "IMAGE". This parameter is ignored
; when writing the primary data unit.
; EXTNAME: EXTNAME for the extension. If not supplied, it will be taken
; from HEADER_IN. If not supplied and not in HEADER_IN, no
; EXTNAME will be written into the output extension.
; EXTVER: EXTVER for the extension. If not supplied, it will be taken
; from HEADER_IN. If not supplied and not in HEADER_IN, no
; EXTVER will be written into the output extension.
; EXTLEVEL: EXTLEVEL for the extension. If not supplied, it will be taken
; from HEADER_IN. If not supplied and not in HEADER_IN, no
; EXTLEVEL will be written into the output extension.
; /NO_ABORT: Set to return to calling program instead of a RETALL
; when an I/O error is encountered. If set, the routine will
; return a non-null string (containing the error message) in the
; keyword MESSAGE. (For backward compatibility, the obsolete
; system variable !ERR is also set to -1 in case of an error.)
; If /NO_ABORT not set, then FITS_WRITE will print the message and
; issue a RETALL
; /NO_DATA: Set if you only want FITS_WRITE to write a header. The
; header supplied will be written without modification and
; the user is expected to write the data using WRITEU to unit
; FCB.UNIT. When FITS_WRITE is called with /NO_DATA, the user is
; responsible for the validity of the header, and must write
; the correct amount and format of the data. When FITS_WRITE
; is used in this fashion, it will pad the data from a previously
; written extension to 2880 blocks before writting the header.
;
;*OUTPUT KEYWORD PARAMETERS:
; MESSAGE: value of the error message for use with /NO_ABORT
; HEADER: actual output header written to the FITS file.
;
;*NOTES:
; If the first call to FITS_WRITE is an extension, FITS_WRITE will
; automatically write a null image as the primary data unit.
;
; Keywords and history in the input header will be properly separated
; into the primary data unit and extension portions when constructing
; the output header (See FITS_READ for information on the internal
; Header format which separates the extension and PDU header portions).
;
;*EXAMPLES:
; Write an IDL variable to a FITS file with the minimal required header.
; FITS_WRITE,'newfile.fits',ARRAY
;
; Write the same array as an image extension, with a null Primary data
; unit.
; FITS_WRITE,'newfile.fits',ARRAY,xtension='IMAGE'
;
; Write 4 additional image extensions to the same file.
; FITS_OPEN,'newfile.fits',fcb
; FITS_WRITE,fcb,data1,extname='FLUX',extver=1
; FITS_WRITE,fcb,err1,extname'ERR',extver=1
; FITS_WRITE,fcb,data2,extname='FLUX',extver=2
; FITS_WRITE,fcb,err2,extname='ERR',extver=2
; FITS_CLOSE,FCB
;
;*PROCEDURES USED:
; FITS_OPEN, SXADDPAR, SXDELPAR, SXPAR()
;*HISTORY:
; Written by: D. Lindler August, 1995
; Work for variable length extensions W. Landsman August 1997
; Converted to IDL V5.0 W. Landsman September 1997
; PCOUNT and GCOUNT added for IMAGE extensions J. Graham October 1999
; Write unsigned data types W. Landsman December 1999
; Pad data area with zeros not blanks W. McCann/W. Landsman October 2000
; Return Message='' to signal normal operation W. Landsman Nov. 2000
; Ensure that required extension table keywords are in proper order
; W.V. Dixon/W. Landsman March 2001
; Assume since V5.1, remove NaNValue keyword W. Landsman Nov. 2002
; NAME:
; FM_UNRED
; PURPOSE:
; Deredden a flux vector using the Fitzpatrick (1999) parameterization
; EXPLANATION:
; The R-dependent Galactic extinction curve is that of Fitzpatrick & Massa
; (Fitzpatrick, 1999, PASP, 111, 63; astro-ph/9809387 ).
; Parameterization is valid from the IR to the far-UV (3.5 microns to 0.1
; microns). UV extinction curve is extrapolated down to 912 Angstroms.
;
; CALLING SEQUENCE:
; FM_UNRED, wave, flux, ebv, [ funred, R_V = , /LMC2, /AVGLMC, ExtCurve=
; gamma =, x0=, c1=, c2=, c3=, c4= ]
; INPUT:
; WAVE - wavelength vector (Angstroms)
; FLUX - calibrated flux vector, same number of elements as WAVE
; If only 3 parameters are supplied, then this vector will
; updated on output to contain the dereddened flux.
; EBV - color excess E(B-V), scalar. If a negative EBV is supplied,
; then fluxes will be reddened rather than deredenned.
;
; OUTPUT:
; FUNRED - unreddened flux vector, same units and number of elements
; as FLUX
;
; OPTIONAL INPUT KEYWORDS
; R_V - scalar specifying the ratio of total to selective extinction
; R(V) = A(V) / E(B - V). If not specified, then R = 3.1
; Extreme values of R(V) range from 2.3 to 5.3
;
; /AVGLMC - if set, then the default fit parameters c1,c2,c3,c4,gamma,x0
; are set to the average values determined for reddening in the
; general Large Magellanic Cloud (LMC) field by Misselt et al.
; (1999, ApJ, 515, 128)
; /LMC2 - if set, then the fit parameters are set to the values determined
; for the LMC2 field (including 30 Dor) by Misselt et al.
; Note that neither /AVGLMC or /LMC2 will alter the default value
; of R_V which is poorly known for the LMC.
;
; The following five input keyword parameters allow the user to customize
; the adopted extinction curve. For example, see Clayton et al. (2003,
; ApJ, 588, 871) for examples of these parameters in different interstellar
; environments.
;
; x0 - Centroid of 2200 A bump in microns (default = 4.596)
; gamma - Width of 2200 A bump in microns (default =0.99)
; c3 - Strength of the 2200 A bump (default = 3.23)
; c4 - FUV curvature (default = 0.41)
; c2 - Slope of the linear UV extinction component
; (default = -0.824 + 4.717/R)
; c1 - Intercept of the linear UV extinction component
; (default = 2.030 - 3.007*c2
;
; OPTIONAL OUTPUT KEYWORD:
; ExtCurve - Returns the E(wave-V)/E(B-V) extinction curve, interpolated
; onto the input wavelength vector
;
; EXAMPLE:
; Determine how a flat spectrum (in wavelength) between 1200 A and 3200 A
; is altered by a reddening of E(B-V) = 0.1. Assume an "average"
; reddening for the diffuse interstellar medium (R(V) = 3.1)
;
; IDL> w = 1200 + findgen(40)*50 ;Create a wavelength vector
; IDL> f = w*0 + 1 ;Create a "flat" flux vector
; IDL> fm_unred, w, f, -0.1, fnew ;Redden (negative E(B-V)) flux vector
; IDL> plot,w,fnew
;
; NOTES:
; (1) The following comparisons between the FM curve and that of Cardelli,
; Clayton, & Mathis (1989), (see ccm_unred.pro):
; (a) - In the UV, the FM and CCM curves are similar for R < 4.0, but
; diverge for larger R
; (b) - In the optical region, the FM more closely matches the
; monochromatic extinction, especially near the R band.
; (2) Many sightlines with peculiar ultraviolet interstellar extinction
; can be represented with the FM curve, if the proper value of
; R(V) is supplied.
; (3) Use the 4 parameter calling sequence if you wish to save the
; original flux vector.
; PROCEDURE CALLS:
; CSPLINE(), POLY()
; REVISION HISTORY:
; Written W. Landsman Raytheon STX October, 1998
; Based on FMRCurve by E. Fitzpatrick (Villanova)
; Added /LMC2 and /AVGLMC keywords, W. Landsman August 2000
; Added ExtCurve keyword, J. Wm. Parker August 2000
;
; NAME:
; FORPRINT
; PURPOSE:
; Print a set of vectors by looping over each index value.
;
; EXPLANATION:
; If W and F are equal length vectors, then the statement
; IDL> forprint, w, f
; is equivalent to
; IDL> for i = 0L, N_elements(w)-1 do print,w[i],f[i]
;
; CALLING SEQUENCE:
; forprint, v1,[ v2, v3, v4,....v18, FORMAT = , TEXTOUT = ,STARTLINE =,
; NUMLINE =, /SILENT, COMMENT= ]
;
; INPUTS:
; V1,V2,...V18 - Arbitary IDL vectors. If the vectors are not of
; equal length then the number of rows printed will be equal
; to the length of the smallest vector. Up to 18 vectors
; can be supplied.
;
; OPTIONAL KEYWORD INPUTS:
;
; TEXTOUT - Controls print output device, defaults to !TEXTOUT
;
; textout=1 TERMINAL using /more option if available
; textout=2 TERMINAL without /more option
; textout=3 file 'forprint.prt'
; textout=4 file 'laser.tmp'
; textout=5 user must open file
; textout = filename (default extension of .prt)
; textout=7 Append to .prt file if it exists
;
; COMMENT - String to write as the first line of output file if
; TEXTOUT > 2. By default, FORPRINT will write a time stamp
; on the first line. Use /NOCOMMENT if you don't want FORPRINT
; to write anything in the output file.
; FORMAT - Scalar format string as in the PRINT procedure. The use
; of outer parenthesis is optional. Ex. - format="(F10.3,I7)"
; This program will automatically remove a leading "$" from
; incoming format statments. Ex. - "$(I4)" would become "(I4)".
; If omitted, then IDL default formats are used.
; /NOCOMMENT - Set this keyword if you don't want any comment line
; line written as the first line in a harcopy output file.
; /SILENT - Normally, with a hardcopy output (TEXTOUT > 2), FORPRINT will
; print an informational message. If the SILENT keyword
; is set and non-zero, then this message is suppressed.
; STARTLINE - Integer scalar specifying the first line in the arrays
; to print. Default is STARTLINE = 1, i.e. start at the
; beginning of the arrays.
; OUTPUTS:
; None
; SYSTEM VARIABLES:
; If keyword TEXTOUT is not used, the default is the nonstandard
; keyword !TEXTOUT. If you want to use FORPRINT to write more than
; once to the same file, or use a different file name then set
; TEXTOUT=5, and open and close then file yourself (see documentation
; of TEXTOPEN for more info).
;
; One way to add the non-standard system variables !TEXTOUT and !TEXTUNIT
; is to use the procedure ASTROLIB
; EXAMPLE:
; Suppose W,F, and E are the wavelength, flux, and epsilon vectors for
; a spectrum. Print these values to a file 'output.dat' in a nice
; format.
;
; IDL> fmt = '(F10.3,1PE12.2,I7)'
; IDL> forprint, F = fmt, w, f, e, TEXT = 'output.dat'
;
; PROCEDURES CALLED:
; TEXTOPEN, TEXTCLOSE
; REVISION HISTORY:
; Written W. Landsman April, 1989
; Keywords textout and format added, J. Isensee, July, 1990
; Made use of parenthesis in FORMAT optional W. Landsman May 1992
; Added STARTLINE keyword W. Landsman November 1992
; Set up so can handle 18 input vectors. J. Isensee, HSTX Corp. July 1993
; Handle string value of TEXTOUT W. Landsman, HSTX September 1993
; Added NUMLINE keyword W. Landsman, HSTX February 1996
; Added SILENT keyword W. Landsman, RSTX, April 1998
; Converted to IDL V5.0 W. Landsman, RSTX, April, 1998
; Much faster printing to a file W. Landsman, RITSS, August, 2001
; Use SIZE(/TNAME) instead of DATATYPE() W. Landsman SSAI October 2001
; Fix skipping of first line bug introduced Aug 2001 W. Landsman Nov2001
; Added /NOCOMMENT keyword, the SILENT keyword now controls only
; the display of informational messages. W. Landsman June 2002
; NAME:
; FTPUT
; PURPOSE:
; Procedure to add or update a field in an FITS ASCII table
;
; CALLING SEQUENCE:
; FTPUT, htab, tab, field, row, values, [ nulls ]
;
; INPUTS:
; htab - FITS ASCII table header string array
; tab - FITS ASCII table array (e.g. as read by READFITS)
; field - string field name or integer field number
; row - either a non-negative integer scalar giving starting row to
; update, or a non-negative integer vector specifying rows to
; update. FTPUT will append a new row to a table if the value
; of 'row' exceeds the number of rows in the tab array
; values - value(s) to add or update. If row is a vector
; then values must contain the same number of elements.
;
; OPTIONAL INPUT:
; nulls - null value flag of same length as values.
; It should be set to 1 at null value positions
; and 0 elsewhere.
;
; OUTPUTS:
; htab,tab will be updated as specified.
;
; EXAMPLE:
; One has a NAME and RA and Dec vectors for 500 stars with formats A6,
; F9.5 and F9.5 respectively. Write this information to an ASCII table
; named 'star.fits'.
;
; IDL> FTCREATE,24,500,h,tab ;Create table header and (empty) data
; IDL> FTADDCOL,h,tab,'RA',8,'F9.5','DEGREES' ;Explicity define the
; IDL> FTADDCOL,h,tab,'DEC',8,'F9.5','DEGREES' ;RA and Dec columns
; IDL> FTPUT,h,tab,'RA',0,ra ;Insert RA vector into table
; IDL> FTPUT,h,tab,'DEC',0,dec ;Insert DEC vector into table
; IDL> FTPUT, h,tab, 'NAME',0,name ;Insert NAME vector with default
; IDL> WRITEFITS,'stars.fits',tab,h ;Write to a file
;
; Note that (1) explicit formatting has been supplied for the (numeric)
; RA and Dec vectors, but was not needed for the NAME vector, (2) A width
; of 24 was supplied in FTCREATE based on the expected formats (6+9+9),
; though the FT* will adjust this value as necessary, and (3) WRITEFITS
; will create a minimal primary header
; NOTES:
; (1) If the specified field is not already in the table, then FTPUT will
; create a new column for that field using default formatting. However,
; FTADDCOL should be called prior to FTPUT for explicit formatting.
;
; PROCEDURES CALLED
; FSTRING(), FTADDCOL, FTINFO, FTSIZE, SXADDPAR, SXPAR()
; HISTORY:
; version 1 D. Lindler July, 1987
; Allow E format W. Landsman March 1992
; Write in F format if E format will overflow April 1994
; Update documentation W. Landsman January 1996
; Allow 1 element vector W. Landsman March 1996
; Adjust string length to maximum of input string array June 1997
; Work for more than 32767 elements August 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Updated call to the new FTINFO W. Landsman May 2000
; Fix case where header does not have any columns yet W.Landsman Sep 2002
; NAME:
; FXADDPAR
; Purpose :
; Add or modify a parameter in a FITS header array.
; Explanation :
; This version of FXADDPAR will write string values longer than 68
; characters using the FITS continuation convention described at
; http://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/ofwg_recomm/r13.html
; Use :
; FXADDPAR, HEADER, NAME, VALUE, COMMENT
; Inputs :
; HEADER = String array containing FITS header. The maximum string
; length must be equal to 80. If not defined, then FXADDPAR
; will create an empty FITS header array.
;
; NAME = Name of parameter. If NAME is already in the header the
; value and possibly comment fields are modified. Otherwise a
; new record is added to the header. If NAME is equal to
; either "COMMENT" or "HISTORY" then the value will be added to
; the record without replacement. In this case the comment
; parameter is ignored.
;
; VALUE = Value for parameter. The value expression must be of the
; correct type, e.g. integer, floating or string.
; String values of 'T' or 'F' are considered logical
; values. If the value is a string and is "long"
; (more than 69 characters), then it may be continued
; over more than one line using the OGIP CONTINUE
; standard.
;
; Opt. Inputs :
; COMMENT = String field. The '/' is added by this routine. Added
; starting in position 31. If not supplied, or set equal to ''
; (the null string), then any previous comment field in the
; header for that keyword is retained (when found).
; Outputs :
; HEADER = Updated header array.
; Opt. Outputs:
; None.
; Keywords :
; BEFORE = Keyword string name. The parameter will be placed before the
; location of this keyword. For example, if BEFORE='HISTORY'
; then the parameter will be placed before the first history
; location. This applies only when adding a new keyword;
; keywords already in the header are kept in the same position.
;
; AFTER = Same as BEFORE, but the parameter will be placed after the
; location of this keyword. This keyword takes precedence over
; BEFORE.
;
; FORMAT = Specifies FORTRAN-like format for parameter, e.g. "F7.3". A
; scalar string should be used. For complex numbers the format
; should be defined so that it can be applied separately to the
; real and imaginary parts.
;
; /NOCONTINUE = By default, FXADDPAR will break strings longer than 68
; characters into multiple lines using the continuation
; convention. If this keyword is set, then the line will
; instead be truncated to 68 characters. This was the default
; behaviour of FXADDPAR prior to December 1999.
; Calls :
; DETABIFY(), FXPAR(), FXPARPOS()
; Common :
; None.
; Restrictions:
; Warning -- Parameters and names are not checked against valid FITS
; parameter names, values and types.
;
; The required FITS keywords SIMPLE (or XTENSION), BITPIX, NAXIS, NAXIS1,
; NAXIS2, etc., must be entered in order. The actual values of these
; keywords are not checked for legality and consistency, however.
;
; Side effects:
; All HISTORY records are inserted in order at the end of the header.
;
; All COMMENT records are also inserted in order at the end of the
; header, but before the HISTORY records. The BEFORE and AFTER keywords
; can override this.
;
; All records with no keyword (blank) are inserted in order at the end of
; the header, but before the COMMENT and HISTORY records. The BEFORE and
; AFTER keywords can override this.
;
; All other records are inserted before any of the HISTORY, COMMENT, or
; "blank" records. The BEFORE and AFTER keywords can override this.
;
; String values longer than 68 characters will be split into multiple
; lines using the OGIP CONTINUE convention, unless the /NOCONTINUE keyword
; is set. For a description of the CONTINUE convention see
; http://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/ofwg_recomm/r13.htm
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; William Thompson, Jan 1992, from SXADDPAR by D. Lindler and J. Isensee.
; Differences include:
;
; * LOCATION parameter replaced with keywords BEFORE and AFTER.
; * Support for COMMENT and "blank" FITS keywords.
; * Better support for standard FITS formatting of string and
; complex values.
; * Built-in knowledge of the proper position of required
; keywords in FITS (although not necessarily SDAS/Geis) primary
; headers, and in TABLE and BINTABLE extension headers.
;
; William Thompson, May 1992, fixed bug when extending length of header,
; and new record is COMMENT, HISTORY, or blank.
; Written :
; William Thompson, GSFC, January 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 5 September 1997
; Fixed bug replacing strings that contain "/" character--it
; interpreted the following characters as a comment.
; Version 3, Craig Markwardt, GSFC, December 1997
; Allow long values to extend over multiple lines
; Version 4, D. Lindler, March 2000, modified to use capital E instead
; of a lower case e for exponential format.
; Version 4.1 W. Landsman April 2000, make user-supplied format uppercase
; Version 4.2 W. Landsman July 2002, positioning of EXTEND keyword
; Version :
; Version 4.2, July 2002
; NAME:
; FXBADDCOL
; Purpose :
; Adds a column to a binary table extension.
; Explanation :
; Modify a basic FITS binary table extension (BINTABLE) header array to
; define a column.
; Use :
; FXBADDCOL, INDEX, HEADER, ARRAY [, TTYPE [, COMMENT ]]
; Inputs :
; HEADER = String array containing FITS extension header.
; ARRAY = IDL variable used to determine the data size and type
; associated with the column. If the column is defined as
; containing variable length arrays, then ARRAY must be of the
; maximum size to be stored in the column.
; Opt. Inputs :
; TTYPE = Column label.
; COMMENT = Comment for TTYPE
; Outputs :
; INDEX = Index (1-999) of the created column.
; HEADER = The header is modified to reflect the added column.
; Opt. Outputs:
; None.
; Keywords :
; VARIABLE= If set, then the column is defined to contain pointers to
; variable length arrays in the heap area.
; DCOMPLEX= If set, and ARRAY is complex, with the first dimension being
; two (real and imaginary parts), then the column is defined as
; double-precision complex (type "M"). This keyword is
; only needed prior to IDL Version 4.0, when the double
; double complex datatype was unavailable in IDL
; BIT = If passed, and ARRAY is of type byte, then the column is
; defined as containg bit mask arrays (type "X"), with the
; value of BIT being equal to the number of mask bits.
; LOGICAL = If set, and array is of type byte, then the column is defined
; as containing logical arrays (type "L").
; NO_TDIM = If set, then the TDIMn keyword is not written out to the
; header. No TDIMn keywords are written for columns containing
; variable length arrays.
; TUNIT = If passed, then corresponding keyword is added to header.
; TSCAL = Same.
; TZERO = Same.
; TNULL = Same.
; TDISP = Same.
; TDMIN = Same.
; TDMAX = Same.
; TDESC = Same.
; TCUNI = Same.
; TROTA = Same.
; TRPIX = Same.
; TRVAL = Same.
; TDELT = Same.
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBADDCOL, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; FXADDPAR, FXPAR
; Common :
; None.
; Restrictions:
; Warning: No checking is done of any of the parameters defining the
; values of optional FITS keywords.
;
; FXBHMAKE must first be called to initialize the header.
;
; If ARRAY is of type character, then it must be of the maximum length
; expected for this column. If a character string array, then the
; largest string in the array is used to determine the maximum length.
;
; The DCOMPLEX keyword is ignored if ARRAY is not double-precision.
; ARRAY must also have a first dimension of two representing the real and
; imaginary parts.
;
; The BIT and LOGICAL keywords are ignored if ARRAY is not of type byte.
; BIT takes precedence over LOGICAL.
;
; Side effects:
; If the data array is multidimensional, then a TDIM keyword is added to
; the header, unless either NO_TDIM or VARIABLE is set.
;
; No TDIMn keywords are written out for bit arrays (format 'X'), since
; the dimensions would refer to bits, not bytes.
;
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; William Thompson, Jan 1992.
; W. Thompson, Feb 1992, changed from function to procedure.
; W. Thompson, Feb 1992, modified to support variable length arrays.
; Written :
; William Thompson, GSFC, January 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 31 May 1994
; Added ERRMSG keyword.
; Version 3, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 4, William Thompson, GSFC, 30 December 1994
; Added keyword TCUNI.
; Version 5, Wayne Landsman, GSFC, 12 Aug 1997
; Recognize double complex IDL datatype
; Version :
; Version 5, 12 Aug 1997
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; FXBCREATE
; Purpose :
; Open a new binary table at the end of a FITS file.
; Explanation :
; Write a binary table extension header to the end of a disk FITS file,
; and leave it open to receive the data.
;
; The FITS file is opened, and the pointer is positioned just after the
; last 2880 byte record. Then the binary header is appended. Calls to
; FXBWRITE will append the binary data to this file, and then FXBFINISH
; will close the file.
;
; Use :
; FXBCREATE, UNIT, FILENAME, HEADER
; Inputs :
; FILENAME = Name of FITS file to be opened.
; HEADER = String array containing the FITS binary table extension
; header.
; Opt. Inputs :
; None.
; Outputs :
; UNIT = Logical unit number of the opened file.
; EXTENSION= Extension number of newly created extension.
; Opt. Outputs:
; None.
; Keywords :
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBCREATE, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; FXADDPAR, FXBFINDLUN, FXBPARSE, FXFINDEND
; Common :
; Uses common block FXBINTABLE--see "fxbintable.pro" for more
; information.
; Restrictions:
; The primary FITS data unit must already be written to a file. The
; binary table extension header must already be defined (FXBHMAKE), and
; must match the data that will be written to the file.
; Side effects:
; None.
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; W. Thompson, Jan 1992, based on WRITEFITS by J. Woffard and W. Landsman.
; W. Thompson, Feb 1992, changed from function to procedure.
; W. Thompson, Feb 1992, removed all references to temporary files.
; Written :
; William Thompson, GSFC, January 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 21 July 1993.
; Fixed bug with variable length arrays.
; Version 3, William Thompson, GSFC, 21 June 1994
; Added ERRMSG keyword.
; Version 4, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 5, Antony Bird, Southampton, 25 June 1997
; Modified to allow very long tables
; Version :
; Version 5, 25 June 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Added EXTENSION parameter, C. Markwardt 1999 Jul 15
; More efficient zeroing of file, C. Markwardt, 26 Feb 2001
; NAME:
; FXBOPEN
; Purpose :
; Open binary table extension in a disk FITS file for reading.
; Explanation :
; Opens a binary table extension in a disk FITS file for reading. The
; columns are then read using FXBREAD, and the file is closed when done
; with FXBCLOSE.
; Use :
; FXBOPEN, UNIT, FILENAME, EXTENSION [, HEADER ]
; Inputs :
; FILENAME = Name of FITS file to be opened. Optional
; extension *number* may be specified, in either of
; the following formats (using the FTOOLS
; convention): FILENAME[EXT] or FILENAME+EXT, where
; EXT is 1 or higher. Such an extension
; specification takes priority over EXTENSION.
;
; EXTENSION = Either the number of the FITS extension, starting with the
; first extension after the primary data unit being one; or a
; character string containing the value of EXTNAME to search
; for.
; Opt. Inputs :
; None.
; Outputs :
; UNIT = Logical unit number of the opened file.
; Opt. Outputs:
; HEADER = String array containing the FITS binary table extension
; header.
; Keywords :
; NO_TDIM = If set, then any TDIMn keywords found in the header are
; ignored.
;
; ACCESS = A scalar string describing access privileges as
; one of READ ('R') or UPDATE ('RW').
; DEFAULT: 'R'
;
; REOPEN = If set, UNIT must be an already-opened file unit.
; FXBOPEN will treat the file as a FITS file.
;
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBOPEN, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; FXBFINDLUN, FXBPARSE, FXHREAD, FXPAR
; Common :
; Uses common block FXBINTABLE--see "fxbintable.pro" for more
; information.
; Restrictions:
; The file must be a valid FITS file.
; Side effects:
; None.
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; W. Thompson, Feb 1992, based on READFITS by J. Woffard and W. Landsman.
; W. Thompson, Feb 1992, changed from function to procedure.
; W. Thompson, June 1992, fixed up error handling.
; Written :
; William Thompson, GSFC, February 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 27 May 1994
; Added ERRMSG keyword.
; Version 3, William Thompson, GSFC, 21 June 1994
; Extended ERRMSG to call to FXBPARSE
; Version 4, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 4, 23 June 1994
; Converted to IDL V5.0 W. Landsman September 1997
;
; Added ACCESS, REOPEN keywords, and FXFILTER package, CM 1999 Feb 03
; Added FILENAME[EXT] and FILENAME+EXT extension parsing, CM 1999 Jun 28
; Some general tidying, CM 1999 Nov 18
;
; NAME:
; FXBPARSE
; Purpose :
; Parse the binary table extension header.
; Explanation :
; Parses the binary table extension header, and store the information
; about the format of the binary table in the FXBINTABLE common
; block--called from FXBCREATE and FXBOPEN.
; Use :
; FXBPARSE, ILUN, UNIT, HEADER
; Inputs :
; ILUN = Index into the arrays in the FXBINTABLE common block.
; HEADER = FITS binary table extension header.
; Opt. Inputs :
; None.
; Outputs :
; None.
; Opt. Outputs:
; None.
; Keywords :
; NO_TDIM = If set, then any TDIMn keywords found in the header are
; ignored.
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBPARSE, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; FXBFIND, FXBTDIM, FXBTFORM, FXPAR
; Common :
; Uses common block FXBINTABLE--see "fxbintable.pro" for more
; information.
; Restrictions:
; None.
; Side effects:
; Any TDIMn keywords found for bit arrays (format 'X') are ignored, since
; the dimensions would refer to bits, not bytes.
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; William Thompson, Feb. 1992.
; William Thompson, Jan. 1993, modified for renamed FXBTFORM and FXBTDIM.
; Written :
; William Thompson, GSFC, February 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 21 June 1994
; Added ERRMSG keyword.
; Version 3, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 4, Michael Schubnell, University of Michigan, 22 May 1996
; Change N_DIMS from short to long integer.
; Version 5, W. Landsman, GSFC, 12 Aug 1997
; Use double complex datatype, if needed
; Version 6, W. Landsman GSFC 30 Aug 1997
; Version :
; Version 6, 31 Aug 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Optimized FXPAR; call FXBFIND for speed, CM 1999 Nov 18
; Modify DHEAP(ILUN) when opening table now, CM 2000 Feb 22
; NAME:
; FXBREAD
; Purpose :
; Read a data array from a disk FITS binary table file.
; Explanation :
; Each call to FXBREAD will read the data from one column and one row
; from the FITS data file, which should already have been opened by
; FXBOPEN. One needs to call this routine for every column and every row
; in the binary table. FXBCLOSE will then close the FITS data file.
; Use :
; FXBREAD, UNIT, DATA, COL [, ROW ]
; Inputs :
; UNIT = Logical unit number corresponding to the file containing the
; binary table.
; COL = Column in the binary table to read data from, either as a
; character string containing a column label (TTYPE), or as a
; numerical column index starting from column one.
; Opt. Inputs :
; ROW = Either row number in the binary table to read data from,
; starting from row one, or a two element array containing a
; range of row numbers to read. If not passed, then the entire
; column is read in.
;
; Row must be passed for variable length arrays.
;
; Outputs :
; DATA = IDL data array to be read from the file.
; Opt. Outputs:
; None.
; Keywords :
; NOSCALE = If set, then the output data will not be scaled using the
; optional TSCAL and TZERO keywords in the FITS header.
; Default is to scale.
; NOIEEE = If set, then the output data is not byte-swapped to
; machine order. NOIEEE implies NOSCALE.
; Default is to perform the byte-swap.
; VIRTUAL = If set, and COL is passed as a name rather than a number,
; then if the program can't find a column with that name, it
; will then look for a keyword with that name in the header.
; Such a keyword would then act as a "virtual column", with the
; same value for every row.
; DIMENSIONS = Vector array containing the dimensions to be used to read
; in the data. Bypasses any dimensioning information stored in
; the header. Ignored for bit arrays. If the data type is
; double-precision complex, then an extra dimension of 2 is
; prepended to the dimensions passed by the user.
; NANVALUE= Value signalling data dropout. All points corresponding to
; IEEE NaN (not-a-number) are converted to this number.
; Ignored unless DATA is of type float, double-precision or
; complex.
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBREAD, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; IEEE_TO_HOST, FXPAR, WHERE_NEGZERO, WHERENAN
; Common :
; Uses common block FXBINTABLE--see "fxbintable.pro" for more
; information.
; Restrictions:
; The binary table file must have been opened with FXBOPEN.
;
; The data must be consistent with the column definition in the binary
; table header.
;
; The row number must be consistent with the number of rows stored in the
; binary table header.
;
; The number of elements implied by the dimensions keyword must not
; exceed the number of elements stored in the file.
;
; Side effects:
; If the DIMENSIONS keyword is used, then the number of data points read
; in may be less than the number of points stored in the table.
;
; If there are no elements to read in (the number of elements is zero),
; then the program sets !ERR to -1, and DATA is unmodified.
;
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; W. Thompson, Jan 1992.
; W. Thompson, Feb 1992, modified to support variable length arrays.
; W. Thompson, Jun 1992, modified way that row ranges are read in. No
; longer works reiteratively.
; W. Thompson, Jun 1992, fixed bug where NANVALUE would be modified by
; TSCAL and TZERO keywords.
; W. Thompson, Jun 1992, fixed bug when reading character strings.
; Treats dimensions better when reading multiple
; rows.
; Written :
; William Thompson, GSFC, January 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 30 June 1993.
; Added overwrite keyword to REFORM call to speed up.
; Version 3, William Thompson, GSFC, 21 July 1993.
; Fixed bug with variable length arrays.
; Version 4, William Thompson, GSFC, 29 October 1993.
; Added error message for not finding column by name.
; Version 5, William Thompson, GSFC, 31 May 1994
; Added ERRMSG keyword.
; Version 6, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 7, William Thompson, GSFC, 29 December 1994
; Fixed bug where single element dimensions were lost.
; Version 8, William Thompson, GSFC, 20 March 1995
; Fixed bug introduced in version 7.
; Version 9, Wayne Landsman, GSFC, 3 July 1996
; Fixed bug involving use of virtual keyword.
; Version 10, William Thompson, GSFC, 31-Jan-1997
; Added call to WHERE_NEGZERO.
; Version 11, Wayne Landsman, GSFC, 12 Aug, 1997
; Use IDL dcomplex datatype if needed
; Version 12, Wayne Landmsan, GSFC, 20 Feb, 1998
; Remove call to WHERE_NEGZERO (now part of IEEE_TO_HOST)
; Version 13, 18 Nov 1999, CM, Add NOIEEE keyword
; Version 14, 21 Aug 2000, William Thompson, GSFC
; Catch I/O errors
; Version :
; Version 14, 21 Aug 2000
; NAME:
; FXBREADM
; PURPOSE:
; Read multiple columns/rows from a disk FITS binary table file.
; EXPLANATION :
; A call to FXBREADM will read data from multiple rows and
; multiple columns in a single procedure call. Up to forty-nine
; columns may be read in a single pass; the number of rows is
; limited essentially by available memory. The file should have
; already been opened with FXBOPEN. FXBREADM optimizes reading
; multiple columns by first reading a large chunk of data from
; the FITS file directly, and then slicing the data into columns
; within memory. FXBREADM can read variable-length arrays (see
; below).
;
; The number of columns is limited to 49 if data are passed by
; positional argument. However, this limitation can be overcome
; by having FXBREADM return the data in an array of pointers or
; handles. The user should set the PASS_METHOD keyword to
; 'POINTER' or 'HANDLE' as appropriate, and an array of pointers
; to the data will be returned in the POINTERS keyword. The
; user is responsible for freeing the pointers; however,
; FXBREADM will reuse any pointers or handles passed into the
; procedure, and hence any pointed-to data will be destroyed.
;
; FXBREADM can also read variable-length columns from FITS
; binary tables. Since such data is not of a fixed size, it is
; returned as a structure. The structure has the following
; elements:
;
; VARICOL: ;; Flag: variable length column (= 1)
; N_ELEMENTS: ;; Total number of elements returned
; TYPE: ;; IDL data type code (integer)
; N_ROWS: ;; Number of rows read from table (integer)
; INDICES: ;; Indices of each row's data (integer array)
; DATA: ;; Raw data elements (variable type array)
;
; In order to gain access to the Ith row's data, one should
; examine DATA(INDICES(I):INDICES(I+1)-1), which is similar in
; construct to the REVERSE_INDICES keyword of the HISTOGRAM
; function.
;
; CALLING SEQUENCE:
; FXBREADM, UNIT, COL, DATA1, [ DATA2, ... DATA48, ROW=, BUFFERSIZE = ]
; /NOIEEE, /NOSCALE, /VIRTUAL, NANVALUE=, PASS_METHOD = POINTERS=,
; ERRMSG = , WARNMSG = , STATUS = ]
;
; INPUT PARAMETERS :
; UNIT = Logical unit number corresponding to the file containing the
; binary table.
; COL = An array of columns in the binary table to read data
; from, either as character strings containing column
; labels (TTYPE), or as numerical column indices
; starting from column one.
; Outputs :
; DATA1, DATA2...DATA48 = A named variable to accept the data values, one
; for each column. The columns are stored in order of the
; list in COL. If the read operation fails for a
; particular column, then the corresponding output Dn
; variable is not altered. See the STATUS keyword.
; Ignored if PASS_METHOD is 'POINTER' or 'HANDLE'
;
; OPTIONAL INPUT KEYWORDS:
; ROW = Either row number in the binary table to read data from,
; starting from row one, or a two element array containing a
; range of row numbers to read. If not passed, then the entire
; column is read in.
; /NOIEEE = If set, then then IEEE floating point data will not
; be converted to the host floating point format (and
; this by definition implies NOSCALE). The user is
; responsible for their own floating point conversion.
; /NOSCALE = If set, then the output data will not be scaled using the
; optional TSCAL and TZERO keywords in the FITS header.
; Default is to scale.
; VIRTUAL = If set, and COL is passed as a name rather than a number,
; then if the program can't find a column with that name, it
; will then look for a keyword with that name in the header.
; Such a keyword would then act as a "virtual column", with the
; same value for every row.
; DIMENSIONS = FXBREADM ignores this keyword. It is here for
; compatibility only.
; NANVALUE= Value signalling data dropout. All points corresponding to
; IEEE NaN (not-a-number) are converted to this number.
; Ignored unless DATA is of type float, double-precision or
; complex.
; PASS_METHOD = A scalar string indicating method of passing
; data from FXBREADM. One of 'ARGUMENT' (indicating
; pass by positional argument), 'POINTER' (indicating
; passing an array of pointers by the POINTERS
; keyword), or 'HANDLE' (indicating passing an array
; of handles by the POINTERS keyword).
; Default: 'ARGUMENT'
; POINTERS = If PASS_METHOD is 'POINTER' then an array of IDL
; pointers is returned in this keyword, one for each
; requested column. If PASS_METHOD is 'HANDLE' then
; an array of IDL handles is returned instead. Any
; handles or pointers passed into FXBREADM will have
; their pointed-to data destroyed. Ultimately the
; user is responsible for deallocating pointers and
; handles.
; BUFFERSIZE = Raw data are transferred from the file in chunks
; to conserve memory. This is the size in bytes of
; each chunk. If a value of zero is given, then all
; of the data are transferred in one pass. Default is
; 32768 (32 kB).
; OPTIONAL OUTPUT KEYWORDS:
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBREAD, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
; WARNMSG = Messages which are considered to be non-fatal
; "warnings" are returned in this output string.
; Note that if some but not all columns are
; unreadable, this is considered to be non-fatal.
; STATUS = An output array containing the status for each
; column read, 1 meaning success and 0 meaning failure.
;
; Calls :
; IEEE_TO_HOST, FXPAR(), WHERENAN()
; Common :
; Uses common block FXBINTABLE--see "fxbintable.pro" for more
; information.
; Restrictions:
; The binary table file must have been opened with FXBOPEN.
;
; The data must be consistent with the column definition in the binary
; table header.
;
; The row number must be consistent with the number of rows stored in the
; binary table header.
;
; Generaly speaking, FXBREADM will be faster than iterative
; calls to FXBREAD when (a) a large number of columns is to be
; read or (b) the size in bytes of each cell is small, so that
; the overhead of the FOR loop in FXBREAD becomes significant.
;
; SIDE EFFECTS:
; If there are no elements to read in (the number of elements is zero),
; then the program sets !ERR to -1, and DATA is unmodified.
;
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; C. Markwardt, based in concept on FXBREAD version 12 from
; IDLASTRO, but with significant and
; major changes to accomodate the
; multiple row/column technique. Mostly
; the parameter checking and general data
; flow remain.
; C. Markwardt, updated to read variable length arrays, and to
; pass columns by handle or pointer.
; 20 Jun 2001
; C. Markwardt, try to conserve memory when creating the arrays
; 13 Oct 2001
; Handle case of GE 50 columns, C. Markwardt, 18 Apr 2002
; Handle case where TSCAL/TZERO changes type of column,
; C. Markwardt, 23 Feb 2003
; Fix bug in handling of FOUND and numeric columns,
; C. Markwardt 12 May 2003
;
;
; NAME:
; FXBWRITM
; PURPOSE:
; Write multiple columns/rows to a disk FITS binary table file.
; EXPLANATION :
; A call to FXBWRITM will write multiple rows and multiple
; columns to a binary table in a single procedure call. Up to
; fifty columns may be read in a single pass. The file should
; have already been opened with FXBOPEN (with write access) or
; FXBCREATE. FXBWRITM optimizes writing multiple columns by
; first writing a large chunk of data to the FITS file all at
; once. FXBWRITM cannot write variable-length arrays; use
; FXBWRITE instead.
;
; The number of columns is limited to 50 if data are passed by
; positional argument. However, this limitation can be overcome
; by passing either pointers or handles to FXBWRITM. The user
; should set the PASS_METHOD keyword to 'POINTER' or 'HANDLE' as
; appropriate, and an array of pointers or handles to the data
; in the POINTERS keyword. The user is responsible for freeing
; the pointers.
;
; CALLING SEQUENCE:
; FXBWRITM, UNIT, COL, D0, D1, D2, ..., [ ROW= ]
;
; INPUT PARAMETERS:
; UNIT = Logical unit number corresponding to the file containing the
; binary table.
; D0,..D49= An IDL data array to be written to the file, one for
; each column.
; COL = Column in the binary table to place data in, starting from
; column one.
; OPTIONAL INPUT KEYWORDS:
; ROW = Either row number in the binary table to write data to,
; starting from row one, or a two element array containing a
; range of row numbers to write. If not passed, then
; the entire column is written.
; NANVALUE= Value signalling data dropout. All points corresponding to
; this value are set to be IEEE NaN (not-a-number). Ignored
; unless DATA is of type float, double-precision or complex.
; PASS_METHOD = A scalar string indicating method of passing
; data to FXBWRITM. One of 'ARGUMENT' (indicating
; pass by positional argument), 'POINTER' (indicating
; passing an array of pointers by the POINTERS
; keyword), or 'HANDLE' (indicating passing an array
; of handles by the POINTERS keyword).
; Default: 'ARGUMENT'
; POINTERS = If PASS_METHOD is 'POINTER' then the user must pass
; an array of IDL pointers to this keyword, one for
; each column. If PASS_METHOD is 'HANDLE' then pass
; an array of IDL handles instead. Ultimately the
; user is responsible for deallocating pointers and
; handles.
; BUFFERSIZE = Data are transferred in chunks to conserve
; memory. This is the size in bytes of each chunk.
; If a value of zero is given, then all of the data
; are transferred in one pass. Default is 32768 (32
; kB).
; OPTIONAL OUTPUT KEYWORDS:
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXBWRITE, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
; WARNMSG = Messages which are considered to be non-fatal
; "warnings" are returned in this output string.
; STATUS = An output array containing the status for each
; read, 1 meaning success and 0 meaning failure.
;
; PROCEDURE CALLS:
; HOST_TO_IEEE, PRODUCT()
; EXAMPLE:
; Write a binary table 'sample.fits' giving 43 X,Y positions and a
; 21 x 21 PSF at each position:
;
; (1) First, create sample values
; x = findgen(43) & y = findgen(43)+1 & psf = randomn(seed,21,21,32)
;
; (2) Create primary header, write it to disk, and make extension header
; fxhmake,header,/initialize,/extend,/date
; fxwrite,'sample.fits',header
; fxbhmake,header,43,'TESTEXT','Test binary table extension'
;
; (3) Fill extension header with desired column names
; fxbaddcol,1,header,x[0],'X' ;Use first element in each array
; fxbaddcol,2,header,y[0],'Y' ;to determine column properties
; fxbaddcol,3,header,psf[*,*,0],'PSF'
;
; (4) Write extension header to FITS file
; fxbcreate,unit,'sample.fits',header
;
; (5) Use FXBWRITM to write all data to the extension in a single call
; fxbwritm,unit,['X','Y','PSF'], x, y, psf
; fxbfinish,unit ;Close the file
;
; COMMON BLOCKS:
; Uses common block FXBINTABLE--see "fxbintable.pro" for more
; information.
; RESTRICTIONS:
; The binary table file must have been opened with FXBCREATE or
; FXBOPEN (with write access).
;
; The data must be consistent with the column definition in the binary
; table header.
;
; The row number must be consistent with the number of rows stored in the
; binary table header.
;
; CATEGORY:
; Data Handling, I/O, FITS, Generic.
; PREVIOUS HISTORY:
; C. Markwardt, based on FXBWRITE and FXBREADM (ver 1), Jan 1999
; WRITTEN:
; Craig Markwardt, GSFC, January 1999.
; MODIFIED:
; Version 1, Craig Markwardt, GSFC 18 January 1999.
; Documented this routine, 18 January 1999.
; C. Markwardt, added ability to pass by handle or pointer.
; Some bug fixes, 20 July 2001
; W. Landsman/B.Schulz Allow more than 50 arguments when using pointers
; NAME:
; FXHMODIFY
; Purpose :
; Modify a FITS header in a file on disk.
; Explanation :
; Opens a FITS file, and adds or modifies a parameter in the FITS header.
; Can be used for either the main header, or for an extension header.
; The modification is performed directly on the disk file.
; Use :
; FXHMODIFY, FILENAME, NAME, VALUE, COMMENT
; Inputs :
; FILENAME = String containing the name of the file to be read.
;
; NAME = Name of parameter. If NAME is already in the header the
; value and possibly comment fields are modified. Otherwise a
; new record is added to the header. If NAME is equal to
; either "COMMENT" or "HISTORY" then the value will be added to
; the record without replacement. In this case the comment
; parameter is ignored.
;
; VALUE = Value for parameter. The value expression must be of the
; correct type, e.g. integer, floating or string. String
; values of 'T' or 'F' are considered logical values.
;
; Opt. Inputs :
; COMMENT = String field. The '/' is added by this routine. Added
; starting in position 31. If not supplied, or set equal to ''
; (the null string), then any previous comment field in the
; header for that keyword is retained (when found).
; Outputs :
; None.
; Opt. Outputs:
; None.
; Keywords :
; EXTENSION = Either the number of the FITS extension, starting with the
; first extension after the primary data unit being one; or a
; character string containing the value of EXTNAME to search
; for. If not passed, then the primary FITS header is
; modified.
;
; BEFORE = Keyword string name. The parameter will be placed before the
; location of this keyword. For example, if BEFORE='HISTORY'
; then the parameter will be placed before the first history
; location. This applies only when adding a new keyword;
; keywords already in the header are kept in the same position.
;
; AFTER = Same as BEFORE, but the parameter will be placed after the
; location of this keyword. This keyword takes precedence over
; BEFORE.
;
; FORMAT = Specifies FORTRAN-like format for parameter, e.g. "F7.3". A
; scalar string should be used. For complex numbers the format
; should be defined so that it can be applied separately to the
; real and imaginary parts.
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXHMODIFY, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; FXHREAD, FXPAR, FXADDPAR, BLKSHIFT
; Common :
; None.
; Restrictions:
; This routine can not be used to modify any of the keywords that control
; the structure of the FITS file, e.g. BITPIX, NAXIS, PCOUNT, etc. Doing
; so could corrupt the readability of the FITS file.
;
; Side effects:
; If adding a record to the FITS header would increase the
; number of 2880 byte records stored on disk, then the file is
; enlarged before modification, unless the NOGROW keyword is passed.
;
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; None.
; Written :
; William Thompson, GSFC, 3 March 1994.
; Modified :
; Version 1, William Thompson, GSFC, 3 March 1994.
; Version 2, William Thompson, GSFC, 31 May 1994
; Added ERRMSG keyword.
; Version 3, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version :
; Version 3, 23 June 1994
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; FXPAR()
; PURPOSE:
; Obtain the value of a parameter in a FITS header.
; EXPLANATION:
; The first 8 chacters of each element of HDR are searched for a match to
; NAME. If the keyword is one of those allowed to take multiple values
; ("HISTORY", "COMMENT", or " " (blank)), then the value is taken
; as the next 72 characters. Otherwise, it is assumed that the next
; character is "=", and the value (and optional comment) is then parsed
; from the last 71 characters. An error occurs if there is no parameter
; with the given name.
;
; If the value is too long for one line, it may be continued on to the
; the next input card, using the OGIP CONTINUE convention. For more info,
; http://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/ofwg_recomm/r13.html
;
; Complex numbers are recognized as two numbers separated by one or more
; space characters.
;
; If a numeric value has no decimal point (or E or D) it is returned as
; type LONG. If it contains more than 8 numerals, or contains the
; character 'D', then it is returned as type DOUBLE. Otherwise it is
; returned as type FLOAT. If an integer is too large to be stored as
; type LONG, then it is returned as DOUBLE.
;
; CALLING SEQUENCE:
; Result = FXPAR( HDR, NAME [, ABORT, COUNT=, COMMENT=, /NOCONTINUE ] )
;
; Result = FXPAR(HEADER,'DATE') ;Finds the value of DATE
; Result = FXPAR(HEADER,'NAXIS*') ;Returns array dimensions as
; ;vector
; REQUIRED INPUTS:
; HDR = FITS header string array (e.g. as returned by FXREAD). Each
; element should have a length of 80 characters
; NAME = String name of the parameter to return. If NAME is of the
; form 'keyword*' then an array is returned containing values
; of keywordN where N is an integer. The value of keywordN
; will be placed in RESULT(N-1). The data type of RESULT will
; be the type of the first valid match of keywordN found.
; OPTIONAL INPUT:
; ABORT = String specifying that FXPAR should do a RETALL if a
; parameter is not found. ABORT should contain a string to be
; printed if the keyword parameter is not found. If not
; supplied, FXPAR will return with a negative !err if a keyword
; is not found.
; START = A best-guess starting position of the sought-after
; keyword in the header. If specified, then FXPAR
; first searches for scalar keywords in the header in
; the index range bounded by START-PRECHECK and
; START+POSTCHECK. This can speed up keyword searches
; in large headers. If the keyword is not found, then
; FXPAR searches the entire header.
;
; If not specified then the entire header is searched.
; Searches of the form 'keyword*' also search the
; entire header and ignore START.
;
; Upon return START is changed to be the position of
; the newly found keyword. Thus the best way to
; search for a series of keywords is to search for
; them in the order they appear in the header like
; this:
;
; START = 0L
; P1 = FXPAR('P1', START=START)
; P2 = FXPAR('P2', START=START)
; PRECHECK = If START is specified, then PRECHECK is the number
; of keywords preceding START to be searched.
; Default: 5
; POSTCHECK = If START is specified, then POSTCHECK is the number
; of keywords after START to be searched.
; Default: 20
; OUTPUT:
; The returned value of the function is the value(s) associated with the
; requested keyword in the header array.
;
; If the parameter is complex, double precision, floating point, long or
; string, then the result is of that type. Apostrophes are stripped from
; strings. If the parameter is logical, 1 is returned for T, and 0 is
; returned for F.
;
; If NAME was of form 'keyword*' then a vector of values are returned.
;
; OPTIONAL INPUT KEYWORDS:
; /NOCONTINUE = If set, then continuation lines will not be read, even
; if present in the header
; OPTIONAL OUTPUT KEYWORD:
; COUNT = Optional keyword to return a value equal to the number of
; parameters found by FXPAR.
; COMMENTS= Array of comments associated with the returned values.
;
; PROCEDURE CALLS:
; GETTOK(), VALID_NUM
; SIDE EFFECTS:
;
; The system variable !err is set to -1 if parameter not found, 0 for a
; scalar value returned. If a vector is returned it is set to the number
; of keyword matches found.
;
; If a keyword occurs more than once in a header, a warning is given,
; and the first occurence is used. However, if the keyword is "HISTORY",
; "COMMENT", or " " (blank), then multiple values are returned.
;
; NOTES:
; The functions SXPAR() and FXPAR() are nearly identical, although
; FXPAR() has slightly more sophisticated parsing. There is no
; particular reason for having two nearly identical procedures, but
; both are too widely used to drop either one.
;
; REVISION HISTORY:
; Version 1, William Thompson, GSFC, 12 April 1993.
; Adapted from SXPAR
; Version 2, William Thompson, GSFC, 14 October 1994
; Modified to use VALID_NUM instead of STRNUMBER. Inserted
; additional call to VALID_NUM to trap cases where character
; strings did not contain quotation marks.
; Version 3, William Thompson, GSFC, 22 December 1994
; Fixed bug with blank keywords, following suggestion by Wayne
; Landsman.
; Version 4, Mons Morrison, LMSAL, 9-Jan-98
; Made non-trailing ' for string tag just be a warning (not
; a fatal error). It was needed because "sxaddpar" had an
; error which did not write tags properly for long strings
; (over 68 characters)
; Version 5, Wayne Landsman GSFC, 29 May 1998
; Fixed potential problem with overflow of LONG values
; Version 6, Craig Markwardt, GSFC, 28 Jan 1998,
; Added CONTINUE parsing
; Version 7, Craig Markwardt, GSFC, 18 Nov 1999,
; Added START, PRE/POSTCHECK keywords for better performance
; NAME:
; FXREAD
; Purpose :
; Read basic FITS files.
; Explanation :
; Read the primary array from a disk FITS file. Optionally allows the
; user to read in only a subarray and/or every Nth pixel.
; Use :
; FXREAD, FILENAME, DATA [, HEADER [, I1, I2 [, J1, J2 ]] [, STEP]]
; Inputs :
; FILENAME = String containing the name of the file to be read.
; Opt. Inputs :
; I1,I2 = Data range to read in the first dimension. If passed, then
; HEADER must also be passed. If not passed, or set to -1,-1,
; then the entire range is read.
; J1,J2 = Data range to read in the second dimension. If passed, then
; HEADER and I1,J2 must also be passed. If not passed, or set
; to -1,-1, then the entire range is read.
; STEP = Step size to use in reading the data. If passed, then
; HEADER must also be passed. Default value is 1. Ignored if
; less than 1.
; Outputs :
; DATA = Data array to be read from the file.
; Opt. Outputs:
; HEADER = String array containing the header for the FITS file.
; Keywords :
; NANVALUE = Value signalling data dropout. All points corresponding to
; IEEE NaN (not-a-number) are set to this value. Ignored
; unless DATA is of type float or double-precision.
; PROMPT = If set, then the optional parameters are prompted for at the
; keyboard.
; AVERAGE = If set, then the array size is reduced by averaging pixels
; together rather than by subselecting pixels. Ignored unless
; STEP is nontrivial. Note: this is much slower.
; YSTEP = If passed, then STEP is the step size in the 1st dimension,
; and YSTEP is the step size in the 2nd dimension. Otherwise,
; STEP applies to both directions.
; NOSCALE = If set, then the output data will not be scaled using the
; optional BSCALE and BZERO keywords in the FITS header.
; Default is to scale, if and only if BSCALE and BZERO are
; present and nontrivial.
; NOUPDATE = If set, then the optional BSCALE and BZERO keywords in the
; optional HEADER array will not be changed. The default is
; to reset these keywords to BSCALE=1, BZERO=0. Ignored if
; NOSCALE is set.
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXREAD, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; GET_DATE, IEEE_TO_HOST, FXADDPAR, FXHREAD, FXPAR, WHERENAN
; Common :
; None.
; Restrictions:
; Groups are not supported.
;
; The optional parameters I1, I2, and STEP only work with one or
; two-dimensional arrays. J1 and J2 only work with two-dimensional
; arrays.
;
; Use of the AVERAGE keyword is not compatible with arrays with missing
; pixels.
;
; Side effects:
; If the keywords BSCALE and BZERO are present in the FITS header, and
; have non-trivial values, then the returned array DATA is formed by the
; equation
;
; DATA = BSCALE*original + BZERO
;
; However, this behavior can overridden by using the /NOSCALE keyword.
;
; If the data is scaled, then the optional HEADER array is changed so
; that BSCALE=1 and BZERO=0. This is so that these scaling parameters
; are not applied to the data a second time by another routine. Also,
; history records are added storing the original values of these
; constants. Note that only the returned array is modified--the header
; in the FITS file itself is untouched.
;
; If the /NOUPDATE keyword is set, however, then the BSCALE and BZERO
; keywords are not changed. It is then the user's responsibility to
; ensure that these parameters are not reapplied to the data. In
; particular, these keywords should not be present in any header when
; writing another FITS file, unless the user wants their values to be
; applied when the file is read back in. Otherwise, FITS readers will
; read in the wrong values for the data array.
;
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; W. Thompson, May 1992, based in part on READFITS by W. Landsman, and
; STSUB by M. Greason and K. Venkatakrishna.
; W. Thompson, Jun 1992, added code to interpret BSCALE and BZERO
; records, and added NOSCALE and NOUPDATE
; keywords.
; W. Thompson, Aug 1992, changed to call FXHREAD, and to add history
; records for BZERO, BSCALE.
; Written :
; William Thompson, GSFC, May 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 17 November 1993.
; Corrected bug with AVERAGE keyword on non-IEEE compatible
; machines.
; Corrected bug with subsampling on VAX machines.
; Version 3, William Thompson, GSFC, 31 May 1994
; Added ERRMSG keyword.
; Version 4, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 5, Zarro (SAC/GSFC), 14 Feb 1997
; Added I/O error checking
; Version 6, 20-May-1998, David Schlegel/W. Thompson
; Allow a single pixel to be read in.
; Change the signal to read in the entire array to be -1
; Version :
; Version 6, 20-May-1998
; NAME:
; FXWRITE
; Purpose :
; Write a disk FITS file.
; Explanation :
; Creates a disk FITS file and writes a FITS primary header, and
; optionally a primary data array.
; Use :
; FXWRITE, FILENAME, HEADER [, DATA ]
; Inputs :
; FILENAME = String containing the name of the file to be written.
; HEADER = String array containing the header for the FITS file.
; Opt. Inputs :
; DATA = IDL data array to be written to the file. If not passed,
; then it is assumed that extensions will be added to the
; file.
; Outputs :
; None.
; Opt. Outputs:
; None.
; Keywords :
; NANVALUE = Value signalling data dropout. All points corresponding to
; this value are set to be IEEE NaN (not-a-number). Ignored
; unless DATA is of type float, double-precision or complex.
; NOUPDATE = If set, then the optional BSCALE and BZERO keywords in the
; HEADER array will not be changed. The default is to reset
; these keywords to BSCALE=1, BZERO=0.
; ERRMSG = If defined and passed, then any error messages will be
; returned to the user in this parameter rather than
; depending on the MESSAGE routine in IDL. If no errors are
; encountered, then a null string is returned. In order to
; use this feature, ERRMSG must be defined first, e.g.
;
; ERRMSG = ''
; FXWRITE, ERRMSG=ERRMSG, ...
; IF ERRMSG NE '' THEN ...
;
; Calls :
; CHECK_FITS, GET_DATE, HOST_TO_IEEE, FXADDPAR, FXPAR
; Common :
; None.
; Restrictions:
; If DATA is passed, then HEADER must be consistent with it. If no data
; array is being written to the file, then HEADER must also be consistent
; with that. The routine FXHMAKE can be used to create a FITS header.
;
; If found, then the optional keywords BSCALE and BZERO in the HEADER
; array is changed so that BSCALE=1 and BZERO=0. This is so that these
; scaling parameters are not applied to the data a second time by another
; routine. Also, history records are added storing the original values
; of these constants.
;
; If the /NOUPDATE keyword is set, however, then the BSCALE and BZERO
; keywords are not changed. The user should then be aware that FITS
; readers will apply these numbers to the data, even if the data is
; already converted to floating point form.
;
; Groups are not supported.
;
; Side effects:
; None.
; Category :
; Data Handling, I/O, FITS, Generic.
; Prev. Hist. :
; W. Thompson, Jan 1992, from WRITEFITS by J. Woffard and W. Landsman.
; Differences include:
;
; * Made DATA array optional, and HEADER array mandatory.
; * Changed order of HEADER and DATA parameters.
; * No attempt made to fix HEADER array.
;
; W. Thompson, May 1992, changed open statement to force 2880 byte fixed
; length records (VMS). The software here does not
; depend on this file configuration, but other
; FITS readers might.
; W. Thompson, Aug 1992, added code to reset BSCALE and BZERO records,
; and added the NOUPDATE keyword.
; Written :
; William Thompson, GSFC, January 1992.
; Modified :
; Version 1, William Thompson, GSFC, 12 April 1993.
; Incorporated into CDS library.
; Version 2, William Thompson, GSFC, 31 May 1994
; Added ERRMSG keyword.
; Version 3, William Thompson, GSFC, 23 June 1994
; Modified so that ERRMSG is not touched if not defined.
; Version 4, William Thompson, GSFC, 12 August 1999
; Catch error if unable to open file.
; Version 4.1 Wayne Landsman, GSFC, 02 May 2000
; Remove !ERR in call to CHECK_FITS, Use ARG_PRESENT()
; Version :
; Version 4.1, 2 May 2000
; NAME:
; GAL_UVW
; PURPOSE:
; Calculate the Galactic space velocity (U,V,W) of star
; EXPLANATION:
; Calculates the Galactic space velocity U, V, W of star given its
; (1) coordinates, (2) proper motion, (3) distance (or parallax), and
; (4) radial velocity.
; CALLING SEQUENCE:
; GAL_UVW, U, V, W, [/LSR, RA=, DEC=, PMRA= ,PMDEC=, VRAD= , DISTANCE=
; PLX= ]
; OUTPUT PARAMETERS:
; U - Velocity (km/s) positive toward the Galactic *anti*center
; V - Velocity (km/s) positive in the direction of Galactic rotation
; W - Velocity (km/s) positive toward the North Galactic Pole
; REQUIRED INPUT KEYWORDS:
; User must supply a position, proper motion,radial velocity and distance
; (or parallax). Either scalars or vectors can be supplied.
; (1) Position:
; RA - Right Ascension in *Degrees*
; Dec - Declination in *Degrees*
; (2) Proper Motion
; PMRA = Proper motion in RA in arc units (typically milli-arcseconds/yr)
; PMDEC = Proper motion in Declination (typically mas/yr)
; (3) Radial Velocity
; VRAD = radial velocity in km/s
; (4) Distance or Parallax
; DISTANCE - distance in parsecs
; or
; PLX - parallax with same distance units as proper motion measurements
; typically milliarcseconds (mas)
;
; OPTIONAL INPUT KEYWORD:
; /LSR - If this keyword is set, then the output velocities will be
; corrected for the solar motion (U,V,W)_Sun = (-9,+12,+7)
; (Mihalas & Binney, 1981) to the local standard of rest
; EXAMPLE:
; (1) Compute the U,V,W coordinates for the halo star HD 6755.
; Use values from Hipparcos catalog, and correct to the LSR
; ra = ten(1,9,42.3)*15. & dec = ten(61,32,49.5)
; pmra = 627.89 & pmdec = 77.84 ;mas/yr
; dis = 144 & vrad = -321.4
; gal_uvw,u,v,w,ra=ra,dec=dec,pmra=pmra,pmdec=pmdec,vrad=vrad,dis=dis,/lsr
; ===> u=154 v = -493 w = 97 ;km/s
;
; (2) Use the Hipparcos Input and Output Catalog IDL databases (see
; http://idlastro.gsfc.nasa.gov/ftp/zdbase/) to obtain space velocities
; for all stars within 10 pc with radial velocities > 10 km/s
;
; dbopen,'hipparcos,hic' ;Need Hipparcos output and input catalogs
; list = dbfind('plx>100,vrad>10') ;Plx > 100 mas, Vrad > 10 km/s
; dbext,list,'pmra,pmdec,vrad,ra,dec,plx',pmra,pmdec,vrad,ra,dec,plx
; ra = ra*15. ;Need right ascension in degrees
; GAL_UVW,u,v,w,ra=ra,dec=dec,pmra=pmra,pmdec=pmdec,vrad=vrad,plx = plx
; forprint,u,v,w ;Display results
; METHOD:
; Follows the general outline of Johnson & Soderblom (1987, AJ, 93,864)
; except that U is positive outward toward the Galactic *anti*center, and
; the J2000 transformation matrix to Galactic coordinates is taken from
; the introduction to the Hipparcos catalog.
; REVISION HISTORY:
; Written, W. Landsman December 2000
; NAME:
; GEO2GEODETIC
;
; PURPOSE:
; Convert from geographic/planetographic to geodetic coordinates
; EXPLANATION:
; Converts from geographic (latitude, longitude, altitude) to geodetic
; (latitude, longitude, altitude). In geographic coordinates, the
; Earth is assumed a perfect sphere with a radius equal to its equatorial
; radius. The geodetic (or ellipsoidal) coordinate system takes into
; account the Earth's oblateness.
;
; Geographic and geodetic longitudes are identical.
; Geodetic latitude is the angle between local zenith and the equatorial plane.
; Geographic and geodetic altitudes are both the closest distance between
; the satellite and the ground.
;
; The PLANET keyword allows a similar transformation for the other
; planets (planetographic to planetodetic coordinates).
;
; The EQUATORIAL_RADIUS and POLAR_RADIUS keywords allow the
; transformation for any ellipsoid.
;
; Latitudes and longitudes are expressed in degrees, altitudes in km.
;
; REF: Stephen P. Keeler and Yves Nievergelt, "Computing geodetic
; coordinates", SIAM Rev. Vol. 40, No. 2, pp. 300-309, June 1998
;
; Planterary constants from "Allen's Astrophysical Quantities",
; Fourth Ed., (2000)
;
; CALLING SEQUENCE:
; ecoord=geo2geodetic(gcoord,[ PLANET=,EQUATORIAL_RADIUS=, POLAR_RADIUS=])
;
; INPUT:
; gcoord = a 3-element array of geographic [latitude,longitude,altitude],
; or an array [3,n] of n such coordinates.
;
;
; OPTIONAL KEYWORD INPUT:
; PLANET = keyword specifying planet (default is Earth). The planet
; may be specified either as an integer (1-9) or as one of the
; (case-independent) strings 'mercury','venus','earth','mars',
; 'jupiter','saturn','uranus','neptune', or 'pluto'
;
; EQUATORIAL_RADIUS : Self-explanatory. In km. If not set, PLANET's
; value is used.
; POLAR_RADIUS : Self-explanatory. In km. If not set, PLANET's value is
; used.
;
; OUTPUT:
; a 3-element array of geodetic/planetodetic [latitude,longitude,altitude],
; or an array [3,n] of n such coordinates, double precision.
;
; COMMON BLOCKS:
; None
;
; RESTRICTIONS:
;
; Whereas the conversion from geodetic to geographic coordinates is given
; by an exact, analytical formula, the conversion from geographic to
; geodetic isn't. Approximative iterations (as used here) exist, but tend
; to become less good with increasing eccentricity and altitude.
; The formula used in this routine should give correct results within
; six digits for all spatial locations, for an ellipsoid (planet) with
; an eccentricity similar to or less than Earth's.
; More accurate results can be obtained via calculus, needing a
; non-determined amount of iterations.
; In any case,
; IDL> PRINT,geodetic2geo(geo2geodetic(gcoord)) - gcoord
; is a pretty good way to evaluate the accuracy of geo2geodetic.pro.
;
; EXAMPLES:
;
; Locate the geographic North pole, altitude 0., in geodetic coordinates
; IDL> geo=[90.d0,0.d0,0.d0]
; IDL> geod=geo2geodetic(geo); convert to equivalent geodetic coordinates
; IDL> PRINT,geod
; 90.000000 0.0000000 21.385000
;
; As above, but for the case of Mars
; IDL> geod=geo2geodetic(geo,PLANET='Mars')
; IDL> PRINT,geod
; 90.000000 0.0000000 18.235500
;
; MODIFICATION HISTORY:
; Written by Pascal Saint-Hilaire (shilaire@astro.phys.ethz.ch), May 2002
; Generalized for all solar system planets by Robert L. Marcialis
; (umpire@lpl.arizona.edu), May 2002
; Modified 2002/05/18, PSH: added keywords EQUATORIAL_RADIUS and
; POLAR_RADIUS
; NAME:
; GEODETIC2GEO
;
; PURPOSE:
; Convert from geodetic (or planetodetic) to geographic coordinates
; EXPLANATION:
; Converts from geodetic (latitude, longitude, altitude) to geographic
; (latitude, longitude, altitude). In geographic coordinates, the
; Earth is assumed a perfect sphere with a radius equal to its equatorial
; radius. The geodetic (or ellipsoidal) coordinate system takes into
; account the Earth's oblateness.
;
; Geographic and geodetic longitudes are identical.
; Geodetic latitude is the angle between local zenith and the equatorial
; plane. Geographic and geodetic altitudes are both the closest distance
; between the satellite and the ground.
;
; The PLANET keyword allows a similar transformation for the other
; planets (planetodetic to planetographic coordinates).
;
; The EQUATORIAL_RADIUS and POLAR_RADIUS keywords allow the
; transformation for any ellipsoid.
;
; Latitudes and longitudes are expressed in degrees, altitudes in km.
;
; REF: Stephen P. Keeler and Yves Nievergelt, "Computing geodetic
; coordinates", SIAM Rev. Vol. 40, No. 2, pp. 300-309, June 1998
; Planetary constants from "Allen's Astrophysical Quantities",
; Fourth Ed., (2000)
;
; CALLING SEQUENCE:
; gcoord = geodetic2geo(ecoord, [ PLANET= ] )
;
; INPUT:
; ecoord = a 3-element array of geodetic [latitude,longitude,altitude],
; or an array [3,n] of n such coordinates.
;
; OPTIONAL KEYWORD INPUT:
; PLANET = keyword specifying planet (default is Earth). The planet
; may be specified either as an integer (1-9) or as one of the
; (case-independent) strings 'mercury','venus','earth','mars',
; 'jupiter','saturn','uranus','neptune', or 'pluto'
;
; EQUATORIAL_RADIUS : Self-explanatory. In km. If not set, PLANET's value
; is used. Numeric scalar
; POLAR_RADIUS : Self-explanatory. In km. If not set, PLANET's value is
; used. Numeric scalar
;
; OUTPUT:
; a 3-element array of geographic [latitude,longitude,altitude], or an
; array [3,n] of n such coordinates, double precision
;
; The geographic and geodetic longitudes will be identical.
; COMMON BLOCKS:
; None
;
; EXAMPLES:
;
; IDL> geod=[90,0,0] ; North pole, altitude 0., in geodetic coordinates
; IDL> geo=geodetic2geo(geod)
; IDL> PRINT,geo
; 90.000000 0.0000000 -21.385000
;
; As above, but the equivalent planetographic coordinates for Mars
; IDL> geod=geodetic2geo(geod,PLANET='Mars');
; IDL> PRINT,geod
; 90.000000 0.0000000 -18.235500
;
; MODIFICATION HISTORY:
; Written by Pascal Saint-Hilaire (shilaire@astro.phys.ethz.ch),
; May 2002
;
; Generalized for all solar system planets by Robert L. Marcialis
; (umpire@lpl.arizona.edu), May 2002
;
; Modified 2002/05/18, PSH: added keywords EQUATORIAL_RADIUS and
; POLAR_RADIUS
;
; NAME:
; GLACTC
; PURPOSE:
; Convert between celestial and Galactic (or Supergalactic) coordinates.
; EXPLANATION:
; Program to convert right ascension (ra) and declination (dec) to
; Galactic longitude (gl) and latitude (gb) (j=1) or vice versa (j=2).
;
; CALLING SEQUENCE:
; GLACTC, ra, dec, year, gl, gb, j, [ /DEGREE, /FK4, /SuperGalactic ]
;
; INPUT PARAMETERS:
; year equinox of ra and dec, scalar (input)
; j direction of conversion (input)
; 1: ra,dec --> gl,gb
; 2: gl,gb --> ra,dec
;
; INPUTS OR OUTPUT PARAMETERS: ( depending on argument J )
; ra Right ascension, hours (or degrees if /DEGREES is set),
; scalar or vector
; dec Declination, degrees,scalar or vector
; gl Galactic longitude, degrees, scalar or vector
; gb Galactic latitude, degrees, scalar or vector
;
; All results forced double precision floating.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
; /DEGREE - If set, then the RA parameter (both input and output) is
; given in degrees rather than hours.
; /FK4 - If set, then the celestial (RA, Dec) coordinates are assumed
; to be input/output in the FK4 system. By default, coordinates
; are assumed to be in the FK5 system. For B1950 coordinates,
; set the /FK4 keyword *and* set the year to 1950.
; /SuperGalactic - If set, the GLACTC returns SuperGalactic coordinates
; as defined by deVaucouleurs et al. (1976) to account for the
; local supercluster. The North pole in SuperGalactic coordinates
; has Galactic coordinates l = 47.47, b = 6.32, and the origin is
; at Galactic coordinates l = 137.37, b= 0
;
; EXAMPLES:
; Find the Galactic coordinates of Altair (RA (J2000): 19 50 47
; Dec (J2000): 08 52 06)
;
; IDL> glactc, ten(19,50,47),ten(8,52,6),2000,gl,gb,1
; ==> gl = 47.74, gb = -8.91
;
; PROCEDURE CALLS:
; BPRECESS, JPRECESS, PRECESS
; HISTORY:
; FORTRAN subroutine by T. A. Nagy, 21-MAR-78.
; Conversion to IDL, R. S. Hill, STX, 19-OCT-87.
; Modified to handle vector input, E. P. Smith, GSFC, 14-OCT-94
; Converted to IDL V5.0 W. Landsman September 1997
; Added DEGREE keyword, C. Markwardt, Nov 1999
; Major rewrite, default now FK5 coordinates, added /FK4 keyword
; use external precession routines W. Landsman April 2002
; Add /Supergalactic keyword W. Landsman September 2002
; Fix major bug when year not 2000 and /FK4 not set W. Landsman July 2003
; NAME:
; HASTROM
; PURPOSE:
; Linear transformation of an image to align it with a reference image
; EXPLANATION:
; A linear transformation is applied (using POLY_2D) to an image so that
; its astrometry is identical with that in a reference header. This
; procedure can be used to align two images.
;
; CALLING SEQUENCE:
; HASTROM, oldim, oldhd, newim, newhd, refhd, [MISSING =, INTERP = ]
; or
; HASTROM, oldim, oldhd, refhd, [MISSING =, INTERP ={0,1,2}, NGRID =,
; CUBIC =, DEGREE = ]
;
; INPUTS:
; OLDIM - Image array to be manipulated. If only 3 parameters are
; supplied then OLDIM and OLDHD will be modified to contain
; the output image array and header
; OLDHD - FITS header array for OLDIM, containing astrometry parameters
; REFHD - Reference header, containing astrometry parameters. OLDIM
; will be rotated, shifted, and compressed or expanded until
; its astrometry matches that in REFHD.
; OUTPUTS:
; NEWIM - Image array after linear tranformation has been performed.
; The dimensions of NEWIM will be identical to the NAXIS1 and
; NAXIS2 keywords specified in REFHD. Regions on the reference
; image that do not exist in OLDIM can be assigned a value with
; the MISSING keyword.
; NEWHD - Updated FITS image header associated with NEWIM
;
; OPTIONAL INPUT KEYWORDS:
; MISSING - Set this keyword to a scalar value which will be assigned
; to pixels in the output image which are out of range of the
; supplied imput image. If not supplied, then linear
; extrapolation is used. See the IDL manual on POLY_2D.
; ***NOTE: A bug was introduced into the POLY_2D function in IDL
; V5.5 (still present in V5.6) such that the MISSING keyword
; may not work properly.***
; INTERP - Scalar, one of 0, 1, or 2 determining type of interpolation
; 0 nearest neighbor, 1 (default) bilinear interpolation,
; 2 cubic interpolation.
; CUBIC - a scalar value between -1 and 0 specifying cubic interpolation
; with the specified value as the cubic interpolation parameter.
; (see poly_2d for info). Setting CUBIC to a value greater
; than zero is equivalent to setting CUBIC = -1.
; NGRID - Integer scalar specifying the number of equally spaced grid
; points on each axis to use to specify the transformation.
; Default is NGRID = 3 (9 total grid points). The value of
; NGRID must always be greater than DEGREE + 1
; DEGREE - Integer scalar specifying the degree of the transformation.
; See the routine POLYWARP for more info. Default = 1
; (linear transformation).
; OPTIONAL OUTPUT KEYWORD:
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
; NOTES:
; (1) The 3 parameter calling sequence is less demanding on virtual
; memory.
; (2) The astrometry in OLDHD will be precessed to match the equinox
; given in REFHD.
; (3) If an ST Guidestar image is used for the reference header, then the
; output header will be converted to standard astrometry.
; EXAMPLE:
; Suppose one has an image array, IM, and an associated FITS header H.
; One desires to warp the image array so that it is aligned with another
; image with a FITS header, HREF. Both headers contain astrometry info.
; Set pixel values to 0 where there is no overlap between the input and
; reference image, and use linear interpolation (default)
;
; IDL> hastrom, IM, H, HREF, MISSING = 0
;
; PROCEDURES USED:
; ad2xy, check_FITS, extast, get_EQUINOX(), gsssextast, hprecess,
; putast, sxaddpar, sxaddhist, sxpar(), xy2ad, zparcheck
;
; REVISION HISTORY:
; Written W. Landsman, STX Co. Feb, 1989
; Updated to CHECK_FITS Dec, 1991
; New astrometry keywords Mar, 1994
; Recognize GSSS header W. Landsman June, 1994
; Added CUBIC keyword W. Landsman March, 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Accept INTERP=0, Convert output GSS header to standard astrometry
; W. Landsman June 1998
; Remove calls to obsolete !ERR system variable March 2000
; Added ERRMSG output keyword W. Landsman April 2000
; Need to re-extract astrometry after precession W. Landsman Nov. 2000
;
; NAME:
; HCONGRID
; PURPOSE:
; CONGRID an image and update astrometry in a FITS header
; EXPLANATION:
; Expand or contract an image using CONGRID and update the
; associated FITS header array.
;
; CALLING SEQUENCE:
; HCONGRID, oldhd ;Update FITS header only
; HCONGRID, oldim, oldhd, [ newim, newhd, newx, newy, /HALF_HALF
; CUBIC = , INTERP=, OUTSIZE =, ERRMSG = ]
;
; INPUTS:
; OLDIM - the original image array
; OLDHD - the original image FITS header, string array
;
; OPTIONAL INPUTS:
; NEWX - size of the new image in the X direction
; NEWY - size of the new image in the Y direction
; The OUTSIZE keyword can be used instead of the
; NEWX, NEWY parameters
;
; OPTIONAL OUTPUTS:
; NEWIM - the image after expansion or contraction with CONGRID
; NEWHD - header for newim containing updated astrometry info
; If output parameters are not supplied, the program
; will modify the input parameters OLDIM and OLDHD
; to contain the new array and updated header.
;
; OPTIONAL KEYWORD INPUTS:
; CUBIC - If set and non-zero, then cubic interpolation is used. Valid
; ranges are -1 <= Cubic < 0. Setting /CUBIC is equivalent to
; CUBIC = -1 and also equivalent to INTERP = 2. See INTERPOLATE
; for more info. Setting CUBIC = -0.5 is recommended.
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
; /HALF_HALF - Due to edge effects, the default behaviour of CONGRID is
; to introduce a slight shift in the image center. Craig Markwardt
; (http://cow.physics.wisc.edu/~craigm/idl/misc.html) has written
; a modified version of CONGRID called CMCONGRID that when used with
; the /HALF_HALF keyword eliminates any shift. The use of the
; /HALF keyword emulates CMCONGRID and eliminates any shift in the
; image centroid.
; INTERP - 0 for nearest neighbor, 1 for bilinear interpolation
; (default), 2 for cubic (=-1) interpolation.
; OUTSIZE - Two element integer vector which can be used instead of the
; NEWX and NEWY parameters to specify the output image dimensions
; OPTIONAL KEYWORD OUTPUT:
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
; PROCEDURE:
; Expansion or contraction is done using the CONGRID function, unless
; HALF_HALF is set.
;
; The parameters BSCALE, NAXIS1, NAXIS2, CRPIX1, and CRPIX2 and
; the CD (or CDELT) parameters are updated for the new header.
;
; NOTES:
; A FITS header can be supplied as the first parameter without having
; to supply an image array. The astrometry in the FITS header will be
; updated to be appropriate to the specified image size.
;
; If the FITS header contains astrometry from a ST Guide Star image,
; then the astrometry will be converted to an approximately equivalent
; tangent projection before applying CONGRID.
; EXAMPLE:
; Congrid an 512 x 512 image array IM and FITS header H to size 300 x 300
; using cubic interpolation. Use the HALF_HALF keyword to avoid
; a shift of the image centroid
;
; IDL> hcongrid, IM ,H, OUT = [300, 300], CUBIC = -0.5, /HALF
;
; The variables IM and H will be modified to the new image size.
;
; PROCEDURES CALLED:
; CHECK_FITS, CONGRID(), EXTAST, GSSS_STDAST, SXADDHIST,
; SXADDPAR, SXPAR(), ZPARCHECK
; MODIFICATION HISTORY:
; Written, Aug. 1986 W. Landsman, STI Corp.
; Added interp keywords, J. Isensee, July, 1990
; Add cubic interpolation W. Landsman HSTX January 1994
; Recognize a GSSS FITS header W. Landsman June 1994
; Fix case where header but not image supplied W. Landsman May 1995
; Remove call to SINCE_VERSION() W. Landsman March 1996
; Assume since IDL V3.5, add CUBIC keyword W. Landsman March 1997
; Update BSCALE even if no astrometry present W. Landsman May 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Added HALF_HALF keyword W. Landsman February 2000
; Added ERRMSG keyword, use double precision formatting W.L. April 2000
; Recognize PC00n00m astrometry format W. Landsman December 2001
; Now works when both /INTERP and /HALF are set W. Landsman January 2002
; NAME:
; HELIO
; PURPOSE:
; Compute (low-precision) heliocentric coordinates for the planets.
; EXPLANATION:
; The mean orbital elements for epoch J2000 are used. These are derived
; from a 250 yr least squares fit of the DE 200 planetary ephemeris to a
; Keplerian orbit where each element is allowed to vary linearly with
; time. For dates between 1800 and 2050, this solution fits the
; terrestrial planet orbits to ~25" or better, but achieves only ~600"
; for Saturn.
;
; Use PLANET_COORDS (which calls HELIO) to get celestial (RA, Dec)
; coordinates of the planets
; CALLING SEQUENCE:
; HELIO, JD, LIST, HRAD, HLONG, HLAT, [/RADIAN]
; INPUTS:
; JD = Julian date, double precision scalar or vector
; LIST = List of planets array. May be a single number.
; 1 = merc, 2 = venus, ... 9 = pluto.
;
; OUTPUTS:
; HRAD = array of Heliocentric radii (A.U).
; HLONG = array of Heliocentric (ecliptic) longitudes (degrees).
; HLAT = array of Heliocentric latitudes (degrees).
; These output parameters will be dimensioned Nplanet by Ndate,
; where Nplanet is the number of elements of list, and Ndate is
; the number of elements of JD.
;
; OPTIONAL INPUT KEYWORD:
; /RADIAN - If set, then the output longitude and latitude are given in
; radians.
; EXAMPLE:
; (1) Find the current heliocentric positions of all the planets
;
; IDL> GET_JULDATE, jd ;Get current Julian date
; IDL> HELIO,jd,indgen(9)+1,hrad,hlong,hlat ;Get radius, long, and lat
;
; (2) Find heliocentric position of Mars on August 23, 2000
; IDL> JDCNV, 2000,08,23,0,jd
; IDL> HELIO,JD,2,HRAD,HLONG,HLAT
; ===> hrad = 1.6407 AU hlong = 124.3197 hlat = 1.7853
; For comparison, the JPL ephemeris gives
; hrad = 1.6407 AU hlong = 124.2985 hlat = 1.7845
; (3) Find the heliocentric positions of Mars and Venus for every day in
; November 2000
; IDL> JDCNV, 2000, 11, 1, 0, jd ;Julian date of November 1, 2000
; IDL> helio, jd+indgen(30), [4,2], hrad,hlong,hlat ;Mars=4, Venus=2
; hrad, hlong, and hlat will be dimensioned [2,30]
; first column contains Mars data, second column Venus
; COMMON BLOCKS:
; None
; ROUTINES USED:
; CIRRANGE - force angle between 0 and 2*!PI
; NOTES:
; (1) The calling sequence for this procedure was changed in August 2000
; (2) This program is based on the two-body model and thus neglects
; interactions between the planets. This is why the worst results
; are for Saturn. Use the procedure JPLEPHINTERp for more accurate
; positions using the JPL ephemeris. Also see
; http://ssd.jpl.nasa.gov/cgi-bin/eph for a more accurate ephemeris
; generator online.
; (3) The coordinates are given for equinox 2000 and *not* the equinox
; of the supplied date(s)
; MODIFICATION HISTORY:
; R. Sterner. 20 Aug, 1986.
; Code cleaned up a bit W. Landsman December 1992
; Converted to IDL V5.0 W. Landsman September 1997
; Major rewrite, use modern orbital elements, vectorize, more accurate
; solution to Kepler's equation W. Landsman August 2000
; Wasn't working for planet vectors W. Landsman August 2000
; NAME:
; HELIO_RV
;
; PURPOSE:
; Return the heliocentric radial velocity of a spectroscopic binary
;
; EXPLANATION:
; This function will return the heliocentric radial velocity of a
; spectroscopic binary star at a given heliocentric Julian date (HJD)
; given its orbit.
;
; CALLING SEQUENCE:
;
; Result = HELIO_RV ( Reduced_HJD ,T ,Period ,Gamma [,e ,Omega ] )
;
; INPUT:
;
; Reduced_HJD - Reduced_HJD of observation
; T - Reduced_HJD of periastron passage (max. +ve velocity
; for circular orbits)
; Period - the period in days
; Gamma - systemic velocity
; K - velocity semi-amplitude in the same units as Gamma.
; e - eccentricity of the orbit, default is 0.
; Omega - longitude of periastron in degrees. Must be specified for
; eccentric orbits.
;
; OUTPUT:
;
; The predicted heliocentric radial velocity in the same units as Gamma
; for the date(s) specified by Reduced_HJD.
;
; RESTRICTIONS:
;
; To ensure consistency with the routines JULDATE and HELIO_JD, the
; reduced HJD must be used throughtout.
;
; EXAMPLES:
;
; Example 1
;
; What was the heliocentric radial velocity of the primary component of HU Tau
; at 1730 UT 25 Oct 1994?
;
; IDL> juldate ,[94,10,25,17,30],JD ;Get Geocentric julian date
; IDL> hjd = helio_jd(jd,ten(04,38,16)*15.,ten(20,41,05)) ; Convert to HJD
; IDL> print, helio_rv(hjd,46487.5303D,2.0563056D,-6.0,59.3)
; -63.661180
;
; NB. 1. The routines JULDATE and HELIO_JD return a reduced HJD (HJD - 2400000)
; and so T and P must be specified in the same fashion.
; 2. The user should be careful to use double precision format to specify
; T and P to sufficient precision where necessary.
;
; Example 2
;
; Plot two cycles of an eccentric orbit, e=0.6, omega=45 for both
; components of a binary star
;
; IDL> phi=findgen(100)/50.0 ; Generates 100 phase points
; IDL> plot, phi,helio_rv(phi,0,1,0,100,0.6,45),yrange=[-100,150]
; IDL> oplot, phi,helio_rv(phi,0,1,0,50,0.6,45+180)
;
; This illustrates both the use of arrays to perform multiple calculations
; and generating radial velocities for a given phase by setting T=0 and P=1.
; Note also that omega has been changed by 180 degrees for the orbit of the
; second component (the same 'trick' can be used for circular orbits).
;
;
; MODIFICATION HISTORY:
;
; Written by: Pierre Maxted CUOBS, October, 1994
;
; Circular orbits handled by setting e=0 and omega=0 to allow
; binary orbits to be handled using omega and omega+180.
; Pierre Maxted,Feb 95
; BUG - omega was altered by the routine - corrected Feb 95,Pierre Maxted
; Iteration for E changed to that given by Reidel , Feb 95,Pierre Maxted
; /SINGLE keyword removed. May 96,Pierre Maxted
;
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; HEXTRACT
; PURPOSE:
; Extract a subimage from an array and update astrometry in FITS header
; EXPLANATION:
; Extract a subimage from an array and create a new FITS header with
; updated astrometry for the subarray
; CALLING SEQUENCE:
; HEXTRACT, Oldim, Oldhd, [ Newim, Newhd, x0, x1, y0, y1, /SILENT ]
; or
; HEXTRACT, Oldim, Oldhd, [x0, x1, y0, y1, /SILENT, ERRMSG = ]
;
; INPUTS:
; Oldim - the original image array
; Oldhd - the original image header
;
; OPTIONAL INPUTS:
; x0, x1, y0, y1 - respectively, first and last X pixel, and first and
; last Y pixel to be extracted from the original image, integer scalars.
; If omitted, HEXTRACT will prompt for these parameters
;
; OPTIONAL OUTPUTS:
; Newim - the new subarray extracted from the original image
; Newhd - header for newim containing updated astrometry info
; If output parameters are not supplied or set equal to
; -1, then the HEXTRACT will modify the input parameters
; OLDIM and OLDHD to contain the subarray and updated header.
;
; OPTIONAL INPUT KEYWORD:
; /SILENT - If set and non-zero, then a message describing the extraction
; is not printed at the terminal. This message can also be
; suppressed by setting !QUIET.
; OPTIONAL KEYWORD OUTPUT:
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
;
; PROCEDURE:
; The FITS header parameters NAXIS1, NAXIS2, CRPIX1, and CRPIX2 are
; updated for the extracted image.
;
; EXAMPLE:
; Read an image from a FITS file 'IMAGE', extract a 512 x 512 subimage
; with the same origin, and write to a new FITS file 'IMAGENEW'
;
; IDL> im = READFITS( 'IMAGE', hdr ) ;Read FITS files into IDL arrays
; IDL> hextract, im, h, 0, 511, 0, 511 ;Extract 512 x 512 subimage
; IDL> writefits, 'IMAGENEW', im ,h ;Write subimage to a FITS file
;
; PROCEDURES CALLED
; CHECK_FITS, STRN(), SXPAR(), SXADDPAR, SXADDHIST
; MODIFICATION HISTORY:
; Written, Aug. 1986 W. Landsman, STX Corp.
; Use astrometry structure, W. Landsman Jan, 1994
; Minor fix if bad Y range supplied W. Landsman Feb, 1996
; Added /SILENT keyword W. Landsman March, 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Added ERRMSG keyword W. Landsman May 2000
; NAME:
; HOR2EQ
;
; PURPOSE:
; Converts local horizon coords (alt-az) of something to equatorial (ra-dec).
;
; EXPLANATION:
; This is a nice code to calculate equatorial (ra,dec) coordinates from
; horizon (alt,az) coords. It is typically accurate to about 1 arcsecond
; or better (I have checked the output against the publicly available XEPHEM
; software). It preforms precession, nutation, abberation, and refraction
; corrections. The perhaps best thing about it is that it can take arrays
; as inputs, in all variables and keywords EXCEPT Lat, lon, and Altitude
; (the code assumes these aren't changing), and uses vector arithmetic in
; every calculation except when calculating the precession matrices.
;
; CALLING SEQUENCE:
;
; HOR2EQ, alt, az, jd, ra, dec, [ha, LAT= , LON= , /WS, OBSNAME= , $
; /B1950 , PRECESS_= 0, NUTATE_= 0, REFRACT_= 0, $
; ABERRATION_= 0, ALTITUDE= , /VERBOSE, _EXTRA= ]
;
;
; INPUT VARIABLES
; alt : altitude (in degrees) [scalar or vector]
; az : azimuth angle (in degrees, measured EAST from NORTH, but see
; keyword WS below.) [scalar or vector]
; JD : Julian Date [scalar or vector]
; Note: if RA and DEC are arrays, then alt and az will also be arrays.
; If RA and DEC are arrays, JD may be a scalar OR an array of
; the same dimensionality.
;
; OPTIONAL INPUT KEYWORDS:
; lat : north geodetic latitude of location in degrees
; lon : EAST longitude of location in degrees
; (Specify west longitude with a negative sign.)
; /WS : Set this to get the azimuth measured westward from south
; (not East of North).
; obsname : Set this to a valid observatory name to be used by the
; astrolib OBSERVATORY procedure, which will return the latitude
; and longitude to be used by this program.
; /B1950 : Set this if your ra and dec are specified in B1950,
; FK4 coordinates (instead of J2000, FK5)
; precess_ : Set this to 1 to force precession [default], 0 for no
; precession.
; nutate_ : Set this to 1 to force nutation [default], 0 for no nutation.
; aberration_ : Set this to 1 to force aberration correction [default],
; 0 for no correction.
; refract_ : Set to 1 to force refraction correction [default], 0 for
; no correction.
; altitude: The altitude of the observing location, in meters. [default=0].
; /verbose: Set this for verbose output. The default is verbose=0.
; _extra: This is for setting TEMPERATURE or PRESSURE explicity, which are
; used by CO_REFRACT to calculate the refraction effect of the
; atmosphere. If you don't set these, the program will make an
; intelligent guess as to what they are (taking into account your
; altitude). See CO_REFRACT for more details.
;
; OUTPUT VARIABLES
; ra : Right Ascension of object (J2000) in degrees (FK5); scalar or
; vector.
; dec : Declination of object (J2000) in degrees (FK5), scalar or vector.
; ha : hour angle (in degrees) (optional)
;
; DEPENDENCIES:
; NUTATE, PRECESS, ADSTRING(), SUNPOS, OBSERVATORY (from the astrolib)
; CO_NUTATE, CO_ABERRATION, CO_REFRACT, HADEC2ALTAZ
;
; BASIC STEPS
; Precess Ra-Dec to current equinox.
; Nutation Correction to Ra-Dec
; Aberration correction to Ra-Dec
; Calculate Local Mean Sidereal Time
; Calculate Local Apparent Sidereal Time
; Calculate Hour Angle
; Do Spherical Trig to find Apparent Alt-Az
; Apply refraction correction to find observed Alt.
;
;CORRECTIONS I DO NOT MAKE:
; * Deflection of Light by the sun due to GR. (typically milliarcseconds,
; can be arseconds within one degree of the sun)
; * The Effect of Annual Parallax (typically < 1 arcsecond)
; * and more (see below)
;
; TO DO
; * Better Refraction Correction. Need to put in wavelength dependence,
; and integrate through the atmosphere.
; * Topocentric Parallax Correction (will take into account elevation of
; the observatory)
; * Proper Motion (but this will require crazy lookup tables or something).
; * Difference between UTC and UT1 in determining LAST -- is this important?
; * Effect of Annual Parallax (is this the same as topocentric Parallax?)
; * Polar Motion
; * Better connection to Julian Date Calculator.
;
; EXAMPLE:
;
; You are at Kitt Peak National Observatory, looking at a star at azimuth
; angle 264d 55m 06s and elevation 37d 54m 41s (in the visible). Today is
; Dec 25, 2041 and the local time is 10 PM precisely. What is the ra and dec
; (J2000) of the star you're looking at? The temperature here is about 0
; Celsius, and the pressure is 781 millibars. The Julian date for this
; time is 2466879.7083333
;
; IDL> hor2eq, ten(37,54,41), ten(264,55,06), 2466879.7083333d, ra, dec, $
; /verb, obs='kpno', pres=781.0, temp=273.0
;
; The program produces this output (because the VERBOSE keyword was set):
;
; Latitude = +31 57 48.0 Longitude = *** 36 0.0 ; longitude prints weirdly b/c of negative input to ADSTRING!!
; Julian Date = 2466879.708333
; Az, El = 17 39 40.4 +37 54 41.0 (Observer Coords)
; Az, El = 17 39 40.4 +37 53 39.6 (Apparent Coords)
; LMST = +03 53 54.1
; LAST = +03 53 53.6
; Hour Angle = +03 38 30.1 (hh:mm:ss)
; Ra, Dec: 00 15 23.5 +15 25 1.9 (Apparent Coords)
; Ra, Dec: 00 15 24.2 +15 25 0.1 (J2041.9841)
; Ra, Dec: 00 13 14.1 +15 11 0.3 (J2000)
;
; The star is therefore Algenib! Compare the derived Ra, Dec with what XEPHEM
; got:
; Ra, Dec: 00 13 14.2 +15 11 1.0 (J2000)
;
; AUTHOR:
; Chris O'Dell
; Univ. of Wisconsin-Madison
; Observational Cosmology Laboratory
; Email: odell@cmb.physics.wisc.edu
; NAME:
; HROT
; PURPOSE:
; Rotate an image and create new FITS header with updated astrometry.
; EXPLANATION:
; Cubic, bilinear or nearest neighbor interpolation can be used.
;
; CALLING SEQUENCE:
; HROT, oldim, oldhd, [ newim, newhd, angle, xc, yc, int,
; MISSING =, INTERP =, CUBIC = , /PIVOT]
; INPUTS:
; OLDIM - the original image array
; OLDHD - the original FITS image header, string array
;
; OPTIONAL INPUTS:
; NEWIM - If NEWIM is set to -1, then the old image and header will
; be updated
; ANGLE - Rotation angle, degrees clockwise
; XC - X Center of rotation (-1 for center of image)
; YC - Y Center of rotation (-1 for center of image)
; INT - 0 for nearest neighbor, 1 for bilinear interpolation
; 2 for cubic interpolation.
;
; OPTIONAL OUTPUTS:
; NEWIM - the rotated image, with the same dimensions as Oldim
; NEWHD - header for newim containing updated astrometry info
; If output parameters are not supplied, the program
; will modify the input parameters OLDIM and OLDHD
; to contain the rotated image and updated header.
;
; OPTIONAL INPUT KEYWORD:
; MISSING - Set this keyword to a scalar value which will be assigned
; to pixels in the output image which do not correspond to
; existing imput images (e.g if one rotates off-center).
; If not supplied then linear extrapolation is used.
;
; INTERP - scalar set to either 0 (nearest neighbor interpolation),
; 1 (bilinear interpolation), or 2 (cubic interpolation).
; The interpolation type can be specified by either the INTERP
; keyword or the int parameter
;
; CUBIC - If set and non-zero then cubic interpolation is used (see ROT),
; which is equivalent to setting INT = 2. In IDL V5.0 and later,
; this keyword can also be set to a value between -1 and 0.
;
; /PIVOT - Setting this keyword causes the image to pivot around the point
; XC, YC, so that this point maps into the same point in the
; output image. If this keyword is set to 0 or omitted, then the
; point XC, YC in the input image is mapped into the center of
; the output image.
;
; OPTIONAL OUTPUT KEYWORD:
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
; EXAMPLE:
; Rotate an image non-interactively 30 degrees clockwise. Use
; bilinear interpolation, and set missing values to 0.
;
; IDL> HROT, im_old, h_old, im_new, h_new, 30, -1, -1, 1, MIS = 0
;
; As above but update the input image and header and pivot about (100,120)
;
; IDL> HROT, im_old, h_old, -1, -1, 30, 100, 120, 1, MIS = 0, /PIVOT
; RESTRICTIONS:
; Unlike the ROT procedure, HROT cannot be used to magnify or
; or demagnify an image. Use HCONGRID or HREBIN instead.
;
; PROCEDURE:
; The image array is rotated using the ROT procedure.
; The CD (or CROTA) and CRPIX parameters, if present in the FITS header,
; are updated for the new rotation.
; History records are also added to the header
;
; PROCEDURES USED:
; CHECK_FITS, EXTAST, GETOPT(), GETROT, ROT(), STRN(), SXADDPAR
;
; MODIFICATION HISTORY:
; Written, Aug. 1986 W. Landsman, ST Systems Corp.
; Added MISSING keyword, W. Landsman March, 1991
; Added cubic interpolation, use astrometry structure Feb 1994
; Removed call to SINCE_VERSION() W. Landsman March 1996
; Assume at least V3.5, add CUBIC parameter W. Landsman March 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Fix for CROTA2 defined and CDELT1 NE CDELT2, W. Landsman November 1998
; Fix documentation to specify clockwise rotation W. Landsman Dec. 1999
; Added /PIVOT keyword W. Landsman January 2000
; Added ERRMSG, Use double precision formatting, W. Landsman April 2000
; Consistent conversion between CROTA and CD matrix W. Landsman Oct 2000
; Work for both CD001001 and CDELT defined W. Landsman March 2001
; Recognize PC matrix astrometry W. Landsman December 2001
; Update astrometry correctly when /PIVOT applied W. Landsman March 2002
; Update CROTA2 astrometry correctly, approximate GSSS W.L. June 2003
; Work with CD1_1, PC1_1 and CROTA keywords W. L. July 2003
; NAME:
; IMCONTOUR
; PURPOSE:
; Make a contour plot labeled with astronomical coordinates.
; EXPLANATION:
; The type of coordinate display is controlled by the keyword TYPE
; Set TYPE=0 (default) to measure distances from the center of the image
; (IMCONTOUR will decide whether the plotting units will be in
; arc seconds, arc minutes, or degrees depending on image size.)
; Set /TYPE for standard RA and Dec labeling
;
; By using the /NODATA keyword, IMCONTOUR can also be used to simply
; provide astronomical labeling of a previously displayed image.
; CALLING SEQUENCE
; IMCONTOUR, im, hdr,[ /TYPE, /PUTINFO, XDELTA = , YDELTA =, _EXTRA =
; XMID=, YMID= ]
;
; INPUTS:
; IM - 2-dimensional image array
; HDR - FITS header associated with IM, string array, must include
; astrometry keywords. IMCONTOUR will also look for the
; OBJECT and IMAGE keywords, and print these if found and the
; PUTINFO keyword is set.
;
; OPTIONAL PLOTTING KEYWORDS:
; /TYPE - the type of astronomical labeling to be displayed. Either set
; TYPE = 0 (default), distance to center of the image is
; marked in units of Arc seconds, arc minutes, or degrees
;
; TYPE = 1 astronomical labeling with Right ascension and
; declination.
;
; /PUTINFO - If set, then IMCONTOUR will add information about the image
; to the right of the contour plot. Information includes image
; name, object, image center, image center, contour levels, and
; date plot was made
;
; XDELTA, YDELTA - Integer scalars giving spacing of labels for TYPE=1.
; Default is to label every major tick (XDELTA=1) but if
; crowding occurs, then the user might wish to label every other
; tick (XDELTA=2) or every third tick (XDELTA=3)
;
; XMID, YMID - Scalars giving the X,Y position from which offset distances
; will be measured when TYPE=0. By default, offset distances
; are measured from the center of the image.
;
; Any keyword accepted by CONTOUR may also be passed through IMCONTOUR
; since IMCONTOUR uses the _EXTRA facility. IMCONTOUR uses its own
; defaults for the XTITLE, YTITLE XMINOR, YMINOR, and SUBTITLE keywords
; but these may be overridden.
;
; NOTES:
; (1) The contour plot will have the same dimensional ratio as the input
; image array
; (2) To contour a subimage, use HEXTRACT before calling IMCONTOUR
; (3) Use the /NODATA keyword to simply provide astronomical labeling
; of a previously displayed image.
;
; EXAMPLE:
; Overlay the contour of an image, im2, with FITS header, h2, on top
; of the display of a different image, im1. Use RA, Dec labeling, and
; seven equally spaced contour levels. The use of a program like
; David Fanning's TVIMAGE http://www.dfanning.com/programs/tvimage.pro
; is suggested to properly overlay plotting and image coordinates. The
; /Keep_aspect_ratio keyword must be used.
;
; IDL> tvimage,im1,/keep_aspect, position = pos
; IDL> imcontour,im2,h2,nlevels=7,/Noerase,/TYPE,position = pos
;
; PROCEDURES USED:
; CHECK_FITS, EXTAST, GETROT, TICPOS, TICLABEL, TIC_ONE, TICS, XYAD
; CONS_RA(), CONS_DEC(), ADSTRING()
;
; REVISION HISTORY:
; Written W. Landsman STX May, 1989
; Fixed RA,Dec labeling W. Landsman November, 1991
; Fix plottting keywords W.Landsman July, 1992
; Recognize GSSS headers W. Landsman July, 1994
; Removed Channel keyword for V4.0 compatibility June, 1995
; Add _EXTRA CONTOUR plotting keywords W. Landsman August, 1995
; Add XDELTA, YDELTA keywords W. Landsman November, 1995
; Use SYSTIME() instead of !STIME August, 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Remove obsolete !ERR system variable W. Landsman May 2000
; Added XMID, YMID keywords to specify central position (default is still
; center of image) W. Landsman March 2002
; Recognize Galactic coordinates, fix Levels display when /PUTINFO set
; W. Landsman May 2003
; Correct conversion from seconds of RA to arcmin is 4 not 15.
; M. Perrin July 2003
;
; NAME:
; JPLEPHINTERP
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
; UPDATED VERSIONs can be found on my WEB PAGE:
; http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
; Interpolate position and motion of planetary bodies (JPL Ephemeris)
;
; MAJOR TOPICS:
; Planetary Orbits, Interpolation
;
; CALLING SEQUENCE:
; JPLEPHINTERP, INFO, RAWDATA, T, X, Y, Z, [VX, VY, VZ, /EARTH, /SUN,
; OBJECTNAME=, CENTER=, TBASE=, POSUNITS=, VELUNITS= ]
;
; DESCRIPTION:
;
; JPLEPHINTERP interpolates the JPL DE200 or DE405 planetary
; ephemeris to find the positions and motions of planetary bodies.
;
; This routine is the second stage of a two-stage process to
; interpolate the JPL ephemeris. In this first stage, the file is
; opened using JPLEPHREAD, and the relevant portions of the table
; are read and stored into the two variables INFO and RAWDATA. In
; the second stage, the user actually interpolates the ephemeris for
; the desired bodies and to the desired ephemeris time using
; JPLEPHINTERP.
;
; The only independent variable which must be specified is T, the
; ephemeris time. For low to moderate accuracy applications, T is
; simply the conventional calendar date, expressed in Julian days.
; See below for high precision applications.
;
; Upon output, the position components of the desired body are
; returned in parameters X, Y and Z, and if requested velocity
; components are returned in parameters VX, VY and VZ. Coordinates
; are referred to the ephemeris's coordinate system: FK5 for
; JPL-DE200 and ICRS for JPL-DE405. By default, the origin of
; coordinates is the solar system barycenter (SSB), unless another
; origin is selected using the CENTER keyword.
;
; Users must set the VELOCITY keyword to generate body velocities.
; By default they are not generated.
;
; Users can select the desired body by using either the EARTH or SUN
; keywords, or the OBJECTNAME keyword.
;
; By default, positions are returned in units of KM and velocities
; in units of KM/DAY. However, the output units are selectable by
; setting the POSUNITS and VELUNITS keywords.
;
; High Precision Applications
;
; If the required precision is finer than a few hundred meters, the
; user must be aware that the formal definition of the ephemeris
; time is the coordinate time of a clock placed at the solar system
; barycenter (SSB). If the user's time is measured by a clock
; positioned elsewhere, then various corrections must be applied.
; Usually, the most significant correction is that from the
; geocenter to the SSB (see Fairhead & Bretagnon 1990; Fukushima
; 1995). Not applying this correction creates an error with
; amplitude ~170 nano-light-seconds ( = 50 m) on the earth's
; position. (see TDB2TDT)
;
; For high precision, the user should also specify the TBASE
; keyword. TBASE should be considered a fixed epoch with respect to
; which T is measured; T should be small compared to TBASE.
; Internally, subtraction of large numbers occurs with TBASE first,
; so truncation error is minimized by specifying TBASE.
;
; Nutations and Librations
;
; This routine also provides information about earth nutations and
; lunar librations, which are stored in the JPL ephemeris tables.
; The POSUNITS and VELUNITS keywords do not affect these
; computations.
;
; Lunar librations in the form of three Euler angles are returned in
; X, Y, Z, in units of radians, and their time derivatives are
; returned in VX, VY, and VZ in units of radians per day.
;
; The earth nutation angles psi (nutation in longitude) and epsilon
; (nutation in obliquity) are returned in X and Y, in units of
; radians. Their time derivatives are returned in VX and VY
; respectively. The quantities returned in Z and VZ are undefined.
;
; Verification
;
; The precision routine has been verified using JPLEPHTEST, which is
; similar to the original JPL program EPHTEST. For years 1950 to
; 2050, JPLEPHINTERP reproduces the original JPL ephemeris to within
; 1.5 centimeters.
;
;
; PARAMETERS:
;
; INFO - structure returned by JPLEPHREAD. Users should not modify
; this structure.
;
; RAWDATA - raw data array returned by JPLEPHREAD. Users should not
; modify this data array.
;
; T - ephemeris time(s) of interest. May be a scalar or vector.
; The actual time is (T+TBASE).
;
; X, Y, Z - upon return, the x-, y- and z-components of the body
; position are returned in these parameters. For
; nutations and librations see above.
;
; VX, VY, VZ - upon return, the x-, y- and z-components of the body
; velocity are returned in these parameters, if the
; VELOCITY keyword is set. For nutations and
; librations see above.
;
;
; KEYWORD PARAMETERS:
;
; EARTH, SUN - set one of these keywords if the desired body is the
; earth or the sun. One of EARTH, SUN or OBJECTNAME
; must be specified.
;
; OBJECTNAME - a scalar string or integer, specifies the planetary
; body of interest. May take any one of the following
; integer or string values.
;
; 1 - 'MERCURY' 9 - 'PLUTO'
; 2 - 'VENUS' 10 - 'MOON' (earth's moon)
; 3 - 'EARTH' 11 - 'SUN'
; 4 - 'MARS' 12 - 'SOLARBARY' (solar system barycenter)
; 5 - 'JUPITER' 13 - 'EARTHBARY' (earth-moon barycenter)
; 6 - 'SATURN' 14 - 'NUTATIONS' (see above)
; 7 - 'URANUS' 15 - 'LIBRATIONS' (see above)
; 8 - 'NEPTUNE'
;
; CENTER - a scalar string or integer, specifies the origin of
; coordinates. See OBJECTNAME for allowed values.
; Default: 12 (Solar system barycenter)
;
; VELOCITY - if set, body velocities are generated and returned in
; VX, VY and VZ.
; Default: unset (no velocities)
;
; POSUNITS - a scalar string specifying the desired units for X, Y,
; and Z. Allowed values:
; 'KM' - kilometers (default)
; 'CM' - centimeters
; 'AU' - astronomical units
; 'LT-S' - light seconds
;
; VELUNITS - a scalar string specifying the desired units for VX, VY
; and VZ. Allowed values:
; 'KM/DAY' - kilometers per day (default)
; 'KM/S' - kilometers per second
; 'CM/S' - centimeters per second
; 'LT-S/S' - light seconds per second
; 'AU/DAY' - astronomical units per day
;
; TBASE - a scalar number, specifies a fixed epoch against wich T is
; measured. The ephemeris time will be (T+TBASE). Use this
; keyword for maximum precision.
;
;
; EXAMPLE:
;
; Find position of earth at ephemeris time 2451544.5 JD. Units are
; in Astronomical Units.
;
; JPLEPHREAD, 'JPLEPH.200', pinfo, pdata, [2451544D, 2451545D]
;
; JPLEPHINTERP, pinfo, pdata, 2451544.5D, xearth, yearth, zearth, $
; /EARTH, posunits='AU'
;
;
; REFERENCES:
;
; AXBARY, Arnold Rots.
; ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/bary/
;
; HORIZONS, JPL Web-based ephermis calculator (Ephemeris DE406)
; http://ssd.jpl.nasa.gov/horizons.html
;
; Fairhead, L. & Bretagnon, P. 1990, A&A, 229, 240
;
; Fukushima, T. 1995, A&A, 294, 895
;
; Standish, E.M. 1982, "Orientation of the JPL Ephemerides,
; DE200/LE200, to the Dynamical Equinox of J2000", Astronomy &
; Astrophysics, vol. 114, pp. 297-302.
;
; Standish, E.M.: 1990, "The Observational Basis for JPL's DE200,
; the planetary ephemeris of the Astronomical Almanac", Astronomy
; & Astrophysics, vol. 233, pp. 252-271.
;
; SEE ALSO
; JPLEPHREAD, JPLEPHINTERP, JPLEPHTEST, TDB2TDT
;
; MODIFICATION HISTORY:
; Written and Documented, CM, Jun 2001
; Corrected bug in name conversion of NUTATIONS and LIBRATIONS, 18
; Oct 2001, CM
;
; $Id: jplephinterp.pro,v 1.9 2001/10/18 19:17:57 craigm Exp $
;
; NAME:
; JPLEPHREAD
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
; UPDATED VERSIONs can be found on my WEB PAGE:
; http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
; Open and read JPL DE200 or DE405 Ephemeride FITS File
;
; MAJOR TOPICS:
; Planetary Orbits, Interpolation
;
; CALLING SEQUENCE:
; JPLEPHREAD, FILENAME, INFO, RAWDATA, JDLIMITS, STATUS=, ERRMSG=
;
; DESCRIPTION:
;
; JPLEPHREAD opens and reads the JPL DE200 or DE405 planetary
; ephemerides, as available in FITS format. The user must have the
; IDL Astronomy Library installed to use this routine.
;
; This routine is the initialization stage of a two-stage process to
; interpolate the JPL ephemeris. In this first stage, the file is
; opened, and the relevant portions of the table are read and stored
; into the two variables INFO and RAWDATA. In the second stage, the
; user actually interpolates the ephemeris for the desired bodies
; and to the desired ephemeris time using JPLEPHINTERP.
;
; Users must decide ahead of time the approximate dates of interest,
; and pass this range in the JDLIMITS parameter. Any date covered
; by the ephemeris is valid.
;
; JPLEPHREAD is able to read files of the following format:
; DE200 - Chebyshev - FITS format - Note 1
; DE405 - Chebyshev - FITS format - Note 1
; DE200 - Taylor - FITS format - Note 2
;
; Note 1 - Chebyshev formatted FITS files are available in the
; AXBARY package by Arnold Rots, found here:
; ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/bary/
; or at the Markwardt FTP site:
; ftp://cow.physics.wisc.edu/pub/craigm/bary/
;
; Note 2 - Taylor-series based ephemerides have been available for
; years in the FTOOLS / LHEASOFT package produced by NASA's
; Goddard Space Flight Center. The original file is
; de200_new.fits, which covers the years 1959-2000,
; inclusive. A newer file is named
; de200_1950-2050_v2.fits, and covers the years 1959-2050.
; See Markwardt FTP site for these files.
;
; PARAMETERS:
;
; FILENAME - name of ephemeris file (scalar string).
;
; INFO - upon completion, information about the ephemeris data is
; returned in this parameter in the form of a structure.
; Users must not modify INFO, although several fields are
; useful and may be accessed read-only:
; TSTART/TSTOP (start and stop time of data in Julian
; days);
; C (speed of light in km/s);
; DENUM (development ephemeris number [200 or 405])
; AU (1 astronomical unit, in units of light-seconds)
;
; RAWDATA - upon completion, raw ephemeris data is returned in this
; parameter. Users are not meant to access this data
; directly, but rather to pass it to JPLEPHINTERP.
;
; JDLIMITS - a two-element vector (optional), describing the desired
; time range of interest. The vector should have the
; form [TSTART, TSTOP], where TSTART and TSTOP are the
; beginning and ending times of the range, expressed in
; Julian days.
; Default: entire table is read (note, this can be
; several megabytes)
;
;
; KEYWORD PARAMETERS:
;
; STATUS - upon completion, a value of 1 indicates success, and 0
; indicates failure.
;
; ERRMSG - upon completion, an error message is returned in this
; keyword. If there were no errors, then the returned
; value is the empty string, ''.
;
;
; EXAMPLE:
;
; Find position of earth at ephemeris time 2451544.5 JD. Units are
; in Astronomical Units.
;
; JPLEPHREAD, 'JPLEPH.200', pinfo, pdata, [2451544D, 2451545D]
;
; JPLEPHINTERP, pinfo, pdata, 2451544.5D, xearth, yearth, zearth, $
; /EARTH, posunits='AU'
;
;
; REFERENCES:
;
; AXBARY, Arnold Rots.
; ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/bary/
;
; HORIZONS, JPL Web-based ephermis calculator (Ephemeris DE406)
; http://ssd.jpl.nasa.gov/horizons.html
;
; JPL Export Ephmeris FTP Site
; ftp://navigator.jpl.nasa.gov/pub/ephem/export/
; (ephemeris files are available here, however, they must be
; converted to FITS format using the "bin2eph" utility found in
; AXBARY)
;
; JPL Export Ephemeris CD-ROM - Ordering Information
; http://www.willbell.com/software/jpl.htm
;
; Standish, E.M. 1982, "Orientation of the JPL Ephemerides,
; DE200/LE200, to the Dynamical Equinox of J2000", Astronomy &
; Astrophysics, vol. 114, pp. 297-302.
;
; Standish, E.M.: 1990, "The Observational Basis for JPL's DE200,
; the planetary ephemeris of the Astronomical Almanac", Astronomy
; & Astrophysics, vol. 233, pp. 252-271.
;
; SEE ALSO
; JPLEPHREAD, JPLEPHINTERP, JPLEPHTEST
; PROCEDURES USED:
; FXBCLOSE, FXBOPEN, FXPAR(),
;
; MODIFICATION HISTORY:
; Written and Documented, CM, Jun 2001
; Use GETTOK() instead of STR_SEP() W. Landsman July 2002
;
; $Id: jplephread.pro,v 1.6 2001/07/01 03:32:02 craigm Exp $
;
; NAME:
; JPRECESS
; PURPOSE:
; Precess astronomical coordinates from B1950 to J2000
; EXPLANATION:
; Calculate the mean place of a star at J2000.0 on the FK5 system from the
; mean place at B1950.0 on the FK4 system.
;
; Use BPRECESS for the reverse direction J2000 ==> B1950
; CALLING SEQUENCE:
; jprecess, ra, dec, ra_2000, dec_2000, [ MU_RADEC = , PARALLAX =
; RAD_VEL =, EPOCH = ]
;
; INPUTS:
; RA,DEC - input B1950 right ascension and declination in *degrees*.
; Scalar or vector
;
; OUTPUTS:
; RA_2000, DEC_2000 - the corresponding J2000 right ascension and
; declination in *degrees*. Same number of elements as RA,DEC
; but always double precision.
;
; OPTIONAL INPUT-OUTPUT KEYWORDS
; MU_RADEC - 2xN element double precision vector containing the proper
; motion in seconds of arc per tropical *century* in right
; ascension and declination.
; PARALLAX - N_element vector giving stellar parallax (seconds of arc)
; RAD_VEL - N_element vector giving radial velocity in km/s
;
; The values of MU_RADEC, PARALLAX, and RADVEL will all be modified
; upon output to contain the values of these quantities in the
; J2000 system. Values will also be converted to double precision.
; The parallax and radial velocity will have a very minor influence on
; the J2000 position.
;
; EPOCH - scalar giving epoch of original observations, default 1950.0d
; This keyword value is only used if the MU_RADEC keyword is not set.
; NOTES:
; The algorithm is taken from the Explanatory Supplement to the
; Astronomical Almanac 1992, page 184.
; Also see Aoki et al (1983), A&A, 128,263
;
; JPRECESS distinguishes between the following two cases:
; (1) The proper motion is known and non-zero
; (2) the proper motion is unknown or known to be exactly zero (i.e.
; extragalactic radio sources). In this case, the algorithm
; in Appendix 2 of Aoki et al. (1983) is used to ensure that
; the output proper motion is exactly zero. Better precision
; can be achieved in this case by inputting the EPOCH of the
; original observations.
;
; The error in using the IDL procedure PRECESS for converting between
; B1950 and J2000 can be up to 1.5", mainly in right ascension. If
; better accuracy than this is needed then JPRECESS should be used.
;
; EXAMPLE:
; The SAO catalogue gives the B1950 position and proper motion for the
; star HD 119288. Find the J2000 position.
;
; RA(1950) = 13h 39m 44.526s Dec(1950) = 8d 38' 28.63''
; Mu(RA) = -.0259 s/yr Mu(Dec) = -.093 ''/yr
;
; IDL> mu_radec = 100D* [ -15D*.0259, -0.093 ]
; IDL> ra = ten(13,39,44.526)*15.D
; IDL> dec = ten(8,38,28.63)
; IDL> jprecess, ra, dec, ra2000, dec2000, mu_radec = mu_radec
; IDL> print, adstring(ra2000, dec2000,2)
; ===> 13h 42m 12.740s +08d 23' 17.69"
;
; RESTRICTIONS:
; "When transferring individual observations, as opposed to catalog mean
; place, the safest method is to tranform the observations back to the
; epoch of the observation, on the FK4 system (or in the system that was
; used to to produce the observed mean place), convert to the FK5 system,
; and transform to the the epoch and equinox of J2000.0" -- from the
; Explanatory Supplement (1992), p. 180
;
; REVISION HISTORY:
; Written, W. Landsman September, 1992
; Corrected a couple of typos in M matrix October, 1992
; Vectorized, W. Landsman February, 1994
; Implement Appendix 2 of Aoki et al. (1983) for case where proper
; motion unknown or exactly zero W. Landsman November, 1994
; Converted to IDL V5.0 W. Landsman September 1997
; Fixed typo in updating proper motion W. Landsman April 1999
; Make sure proper motion is floating point W. Landsman December 2000
; NAME:
; LEGEND
; PURPOSE:
; Create an annotation legend for a plot.
; EXPLANATION:
; This procedure makes a legend for a plot. The legend can contain
; a mixture of symbols, linestyles, Hershey characters (vectorfont),
; and filled polygons (usersym). A test procedure, legendtest.pro,
; shows legend's capabilities. Placement of the legend is controlled
; with keywords like /right, /top, and /center or by using a position
; keyword for exact placement (position=[x,y]) or via mouse (/position).
; CALLING SEQUENCE:
; LEGEND [,items][,keyword options]
; EXAMPLES:
; The call:
; legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
; produces:
; -----------------
; | |
; | + Plus sign |
; | * Asterisk |
; | . Period |
; | |
; -----------------
; Each symbol is drawn with a plots command, so they look OK.
; Other examples are given in optional output keywords.
;
; lines = indgen(6) ; for line styles
; items = 'linestyle '+strtrim(lines,2) ; annotations
; legend,items,linestyle=lines ; vertical legend---upper left
; items = ['Plus sign','Asterisk','Period']
; sym = [1,2,3]
; legend,items,psym=sym ; ditto except using symbols
; legend,items,psym=sym,/horizontal ; horizontal format
; legend,items,psym=sym,box=0 ; sans border
; legend,items,psym=sym,delimiter='=' ; embed '=' betw psym & text
; legend,items,psym=sym,margin=2 ; 2-character margin
; legend,items,psym=sym,position=[x,y] ; upper left in data coords
; legend,items,psym=sym,pos=[x,y],/norm ; upper left in normal coords
; legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
; legend,items,psym=sym,/position ; interactive position
; legend,items,psym=sym,/right ; at upper right
; legend,items,psym=sym,/bottom ; at lower left
; legend,items,psym=sym,/center ; approximately near center
; legend,items,psym=sym,number=2 ; plot two symbols, not one
; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
; INPUTS:
; items = text for the items in the legend, a string array.
; For example, items = ['diamond','asterisk','square'].
; You can omit items if you don't want any text labels.
; OPTIONAL INPUT KEYWORDS:
;
; linestyle = array of linestyle numbers If linestyle[i] < 0, then omit
; ith symbol or line to allow a multi-line entry. If
; linestyle = -99 then text will be left-justified.
; psym = array of plot symbol numbers. If psym[i] is negative, then a
; line connects pts for ith item. If psym[i] = 8, then the
; procedure usersym is called with vertices define in the
; keyword usersym. If psym[i] = 88, then use the previously
; defined user symbol
; vectorfont = vector-drawn characters for the sym/line column, e.g.,
; ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
; and a partial derivative, which might have accompanying items
; ['BOX','CHECK','PARTIAL DERIVATIVE'].
; There is no check that !p.font is set properly, e.g., -1 for
; X and 0 for PostScript. This can produce an error, e.g., use
; !20 with PostScript and !p.font=0, but allows use of Hershey
; *AND* PostScript fonts together.
; N. B.: Choose any of linestyle, psym, and/or vectorfont. If none is
; present, only the text is output. If more than one
; is present, all need the same number of elements, and normal
; plot behaviour occurs.
; By default, if psym is positive, you get one point so there is
; no connecting line. If vectorfont[i] = '',
; then plots is called to make a symbol or a line, but if
; vectorfont[i] is a non-null string, then xyouts is called.
; /help = flag to print header
; /horizontal = flag to make the legend horizontal
; /vertical = flag to make the legend vertical (D=vertical)
; box = flag to include/omit box around the legend (D=include)
; clear = flag to clear the box area before drawing the legend
; delimiter = embedded character(s) between symbol and text (D=none)
; colors = array of colors for plot symbols/lines (D=!P.color)
; textcolors = array of colors for text (D=!P.color)
; margin = margin around text measured in characters and lines
; spacing = line spacing (D=bit more than character height)
; pspacing = psym spacing (D=3 characters) (when number of symbols is
; greater than 1)
; charsize = just like !p.charsize for plot labels
; charthick = just like !p.charthick for plot labels
; thick = array of line thickness numbers (D = !P.thick), if used, then
; linestyle must also be specified
; position = data coordinates of the /top (D) /left (D) of the legend
; normal = use normal coordinates for position, not data
; device = use device coordinates for position, not data
; number = number of plot symbols to plot or length of line (D=1)
; usersym = 2-D array of vertices, cf. usersym in IDL manual.
; (/USERSYM =square, default is to use existing USERSYM definition)
; /fill = flag to fill the usersym
; /left_legend = flag to place legend snug against left side of plot
; window (D)
; /right_legend = flag to place legend snug against right side of plot
; window. If /right,pos=[x,y], then x is position of RHS and
; text runs right-to-left.
; /top_legend = flag to place legend snug against top of plot window (D)
; /bottom = flag to place legend snug against bottom of plot window
; /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
;
; If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
; supplied as scalars, then the scalar value is set for every line or
; symbol in the legend.
; Outputs:
; legend to current plot device
; OPTIONAL OUTPUT KEYWORDS:
; corners = 4-element array, like !p.position, of the normalized
; coords for the box (even if box=0): [llx,lly,urx,ury].
; Useful for multi-column or multi-line legends, for example,
; to make a 2-column legend, you might do the following:
; c1_items = ['diamond','asterisk','square']
; c1_psym = [4,2,6]
; c2_items = ['solid','dashed','dotted']
; c2_line = [0,2,1]
; legend,c1_items,psym=c1_psym,corners=c1,box=0
; legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
; c = [c1[0]c2[2],c1[3]>c2[3]]
; plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
; Useful also to place the legend. Here's an automatic way to place
; the legend in the lower right corner. The difficulty is that the
; legend's width is unknown until it is plotted. In this example,
; the legend is plotted twice: the first time in the upper left, the
; second time in the lower right.
; legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
; ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
; xydims = [corners[2]-corners[0],corners[3]-corners[1]]
; ; SAVE WIDTH AND HEIGHT
; chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
; ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
; pos = [!x.window[1]-chdim[0]-xydims[0] $
; ,!y.window[0]+chdim[1]+xydims[1]]
; ; CALCULATE POSITION FOR LOWER RIGHT
; plot,findgen(10) ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
; legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
; ; REDO THE LEGEND IN LOWER RIGHT CORNER
; You can modify the pos calculation to place the legend where you
; want. For example to place it in the upper right:
; pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
; Common blocks:
; none
; Procedure:
; If keyword help is set, call doc_library to print header.
; See notes in the code. Much of the code deals with placement of the
; legend. The main problem with placement is not being
; able to sense the length of a string before it is output. Some crude
; approximations are used for centering.
; Restrictions:
; Here are some things that aren't implemented.
; - An orientation keyword would allow lines at angles in the legend.
; - An array of usersyms would be nice---simple change.
; - An order option to interchange symbols and text might be nice.
; - Somebody might like double boxes, e.g., with box = 2.
; - Another feature might be a continuous bar with ticks and text.
; - There are no guards to avoid writing outside the plot area.
; - There is no provision for multi-line text, e.g., '1st line!c2nd line'
; Sensing !c would be easy, but !c isn't implemented for PostScript.
; A better way might be to simply output the 2nd line as another item
; but without any accompanying symbol or linestyle. A flag to omit
; the symbol and linestyle is linestyle[i] = -1.
; - There is no ability to make a title line containing any of titles
; for the legend, for the symbols, or for the text.
; Side Effects:
; Modification history:
; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
; allow omission of items or omission of both psym and linestyle, add
; corners keyword to facilitate multi-column legends, improve place-
; ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
; add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
; add keyword vectorfont to allow characters in the sym/line column,
; 28 Aug 92, FKK
; add /top, /bottom, /left, /right keywords for automatic placement at
; the four corners of the plot window. The /right keyword forces
; right-to-left printing of menu. 18 Jun 93, FKK
; change default position to data coords and add normal, data, and
; device keywords, 17 Jan 94, FKK
; add /center keyword for positioning, but it is not precise because
; text string lengths cannot be known in advance, 17 Jan 94, FKK
; add interactive positioning with /position keyword, 17 Jan 94, FKK
; allow a legend with just text, no plotting symbols. This helps in
; simply describing a plot or writing assumptions done, 4 Feb 94, FKK
; added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
; David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
; allow scalar specification of keywords, Mar 96, W. Landsman HSTX
; added charthick keyword, June 96, W. Landsman HSTX
; Made keyword names left,right,top,bottom,center longer,
; Aug 16, 2000, Kim Tolbert
; Added ability to have regular text lines in addition to plot legend
; lines in legend. If linestyle is -99 that item is left-justified.
; Previously, only option for no sym/line was linestyle=-1, but then text
; was lined up after sym/line column. 10 Oct 2000, Kim Tolbert
; Make default value of thick = !P.thick W. Landsman Jan. 2001
; Don't overwrite existing USERSYM definition W. Landsman Mar. 2002
; NAME:
; MINF_BRACKET
; PURPOSE:
; Bracket a local minimum of a 1-D function with 3 points,
; EXPLANATION:
; Brackets a local minimum of a 1-d function with 3 points,
; thus ensuring that a minimum exists somewhere in the interval.
; This routine assumes that the function has a minimum somewhere....
; Routine can also be applied to a scalar function of many variables,
; for such case the local minimum in a specified direction is bracketed,
; This routine is called by minF_conj_grad, to bracket minimum in the
; direction of the conjugate gradient of function of many variables
; CALLING EXAMPLE:
; xa=0 & xb=1
; minF_bracket, xa,xb,xc, fa,fb,fc, FUNC_NAME="name" ;for 1-D func.
; or:
; minF_bracket, xa,xb,xc, fa,fb,fc, FUNC="name", $
; POINT=[0,1,1], $
; DIRECTION=[2,1,1] ;for 3-D func.
; INPUTS:
; xa = scalar, guess for point bracketing location of minimum.
; xb = scalar, second guess for point bracketing location of minimum.
; KEYWORDS:
; FUNC_NAME = function name (string)
; Calling mechanism should be: F = func_name( px )
; where:
; px = scalar or vector of independent variables, input.
; F = scalar value of function at px.
; POINT_NDIM = when working with function of N variables,
; use this keyword to specify the starting point in N-dim space.
; Default = 0, which assumes function is 1-D.
; DIRECTION = when working with function of N variables,
; use this keyword to specify the direction in N-dim space
; along which to bracket the local minimum, (default=1 for 1-D).
; (xa,xb,xc) are then relative distances from POINT_NDIM.
; OUTPUTS:
; xa,xb,xc = scalars, 3 points which bracket location of minimum,
; that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists.
; When working with function of N variables
; (xa,xb,xc) are then relative distances from POINT_NDIM,
; in the direction specified by keyword DIRECTION,
; with scale factor given by magnitude of DIRECTION.
; OPTIONAL OUTPUT:
; fa,fb,fc = value of function at 3 points which bracket the minimum,
; again note that fb < fa and fb < fc if minimum exists.
; PROCEDURE:
; algorithm from Numerical Recipes (by Press, et al.), sec.10.1 (p.281).
; MODIFICATION HISTORY:
; Written, Frank Varosi NASA/GSFC 1992.
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; MINF_PARABOLIC
; PURPOSE:
; Minimize a function using Brent's method with parabolic interpolation
; EXPLANATION:
; Find a local minimum of a 1-D function up to specified tolerance.
; This routine assumes that the function has a minimum nearby.
; (recommend first calling minF_bracket, xa,xb,xc, to bracket minimum).
; Routine can also be applied to a scalar function of many variables,
; for such case the local minimum in a specified direction is found,
; This routine is called by minF_conj_grad, to locate minimum in the
; direction of the conjugate gradient of function of many variables.
;
; CALLING EXAMPLES:
; minF_parabolic, xa,xb,xc, xmin, fmin, FUNC_NAME="name" ;for 1-D func.
; or:
; minF_parabolic, xa,xb,xc, xmin, fmin, FUNC="name", $
; POINT=[0,1,1], $
; DIRECTION=[2,1,1] ;for 3-D func.
; INPUTS:
; xa,xb,xc = scalars, 3 points which bracket location of minimum,
; that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists.
; When working with function of N variables
; (xa,xb,xc) are then relative distances from POINT_NDIM,
; in the direction specified by keyword DIRECTION,
; with scale factor given by magnitude of DIRECTION.
; INPUT KEYWORDS:
; FUNC_NAME = function name (string)
; Calling mechanism should be: F = func_name( px )
; where:
; px = scalar or vector of independent variables, input.
; F = scalar value of function at px.
;
; POINT_NDIM = when working with function of N variables,
; use this keyword to specify the starting point in N-dim space.
; Default = 0, which assumes function is 1-D.
; DIRECTION = when working with function of N variables,
; use this keyword to specify the direction in N-dim space
; along which to bracket the local minimum, (default=1 for 1-D).
; (xa, xb, xc, x_min are then relative distances from POINT_NDIM)
; MAX_ITER = maximum allowed number iterations, default=100.
; TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7).
; OUTPUTS:
; xmin = estimated location of minimum.
; When working with function of N variables,
; xmin is the relative distance from POINT_NDIM,
; in the direction specified by keyword DIRECTION,
; with scale factor given by magnitude of DIRECTION,
; so that min. Loc. Pmin = Point_Ndim + xmin * Direction.
; fmin = value of function at xmin (or Pmin).
; PROCEDURE:
; Brent's method to minimize a function by using parabolic interpolation.
; Based on function BRENT in Numerical Recipes in FORTRAN (Press et al.
; 1992), sec.10.2 (p. 397).
; MODIFICATION HISTORY:
; Written, Frank Varosi NASA/GSFC 1992.
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; MINF_PARABOL_D
; PURPOSE:
; Minimize a function using a modified Brent's method with derivatives
; EXPLANATION:
; Based on the procedure DBRENT in Numerical Recipes by Press et al.
; Finds a local minimum of a 1-D function up to specified tolerance,
; using the first derivative of function in the algorithm.
; This routine assumes that the function has a minimum nearby.
; (recommend first calling minF_bracket, xa,xb,xc, to bracket minimum).
; Routine can also be applied to a scalar function of many variables,
; for such case the local minimum in a specified direction is found,
; This routine is called by minF_conj_grad, to locate minimum in the
; direction of the conjugate gradient of function of many variables.
;
; CALLING EXAMPLES:
; minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC_NAME="name" ;for 1-D func.
; or:
; minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC="name", $
; POINT=[0,1,1], $
; DIRECTION=[2,1,1] ;for 3-D func.
; INPUTS:
; xa,xb,xc = scalars, 3 points which bracket location of minimum,
; that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists.
; When working with function of N variables
; (xa,xb,xc) are then relative distances from POINT_NDIM,
; in the direction specified by keyword DIRECTION,
; with scale factor given by magnitude of DIRECTION.
; KEYWORDS:
; FUNC_NAME = function name (string)
; Calling mechanism should be: F = func_name( px, gradient )
; where:
; px = scalar or vector of independent variables, input.
; F = scalar value of function at px.
; gradient = derivative of function, a scalar if 1-D,
; a gradient vector if N-D,
; (should only be computed if arg. is present).
;
; POINT_NDIM = when working with function of N variables,
; use this keyword to specify the starting point in N-dim space.
; Default = 0, which assumes function is 1-D.
; DIRECTION = when working with function of N variables,
; use this keyword to specify the direction in N-dim space
; along which to bracket the local minimum, (default=1 for 1-D).
; (xa, xb, xc, x_min are then relative distances from POINT_NDIM)
; MAX_ITER = maximum allowed number iterations, default=100.
; TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7).
;
; OUTPUTS:
; xmin = estimated location of minimum.
; When working with function of N variables,
; xmin is the relative distance from POINT_NDIM,
; in the direction specified by keyword DIRECTION,
; with scale factor given by magnitude of DIRECTION,
; so that min. Loc. Pmin = Point_Ndim + xmin * Direction.
; fmin = value of function at xmin (or Pmin).
; PROCEDURE:
; Brent's method to minimize a function by using parabolic interpolation
; and using first derivative of function,
; from Numerical Recipes (by Press, et al.), sec.10.3 (p.287),
; MODIFICATION HISTORY:
; Written, Frank Varosi NASA/GSFC 1992.
; NAME:
; MKHDR
; PURPOSE:
; Make a minimal primary (or IMAGE extension) FITS header
; EXPLANATION:
; If an array is supplied, then the created FITS header will be
; appropriate to the supplied array. Otherwise, the user can specify
; the dimensions and datatype.
;
; CALLING SEQUENCE:
; MKHDR, header ;Prompt for image size and type
; or
; MKHDR, header, im, [ /IMAGE, /EXTEND ]
; or
; MKHDR, header, type, naxisx, [/IMAGE, /EXTEND ]
;
; OPTIONAL INPUTS:
; IM - If IM is a vector or array then the header will be made
; appropriate to the size and type of IM. IM does not have
; to be the actual data; it can be a dummy array of the same
; type and size as the data. Set IM = '' to create a dummy
; header with NAXIS = 0.
; TYPE - If more than 2 parameters are supplied, then the second parameter
; is interpreted as an integer giving the IDL datatype e.g.
; 1 - LOGICAL*1, 2 - INTEGER*2, 4 - REAL*4, 3 - INTEGER*4
; NAXISX - Vector giving the size of each dimension (NAXIS1, NAXIS2,
; etc.).
;
; OUTPUT:
; HEADER - image header, (string array) with required keywords
; BITPIX, NAXIS, NAXIS1, ... Further keywords can be added
; to the header with SXADDPAR.
;
; OPTIONAL INPUT KEYWORDS:
; /IMAGE = If set, then a minimal header for a FITS IMAGE extension
; is created. An IMAGE extension header is identical to
; a primary FITS header except the first keyword is
; 'XTENSION' = 'IMAGE' instead of 'SIMPLE ' = 'T'
; /EXTEND = If set, then the keyword EXTEND is inserted into the file,
; with the value of "T" (true). The EXTEND keyword must be
; included in a primary header, if the FITS file contains
; extensions.
;
; RESTRICTIONS:
; (1) MKHDR should not be used to make an STSDAS header or a FITS
; ASCII or Binary Table extension header. Instead use
;
; SXHMAKE - to create a minimal STSDAS header
; FXBHMAKE - to create a minimal FITS binary table header
; FTCREATE - to create a minimal FITS ASCII table header
;
; (2) Any data already in the header before calling MKHDR
; will be destroyed.
; EXAMPLE:
; Create a minimal FITS header, Hdr, for a 30 x 40 x 50 INTEGER*2 array
;
; IDL> mkhdr, Hdr, 2, [30,40,50]
;
; Alternatively, if the array already exists as an IDL variable, Array,
;
; IDL> mkhdr, Hdr, Array
;
; PROCEDURES CALLED:
; SXADDPAR, GET_DATE
;
; REVISION HISTORY:
; Written November, 1988 W. Landsman
; May, 1990, Adapted for IDL Version 2.0, J. Isensee
; Aug, 1997, Use SYSTIME(), new DATE format W. Landsman
; Converted to IDL V5.0 W. Landsman September 1997
; Allow unsigned data types W. Landsman December 1999
; Set BZERO = 0 for unsigned integer data W. Landsman January 2000
; EXTEND keyword must immediately follow last NAXISi W. Landsman Sep 2000
; Add FITS definition COMMENT to primary headers W. Landsman Oct. 2001
; Allow (nonstandard) 64 bit integers W. Landsman Feb. 2003
; NAME:
; MODFITS
; PURPOSE:
; Modify a FITS file by updating the header and/or data array.
; EXPLANATION:
; Since August 2003, the size of the supplied FITS header or data array
; does not need to match the size of the existing header or data array.
;
; CALLING SEQUENCE:
; MODFITS, Filename_or_fcb, Data, [ Header, EXTEN_NO =, ERRMSG=]
;
; INPUTS:
; FILENAME/FCB = Scalar string containing either the name of the FITS file
; to be modified, or the IO file control block returned after
; opening the file with FITS_OPEN,/UPDATE. The explicit
; use of FITS_OPEN can save time if many extensions in a
; single file will be updated.
;
; DATA - data array to be inserted into the FITS file. Set DATA = 0
; to leave the data portion of the FITS file unmodified
;
; HEADER - FITS header (string array) to be updated in the FITS file.
;
; OPTIONAL INPUT KEYWORDS:
; EXTEN_NO - scalar integer specifying the FITS extension to modified. For
; example, specify EXTEN = 1 or /EXTEN to modify the first
; FITS extension.
; OPTIONAL OUTPUT KEYWORD:
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
;
; EXAMPLES:
; (1) Modify the value of the DATE keyword in the primary header of a
; file TEST.FITS.
;
; IDL> h = headfits('test.fits') ;Read primary header
; IDL> sxaddpar,h,'DATE','2001-03-23' ;Modify value of DATE
; IDL> modfits,'test.fits',0,h ;Update header only
;
; (2) Replace the values of the primary image array in 'test.fits' with
; their absolute values
;
; IDL> im = readfits('test.fits') ;Read image array
; IDL> im = abs(im) ;Take absolute values
; IDL> modfits,'test.fits',im ;Update image array
;
; (3) Add some HISTORY records to the FITS header in the first extension
; of a file 'test.fits'
;
; IDL> h = headfits('test.fits',/ext) ;Read first extension hdr
; IDL> sxaddhist,['Comment 1','Comment 2'],h
; IDL> modfits,'test.fits',0,h,/ext ;Update extension hdr
;
; (4) Change 'OBSDATE' keyword to 'OBS-DATE' in every extension in a
; FITS file. Explicitly open with FITS_OPEN to save compute time.
;
; fits_open,'test.fits',io,/update ;Faster to explicity open
; for i = 1,nextend do begin ;Loop over extensions
; fits_read,io,0,h,/header_only,exten_no=i,/No_PDU ;Get header
; date= sxpar(h,'OBSDATE') ;Save keyword value
; sxaddpar,h,'OBS-DATE',date,after='OBSDATE'
; sxdelpar,h,'OBSDATE' ;Delete bad keyword
; modfits,io,0,h,exten_no=i ;Update header
; endfor
;
; Note the use of the /No_PDU keyword in the FITS_READ call -- one
; does *not* want to append the primary header, if the STScI
; inheritance convention is adopted.
;
; NOTES:
; Use the BLKSHIFT procedure to shift the contents of the FITS file if
; the new data or header differs in size by more than 2880 bytes from the
; old data or header.
;
; Also see the procedures FXHMODIFY to add a single FITS keyword to a
; header in a FITS files, and FXBGROW to enlarge the size of a binary
; table.
; RESTRICTIONS:
; (1) Cannot be used to modifiy the data in FITS files with random
; groups or variable length binary tables. (The headers in such
; files *can* be modified.)
;
; (2) If a data array but no FITS header is supplied, then MODFITS does
; not check to make sure that the existing header is consistent with
; the new data.
; PROCEDURES USED:
; Functions: IS_IEEE_BIG(), N_BYTES(), SXPAR()
; Procedures: BLKSHIFT, CHECK_FITS, FITS_OPEN, FITS_READ, HOST_TO_IEEE
;
; MODIFICATION HISTORY:
; Written, Wayne Landsman December, 1994
; Converted to IDL V5.0 W. Landsman September 1997
; Fixed possible problem when using WRITEU after READU October 1997
; New and old sizes need only be the same within multiple of 2880 bytes
; Added call to IS_IEEE_BIG() W. Landsman May 1999
; Added ERRMSG output keyword W. Landsman May 2000
; Update tests for incompatible sizes W. Landsman December 2000
; Major rewrite to use FITS_OPEN procedures W. Landsman November 2001
; Add /No_PDU call to FITS_READ call W. Landsman June 2002
; Update CHECKSUM keywords if already present in header, add padding
; if new data size is smaller than old W.Landsman December 2002
; Only check XTENSION value if EXTEN_NO > 1 W. Landsman Feb. 2003
; Correct for unsigned data on little endian machines W. Landsman Apr 2003
; Major rewrite to allow changing size of data or header W.L. Aug 2003
; NAME:
; MOONPOS
; PURPOSE:
; To compute the RA and Dec of the Moon at specified Julian date(s).
;
; CALLING SEQUENCE:
; MOONPOS, jd, ra, dec, dis, geolong, geolat, [/RADIAN ]
;
; INPUTS:
; JD - Julian date, scalar or vector, double precision suggested
;
; OUTPUTS:
; Ra - Apparent right ascension of the moon in DEGREES, referred to the
; true equator of the specified date(s)
; Dec - The declination of the moon in DEGREES
; Dis - The Earth-moon distance in kilometers (between the center of the
; Earth and the center of the Moon).
; Geolong - Apparent longitude of the moon in DEGREES, referred to the
; ecliptic of the specified date(s)
; Geolat - Apparent longitude of the moon in DEGREES, referred to the
; ecliptic of the specified date(s)
;
; The output variables will all have the same number of elements as the
; input Julian date vector, JD. If JD is a scalar then the output
; variables will be also.
;
; OPTIONAL INPUT KEYWORD:
; /RADIAN - If this keyword is set and non-zero, then all output variables
; are given in Radians rather than Degrees
;
; EXAMPLES:
; (1) Find the position of the moon on April 12, 1992
;
; IDL> jdcnv,1992,4,12,0,jd ;Get Julian date
; IDL> moonpos, jd, ra ,dec ;Get RA and Dec of moon
; IDL> print,adstring(ra,dec,1)
; ==> 08 58 45.23 +13 46 6.1
;
; This is within 1" from the position given in the Astronomical Almanac
;
; (2) Plot the Earth-moon distance for every day at 0 TD in July, 1996
;
; IDL> jdcnv,1996,7,1,0,jd ;Get Julian date of July 1
; IDL> moonpos,jd+dindgen(31), ra, dec, dis ;Position at all 31 days
; IDL> plot,indgen(31),dis, /YNOZ
;
; METHOD:
; Derived from the Chapront ELP2000/82 Lunar Theory (Chapront-Touze' and
; Chapront, 1983, 124, 50), as described by Jean Meeus in Chapter 47 of
; ``Astronomical Algorithms'' (Willmann-Bell, Richmond), 2nd edition,
; 1998. Meeus quotes an approximate accuracy of 10" in longitude and
; 4" in latitude, but he does not give the time range for this accuracy.
;
; Comparison of this IDL procedure with the example in ``Astronomical
; Algorithms'' reveals a very small discrepancy (~1 km) in the distance
; computation, but no difference in the position calculation.
;
; This procedure underwent a major rewrite in June 1996, and the new
; calling sequence is *incompatible with the old* (e.g. angles now
; returned in degrees instead of radians).
;
; PROCEDURES CALLED:
; CIRRANGE, ISARRAY(), NUTATE, TEN() - from IDL Astronomy Library
; POLY() - from IDL User's Library
; MODIFICATION HISTORY:
; Written by Michael R. Greason, STX, 31 October 1988.
; Major rewrite, new (incompatible) calling sequence, much improved
; accuracy, W. Landsman Hughes STX June 1996
; Added /RADIAN keyword W. Landsman August 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Use improved expressions for L',D,M,M', and F given in 2nd edition of
; Meeus (very slight change), W. Landsman November 2000
; NAME:
; MULTIPLOT
; PURPOSE:
; Create multiple plots with shared axes.
; EXPLANATION:
; This procedure makes a matrix of plots with *SHARED AXES*, either using
; parameters passed to multiplot or !p.multi in a non-standard way.
; It is good for data with one or two shared axes and retains all the
; versatility of the plot commands (e.g. all keywords and log scaling).
; The plots are connected with the shared axes, which saves space by
; omitting redundant ticklabels and titles. Multiplot does this by
; setting !p.position, !x.tickname and !y.tickname automatically.
; A call (multiplot,/reset) restores original values.
;
; Note: This method may be superseded by future improvements in !p.multi
; by RSI. For now, it's a good way to gang plots together.
; CALLING SEQUENCE:
; multiplot[pmulti][,/help][,/initialize][,/reset][,/rowmajor]
; EXAMPLES:
; multiplot,/help ; print this header.
; ; Then copy & paste, from your xterm, the following lines to test:
;
; x = findgen(100) ; MULTIPLOT
; t=exp(-(x-50)^2/300) ; -------------------------
; erase ; | | |
; u=exp(-x/30) ; | | |
; y = sin(x) ; | UL plot | UR plot |
; r = reverse(y*u) ; | | |
; !p.multi=[0,2,2,0,0] ; | | |
; multiplot ; y-------------------------
; plot,x,y*u,title='MULTIPLOT' ; l| | |
; multiplot & plot,x,r ; a| | |
; multiplot ; b| LL plot | LR plot |
; plot,x,y*t,ytit='ylabels' ; e| | |
; multiplot ; l| | |
; plot,x,y*t,xtit='xlabels' ; s-------------------------
; multiplot,/reset ; xlabels
;
; wait,2 & erase ; TEST
; multiplot,[1,3] ; H------------------------
; plot,x,y*u,title='TEST' ; E| plot #1 |
; multiplot ; I------------------------
; plot,x,y*t,ytit='HEIGHT' ; G| plot #2 |
; multiplot ; H------------------------
; plot,x,r,xtit='PHASE' ; T| plot #3 |
; multiplot,/reset ; ------------------------
; ; PHASE
;
; multiplot,[1,1],/init,/verbose ; one way to return to single plot
; % MULTIPLOT: Initialized for 1x1, plotted across then down (column major).
; OPTIONAL INPUTS:
; pmulti = 2-element or 5-element vector giving number of plots, e.g.,
; multiplot,[1,6] ; 6 plots vertically
; multiplot,[0,4,2,0,0] ; 4 plots along x and 2 along y
; multiplot,[0,4,2,0,1] ; ditto, except rowmajor (down 1st)
; multiplot,[4,2],/rowmajor ; identical to previous line
; OPTIONAL KEYWORDS:
; help = flag to print header
; initialize = flag to begin only---no plotting, just setup,
; e.g., multiplot,[4,2],/init,/verbose & multiplot & plot,x,y
; reset = flag to reset system variables to values prior to /init
; default = flag to restore IDL's default value for system variables
; rowmajor = flag to number plots down column first (D=columnmajor)
; verbose = flag to output informational messages
; Outputs:
; !p.position = 4-element vector to place a plot
; !x.tickname = either '' or else 30 ' ' to suppress ticknames
; !y.tickname = either '' or else 30 ' ' to suppress ticknames
; !p.noerase = 1
; Common blocks:
; multiplot---to hold saved variables and plot counter. See code.
; Side Effects:
; Multiplot sets a number of system variables: !p.position, !p.multi,
; !x.tickname, !y.tickname, !P.noerase---but all can be reset with
; the call: multiplot,/reset
; RESTRICTIONS:
; 1. If you use !p.multi as the method of telling how many plots
; are present, you have to set !p.multi at the beginning each time you
; use multiplot or call multiplot with the /reset keyword.
; 2. There's no way to make an xtitle or ytitle span more than one plot,
; except by adding spaces to shift it or to add it manually with xyouts.
; 3. There is no way to make plots of different sizes; each plot
; covers the same area on the screen or paper.
; PROCEDURE:
; This routine makes a matrix of plots with common axes, as opposed to
; the method of !p.multi where axes are separated to allow labels.
; Here the plots are joined and labels are suppressed, except at the
; left edge and the bottom. You tell multiplot how many plots to make
; using either !p.multi (which is then reset) or the parameter pmulti.
; However, multiplot keeps track of the position by itself because
; !p.multi interacts poorly with !p.position.
; MODIFICATION HISTORY:
; write, 21-23 Mar 94, Fred Knight (knight@ll.mit.edu)
; alter plot command that sets !x.window, etc. per suggestion of
; Mark Hadfield (hadfield@storm.greta.cri.nz), 7 Apr 94, FKK
; add a /default keyword restore IDL's default values of system vars,
; 7 Apr 94, FKK
; modify two more sys vars !x(y).tickformat to suppress user-formatted
; ticknames, per suggestion of Mark Hadfield (qv), 8 Apr 94, FKK
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; MWRFITS
; PURPOSE:
; Write all standard FITS data types from input arrays or structures.
;
; CALLING SEQUENCE:
; MWRFITS, Input, Filename, [Header],
; /LSCALE , /ISCALE, /BSCALE,
; /USE_COLNUM, /Silent, /Create, /No_comment, /Version, $
; Alias=, /ASCII, Separator=, Terminator=, Null=,
; /Logical_cols, /Bit_cols, /Nbit_cols,
; Group=, Pscale=, Pzero=
;
; INPUTS:
; Input = Array or structure to be written to FITS file.
;
; -When writing FITS primary data or image extensions
; input should be an array.
; --If data is to be grouped
; the Group keyword should be specified to point to
; a two dimensional array. The first dimension of the
; Group array will be PCOUNT while the second dimension
; should be the same as the last dimension of Input.
; --If Input is undefined, then a dummy primary dataset
; or Image extension is created [This might be done, e.g.,
; to put appropriate keywords in a dummy primary
; HDU].
;
; -When writing an ASCII table extension, Input should
; be a structure array where no element of the structure
; is a structure or array (except see below).
; --A byte array will be written as A field. No checking
; is done to ensure that the values in the byte field
; are valid ASCII.
; --Complex numbers are written to two columns with '_R' and
; '_I' appended to the TTYPE fields (if present). The
; complex number is enclosed in square brackets in the output.
; --Strings are written to fields with the length adjusted
; to accommodate the largest string. Shorter strings are
; blank padded to the right.
;
; -When writing a binary table extension, the input should
; be a structure array with no element of the structure
; being a substructure.
;
; If a structure is specified on input and the output
; file does not exist or the /CREATE keyword is specified
; a dummy primary HDU is created.
;
; Filename = String containing the name of the file to be written.
; By default MWRFITS appends a new extension to existing
; files which are assumed to be valid FITS. The /CREATE
; keyword can be used to ensure that a new FITS file
; is created even if the file already exists.
;
; OUTPUTS:
;
; OPTIONAL INPUTS:
; Header = Header should be a string array. Each element of the
; array is added as a row in the FITS header. No
; parsing is done of this data. MWRFITS will prepend
; required structural (and, if specified, scaling)
; keywords before the rows specified in Header.
; Rows describing columns in the table will be appended
; to the contents of Header.
; Header lines will be extended or truncated to
; 80 characters as necessary.
; If Header is specified then on return Header will have
; the header generated for the specified extension.
;
; OPTIONAL INPUT KEYWORDS:
; ALias= Set up aliases to convert from the IDL structure
; to the FITS column name. The value should be
; a STRARR(2,*) value where the first element of
; each pair of values corresponds to a column
; in the structure and the second is the name
; to be used in the FITS file.
; The order of the alias keyword is compatible with
; use in MRDFITS.
; ASCII - Creates an ASCII table rather than a binary table.
; This keyword may be specified as:
; /ASCII - Use default formats for columns.
; ASCII='format_string' allows the user to specify
; the format of various data types such using the following
; syntax 'column_type:format, column_type:format'. E.g.,
; ASCII='A:A1,I:I6,L:I10,B:I4,F:G15.9,D:G23.17,C:G15.9,M:G23.17'
; gives the default formats used for each type. The TFORM
; fields for the real and complex types indicate will use corresponding
; E and D formats when a G format is specified.
; Note that the length of the field for ASCII strings and
; byte arrays is automatically determined for each column.
; BIT_COLS= An array of indices of the bit columns. The data should
; comprise a byte array with the appropriate dimensions.
; If the number of bits per row (see NBIT_COLS)
; is greater than 8, then the first dimension of the array
; should match the number of input bytes per row.
; BSCALE Scale floats, longs, or shorts to unsigned bytes (see LSCALE)
; CREATE If this keyword is non-zero, then a new FITS file will
; be created regardless of whether the file currently
; exists. Otherwise when the file already exists,
; a FITS extension will be appended to the existing file
; which is assumed to be a valid FITS file.
; GROUP= This keyword indicates that GROUPed FITS data is to
; be generated.
; Group should be a 2-D array of the appropriate output type.
; The first dimension will set the number of group parameters.
; The second dimension must agree with the last dimension
; of the Input array.
; ISCALE Scale floats or longs to short integer (see LSCALE)
; LOGICAL_COLS= An array of indices of the logical column numbers.
; These should start with the first column having index 0.
; The structure element should be an array of characters
; with the values 'T' or 'F'. This is not checked.
; LSCALE Scale floating point numbers to long integers.
; This keyword may be specified in three ways.
; /LSCALE (or LSCALE=1) asks for scaling to be automatically
; determined. LSCALE=value divides the input by value.
; I.e., BSCALE=value, BZERO=0. Numbers out of range are
; given the value of NULL if specified, otherwise they are given
; the appropriate extremum value. LSCALE=(value,value)
; uses the first value as BSCALE and the second as BZERO
; (or TSCALE and TZERO for tables).
; NBIT_COLS= The number of bits actually used in the bit array.
; This argument must point to an array of the same dimension
; as BIT_COLS.
; NO_TYPES If the NO_TYPES keyword is specified, then no TTYPE
; keywords will be created for ASCII and BINARY tables.
; No_comment Do not write comment keywords in the header
; NULL= Value to be written for integers/strings which are
; undefined or unwritable.
; PSCALE= An array giving scaling parameters for the group keywords.
; It should have the same dimension as the first dimension
; of Group.
; PZERO= An array giving offset parameters for the group keywords.
; It should have the same dimension as the first dimension
; of Group.
; Separator= This keyword can be specified as a string which will
; be used to separate fields in ASCII tables. By default
; fields are separated by a blank.
; SILENT Suppress informative messages. Errors will still
; be reported.
; Terminator= This keyword can be specified to provide a string which
; will be placed at the end of each row of an ASCII table.
; No terminator is used when not specified.
; If a non-string terminator is specified (including
; when the /terminator form is used), a new line terminator
; is appended.
; USE_COLNUM When creating column names for binary and ASCII tables
; MWRFITS attempts to use structure field name
; values. If USE_COLNUM is specified and non-zero then
; column names will be generated as 'C1, C2, ... 'Cn'
; for the number of columns in the table.
; Version Print the version number of MWRFITS.
;
; EXAMPLE:
; Write a simple array:
; a=fltarr(20,20)
; mwrfits,a,'test.fits'
;
; Append a 3 column, 2 row, binary table extension to file just created.
; a={name:'M31', coords:(30., 20.), distance:2}
; a=replicate(a, 2);
; mwrfits,a,'test.fits'
;
; Now add on an image extension:
; a=lonarr(10,10,10)
; hdr=("COMMENT This is a comment line to put in the header", $
; "MYKEY = "Some desired keyword value")
; mwrfits,a,'test.fits',hdr
;
; RESTRICTIONS:
; (1) Variable length columns are not supported for anything
; other than simple types (byte, int, long, float, double).
; NOTES:
; This multiple format FITS writer is designed to provide a
; single, simple interface to writing all common types of FITS data.
; Given the number of options within the program and the
; variety of IDL systems available it is likely that a number
; of bugs are yet to be uncovered. If you find an anomaly
; please send a report to:
; Tom McGlynn
; NASA/GSFC Code 660.2
; tam@silk.gsfc.nasa.gov (or 301-286-7743)
;
; PROCEDURES USED:
; FXPAR(), FXADDPAR, IS_IEEE_BIG(), HOST_TO_IEEE
; MODIfICATION HISTORY:
; Version 0.9: By T. McGlynn 1997-07-23
; Initial beta release.
; Dec 1, 1997, Lindler, Modified to work under VMS.
; Version 0.91: T. McGlynn 1998-03-09
; Fixed problem in handling null primary arrays.
; Reconverted to IDL 5.0 format using IDLv4_to_v5
; Version 0.92: T. McGlynn 1998-09-09
; Add no_comment flag and keep user comments on fields.
; Fix handling of bit fields.
; Version 0.93: T. McGlynn 1999-03-10
; Fix table appends on VMS.
; Version 0.93a W. Landsman/D. Schlegel
; Update keyword values in chk_and_upd if data type has changed
; Version 0.94: T. McGlynn 2000-02-02
; Efficient processing of ASCII tables.
; Use G rather than E formats as defaults for ASCII tables
; and make the default precision long enough that transformations
; binary to/from ASCII are invertible.
; Some loop indices made long.
; Fixed some ends to match block beginnings.
; Version 0.95: T. McGlynn 2000-11-06
; Several fixes to scaling. Thanks to David Sahnow for
; documenting the problems.
; Added PCOUNT,GCOUNT keywords to Image extensions.
; Version numbers shown in SIMPLE/XTENSION comments
; Version 0.96: T. McGlynn 2001-04-06
; Changed how files are opened to handle ~ consistently
; Version 1.0: T. McGlynn 2001-12-04
; Unsigned integers,
; 64 bit integers.
; Aliases
; Variable length arrays
; Some code cleanup
; Version 1.1: T. McGlynn 2002-2-18
; Fixed major bug in processing of unsigned integers.
; (Thanks to Stephane Beland)
; Version 1.2: Stephane Beland 2003-03-17
; Fixed problem in creating dummy dataset when passing undefined
; data, caused by an update to FXADDPAR routine.
; Version 1.2.1 Stephane Beland 2003-09-10
; Exit gracefully if write priveleges unavailable
;
;
; NAME:
; OPLOTERROR
; PURPOSE:
; Over-plot data points with accompanying X or Y error bars.
; EXPLANATION:
; For use instead of PLOTERROR when the plotting system has already been
; defined.
;
; CALLING SEQUENCE:
; oploterror, [ x,] y, [xerr], yerr,
; [ /NOHAT, HATLENGTH= , ERRTHICK =, ERRSTYLE=, ERRCOLOR =,
; /LOBAR, /HIBAR, NSKIP = , NSUM = , ... OPLOT keywords ]
; INPUTS:
; X = array of abcissae, any datatype except string
; Y = array of Y values, any datatype except string
; XERR = array of error bar values (along X)
; YERR = array of error bar values (along Y)
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
; /NOHAT = if specified and non-zero, the error bars are drawn
; without hats.
; HATLENGTH = the length of the hat lines used to cap the error bars.
; Defaults to !D.X_VSIZE / 100).
; ERRTHICK = the thickness of the error bar lines. Defaults to the
; THICK plotting keyword.
; ERRSTYLE = the line style to use when drawing the error bars. Uses
; the same codes as LINESTYLE.
; ERRCOLOR = scalar integer (0 - !D.N_TABLE) specifying the color to
; use for the error bars
; NSKIP = Positive Integer specifying the error bars to be plotted.
; For example, if NSKIP = 2 then every other error bar is
; plotted; if NSKIP=3 then every third error bar is plotted.
; Default is to plot every error bar (NSKIP = 1)
; NSUM = Number of points to average over before plotting, default =
; !P.NSUM The errors are also averaged, and then divided by
; sqrt(NSUM). This approximation is meaningful only when the
; neighboring error bars have similar sizes.
;
; /LOBAR = if specified and non-zero, will draw only the -ERR error bars.
; /HIBAR = if specified and non-zero, will draw only the +ERR error bars.
; If neither LOBAR or HIBAR are set _or_ if both are set,
; you will get both error bars. Just specify one if you
; only want one set.
; Any valid keywords to the OPLOT command (e.g. PSYM, YRANGE) are also
; accepted by OPLOTERROR via the _EXTRA facility.
;
; NOTES:
; If only two parameters are input, they are taken as Y and YERR. If only
; three parameters are input, they will be taken as X, Y and YERR,
; respectively.
;
; EXAMPLE:
; Suppose one has X and Y vectors with associated errors XERR and YERR
; and that a plotting system has already been defined:
;
; (1) Overplot Y vs. X with both X and Y errors and no lines connecting
; the points
; IDL> oploterror, x, y, xerr, yerr, psym=3
;
; (2) Like (1) but overplot only the Y errors bars and omits "hats"
; IDL> oploterror, x, y, yerr, psym=3, /NOHAT
;
; (3) Like (2) but suppose one has a positive error vector YERR1, and
; a negative error vector YERR2 (asymmetric error bars)
; IDL> oploterror, x, y, yerr1, psym=3, /NOHAT,/HIBAR
; IDL> oploterror, x, y, yerr2, psym=3, /NOHAT,/LOBAR
;
; PROCEDURE:
; A plot of X versus Y with error bars drawn from Y - YERR to Y + YERR
; and optionally from X - XERR to X + XERR is written to the output device
;
; WARNING:
; This an enhanced version of the procedure OPLOTERR in the standard RSI
; library. It was renamed to OPLOTERROR in June 1998 in the IDL
; Astronomy library.
;
; MODIFICATION HISTORY:
; Adapted from the most recent version of PLOTERR. M. R. Greason,
; Hughes STX, 11 August 1992.
; Added COLOR keyword option to error bars W. Landsman November 1993
; Add ERRCOLOR, use _EXTRA keyword, W. Landsman, July 1995
; Remove spurious call to PLOT_KEYWORDS W. Landsman, August 1995
; OPLOT more than 32767 error bars W. Landsman, Feb 1996
; Added NSKIP keyword W. Landsman, Dec 1996
; Added HIBAR and LOBAR keywords, M. Buie, Lowell Obs., Feb 1998
; Rename to OPLOTERROR W. Landsman June 1998
; Converted to IDL V5.0 W. Landsman June 1998
; Ignore !P.PSYM when drawing error bars W. Landsman Jan 1999
; Handle NSUM keyword correctly W. Landsman Aug 1999
; Check limits for logarithmic axes W. Landsman Nov. 1999
; Work in the presence of NAN values W. Landsman Dec 2000
; Improve logic when NSUM or !P.NSUM is set W. Landsman Jan 2001
; Remove NSUM keyword from PLOTS call W. Landsman March 2001
; Only draw error bars with in XRANGE (for speed) W. Landsman Jan 2002
; Fix Jan 2002 update to work with log plots W. Landsman Jun 2002
; NAME:
; PLANET_COORDS
; PURPOSE:
; Find low or high precision RA and DEC for the planets given a date
;
; EXPLANATION:
; For low precision this routine uses HELIO to get the heliocentric ecliptic
; coordinates of the planets at the given date, then converts these to
; geocentric ecliptic coordinates ala "Astronomical Algorithms" by Jean
; Meeus (1991, p 209). These are then converted to RA and Dec using EULER.
; The accuracy between the years 1800 and 2050 is better than 1 arcminute
; for the terrestial planets, but reaches 10 arcminutes for Saturn.
; Before 1850 or after 2050 the accuracy can get much worse.
;
; For high precision use the /JPL option ito use the full JPL ephemeris.
; CALLING SEQUENCE:
; PLANET_COORDS, DATE, RA, DEC, [ PLANET = , /JD, /JPL]
;
; INPUTS:
; DATE - If /JD is not set, then date is a 3-6 element vector containing
; year,month (1-12), day, and optionally hour, minute, & second.
; If /JD is set then DATE is a Julian date. An advantage of the
; /JD option is that it allows the use of vector dates.
; OUTPUTS:
; RA - right ascension of planet(s), J2000 degrees, double precision
; DEC - declination of planet(s), J2000 degrees, double precision
;
; OPTIONAL INPUT KEYWORD:
; PLANET - scalar string giving name of a planet, e.g. 'venus'. Default
; is to compute coords for all of them (except Earth).
; /JD - If set, then the date parameter should be supplied as Julian date
; JPL - if /JPL set, then PLANET_COORDS will call the procedure
; JPLEPHINTERP to compute positions using the full JPL ephemeris.
; The JPL ephemeris FITS file JPLEPH.405 must exist in either the
; current directory, or in the directory specified by the
; environment variable ASTRO_DATA. Alternatively, the JPL keyword
; can be set to the full path and name of the ephemeris file.
; A copy of the JPL ephemeris FITS file JPLEPH.405 is available in
; http://idlastro.gsfc.nasa.gov/ftp/data/
; EXAMPLES:
; (1) Find the RA, Dec of Venus on 1992 Dec 20
; IDL> planet_coords, [1992,12,20], ra,dec ;Compute for all planets
; IDL> print,adstring(ra[1],dec[1],1) ;Venus is second planet
; ====> RA = 21 05 2.66 Dec = -18 51 45.7
; This position is 37" from the full DE406 ephemeris position of
; RA = 21 05 5.24 -18 51 43.1
;
; (2) Return the current RA and Dec of all 8 planets using JPL ephemeris
; IDL> get_juldate, jd ;Get current Julian Date
; IDL> planet_coords,jd,ra,dec,/jd,/jpl ;Find positions of all planets
; IDL> forprint,adstring(ra,dec,0) ;Display positions
;
; (3) Plot the declination of Mars for every day in the year 2001
; IDL> jdcnv,2001,1,1,0,jd ;Get Julian date of midnight on Jan 1
; Now get Mars RA,Dec for 365 consecutive days
; IDL> planet_coords,jd+indgen(365),ra,dec,/jd, planet = 'mars'
; IDL> plot,indgen(365)+1,dec
; NOTES:
; HELIO is based on the two-body problem and neglects interactions
; between the planets. This is why the worst results are for
; Saturn. Use the /JPL option or the online ephemeris generator
; http://ssd.jpl.nasa.gov/cgi-bin/eph for more accuracy.
;
; The procedure returns astrometric coordinates, i.e. no correction
; for aberration. A correction for light travel time is applied
; when /JPL is set, but not for the default low-precision calculation.
; PROCEDURES USED:
; JULDATE
; EULER, HELIO - if /JPL is not set
; JPLEPHREAD, JPLEPHINTERP - if /JPL is set
; REVISION HISTORY:
; Written P.Plait & W. Landsman August 2000
; Fixed Julian date conversion W. Landsman August 2000
; Added /JPL keyword W. Landsman July 2001
; Allow vector Julian dates with JPL ephemeris W. Landsman December 2002
; NAME:
; PLOTERROR
; PURPOSE:
; Plot data points with accompanying X or Y error bars.
; EXPLANATION:
; This is a greatly enhanced version of the standard IDL Library routine
; PLOTERR
;
; CALLING SEQUENCE:
; ploterror, [ x,] y, [xerr], yerr [, TYPE=, /NOHAT, HATLENGTH= , NSUM =
; ERRTHICK=, ERRSTYLE=, ERRCOLOR=, NSKIP=, .. PLOT keywords]
;
; INPUTS:
; X = array of abcissae.
; Y = array of Y values.
; XERR = array of error bar values (along X)
; YERR = array of error bar values (along Y)
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
; TYPE = type of plot produced. The possible types are:
; TYPE = 0 : X Linear - Y Linear (default)
; TYPE = 1 : X Linear - Y Log
; TYPE = 2 : X Log - Y Linear
; TYPE = 3 : X Log - Y Log
; Actually, if 0 is specified, the XLOG and YLOG keywords
; are used. If these aren't specified, then a linear-linear
; plot is produced. This keyword is available to maintain
; compatibility with the previous version of PLOTERROR.
; /NOHAT = if specified and non-zero, the error bars are drawn
; without hats.
; HATLENGTH = the length of the hat lines in device units used to cap the
; error bars. Defaults to !D.X_VSIZE / 100).
; ERRTHICK = the thickness of the error bar lines. Defaults to the
; THICK plotting keyword.
; ERRSTYLE = the line style to use when drawing the error bars. Uses
; the same codes as LINESTYLE.
; ERRCOLOR = scalar integer (0 - !D.N_TABLE) specifying the color to
; use for the error bars
; NSKIP = Integer specifying the error bars to be plotted. For example,
; if NSKIP = 2 then every other error bar is plotted; if NSKIP=3
; then every third error bar is plotted. Default is to plot
; every error bar (NSKIP = 1)
; NSUM = Number of points to average over before plotting, default=!P.NSUM
; The errors are also averaged, and then divided by sqrt(NSUM).
; This approximation is meaningful only when the neighboring error
; bars have similar sizes. PLOTERROR does not pass the NSUM
; keyword to the PLOT command, but rather computes the binning
; itself using the FREBIN function.
;
; Any valid keywords to the PLOT command (e.g. PSYM, YRANGE) are also
; accepted by PLOTERROR via the _EXTRA facility.
;
; RESTRICTIONS:
; Arrays must not be of type string. There must be enough points to plot.
; If only three parameters are input, they will be taken as X, Y and
; YERR respectively.
;
; PLOTERROR cannot be used for asymmetric error bars. Instead use
; OPLOTERROR with the /LOBAR and /HIBAR keywords.
;
; Any data points with NAN values in the X, Y, or error vectors are
; ignored.
; EXAMPLE:
; Suppose one has X and Y vectors with associated errors XERR and YERR
;
; (1) Plot Y vs. X with both X and Y errors and no lines connecting
; the points
; IDL> ploterror, x, y, xerr, yerr, psym=3
;
; (2) Like (1) but plot only the Y errors bars and omits "hats"
; IDL> ploterror, x, y, yerr, psym=3, /NOHAT
;
; WARNING:
; This an enhanced version of the procedure PLOTERR in the standard IDL
; distribution. It was renamed from PLOTERR to PLOTERROR in June 1998
; in the IDL Astronomy Library to avoid conflict with the RSI procedure.
;
; PROCEDURE:
; A plot of X versus Y with error bars drawn from Y - YERR to Y + YERR
; and optionally from X - XERR to X + XERR is written to the output device
;
; PROCEDURE CALLS:
; FREBIN - used to compute binning if NSUM keyword is present
; MODIFICATION HISTORY:
; William Thompson Applied Research Corporation July, 1986
; DMS, April, 1989 Modified for Unix
; Michael R. Greason ST Systems
; May, 1991 Added most of the plotting keywords, put hats
; on the error bars.
; K. Venkatakrishna Added option to plot xerr, May, 1992
; Michael R. Greason Corrected handling of reversed axes. Aug. 1992
; W. Landsman Use _EXTRA keyword July 1995
; W. Landsman Plot more than 32767 points Feb 1996
; W. Landsman Fix Y scaling when only XRANGE supplied Nov 1996
; W. Landsman Added NSKIP keyword Dec 1996
; W. Landsman Use XLOG, YLOG instead of XTYPE, YTYPE Jan 1998
; W. Landsman Rename to PLOTERROR, OPLOTERROR Jun 1998
; W. Landsman Convert to IDL V5.0 Jun 1998
; W. Landsman Better default scaling when NSKIP supplied Oct 1998
; W. Landsman Ignore !P.PSYM when drawing error bars Jan 1999
; W. Landsman Handle NSUM keyword correctly Aug 1999
; W. Landsman Fix case of /XLOG but no X error bars Oct 1999
; W. Landsman Work in the presence of NAN values Nov 2000
; W. Landsman Improve logic when NSUM or !P.NSUM is set Jan 2001
; W. Landsman Only draw error bars with in XRANGE (for speed) Jan 2002
; W. Landsman Fix Jan 2002 update to work with log plots Jun 2002
; NAME:
; POSANG
; PURPOSE:
; Computes rigorous position angle of source 2 relative to source 1
;
; EXPLANATION:
; Computes the rigorous position angle of source 2 (with given RA, Dec)
; using source 1 (with given RA, Dec) as the center.
;
; CALLING SEQUENCE:
; POSANG, U, RA1, DC1, RA2, DC2, ANGLE
;
; INPUTS:
; U -- Describes units of inputs and output:
; 0: everything radians
; 1: RAx in decimal hours, DCx in decimal
; degrees, ANGLE in degrees
; RA1 -- Right ascension of point 1
; DC1 -- Declination of point 1
; RA2 -- Right ascension of point 2
; DC2 -- Declination of point 2
;
; OUTPUTS:
; ANGLE-- Angle of the great circle containing [ra2, dc2] from
; the meridian containing [ra1, dc1], in the sense north
; through east rotating about [ra1, dc1]. See U above
; for units.
;
; PROCEDURE:
; The "four-parts formula" from spherical trig (p. 12 of Smart's
; Spherical Astronomy or p. 12 of Green' Spherical Astronomy).
;
; EXAMPLE:
; For the star 56 Per, the Hipparcos catalog gives a position of
; RA = 66.15593384, Dec = 33.94988843 for component A, and
; RA = 66.15646079, Dec = 33.96100069 for component B. What is the
; position angle of B relative to A?
;
; IDL> RA1 = 66.15593384/15.d & DC1 = 33.95988843
; IDL> RA2 = 66.15646079/15.d & DC2 = 33.96100069
; IDL> posang,1,ra1,dc1,ra2,dc2, ang
; will give the answer of ang = 21.4 degrees
; NOTES:
; (1) If RA1,DC1 are scalars, and RA2,DC2 are vectors, then ANGLE is a
; vector giving the position angle between each element of RA2,DC2 and
; RA1,DC1. Similarly, if RA1,DC1 are vectors, and RA2, DC2 are scalars,
; then DIS is a vector giving the position angle of each element of RA1,
; DC1 and RA2, DC2. If both RA1,DC1 and RA2,DC2 are vectors then ANGLE
; is a vector giving the position angle between each element of RA1,DC1
; and the corresponding element of RA2,DC2. If then vectors are not the
; same length, then excess elements of the longer one will be ignored.
;
; (2) Note that POSANG is not commutative -- the position angle between
; A and B is theta, then the position angle between B and A is 180+theta
; PROCEDURE CALLS:
; ISARRAY()
; HISTORY:
; Modified from GCIRC, R. S. Hill, RSTX, 1 Apr. 1998
;
; NAME:
; PRECESS
; PURPOSE:
; Precess coordinates from EQUINOX1 to EQUINOX2.
; EXPLANATION:
; For interactive display, one can use the procedure ASTRO which calls
; PRECESS or use the /PRINT keyword. The default (RA,DEC) system is
; FK5 based on epoch J2000.0 but FK4 based on B1950.0 is available via
; the /FK4 keyword.
;
; Use BPRECESS and JPRECESS to convert between FK4 and FK5 systems
; CALLING SEQUENCE:
; PRECESS, ra, dec, [ equinox1, equinox2, /PRINT, /FK4, /RADIAN ]
;
; INPUT - OUTPUT:
; RA - Input right ascension (scalar or vector) in DEGREES, unless the
; /RADIAN keyword is set
; DEC - Input declination in DEGREES (scalar or vector), unless the
; /RADIAN keyword is set
;
; The input RA and DEC are modified by PRECESS to give the
; values after precession.
;
; OPTIONAL INPUTS:
; EQUINOX1 - Original equinox of coordinates, numeric scalar. If
; omitted, then PRECESS will query for EQUINOX1 and EQUINOX2.
; EQUINOX2 - Equinox of precessed coordinates.
;
; OPTIONAL INPUT KEYWORDS:
; /PRINT - If this keyword is set and non-zero, then the precessed
; coordinates are displayed at the terminal. Cannot be used
; with the /RADIAN keyword
; /FK4 - If this keyword is set and non-zero, the FK4 (B1950.0) system
; will be used otherwise FK5 (J2000.0) will be used instead.
; /RADIAN - If this keyword is set and non-zero, then the input and
; output RA and DEC vectors are in radians rather than degrees
;
; RESTRICTIONS:
; Accuracy of precession decreases for declination values near 90
; degrees. PRECESS should not be used more than 2.5 centuries from
; 2000 on the FK5 system (1950.0 on the FK4 system).
;
; EXAMPLES:
; (1) The Pole Star has J2000.0 coordinates (2h, 31m, 46.3s,
; 89d 15' 50.6"); compute its coordinates at J1985.0
;
; IDL> precess, ten(2,31,46.3)*15, ten(89,15,50.6), 2000, 1985, /PRINT
;
; ====> 2h 16m 22.73s, 89d 11' 47.3"
;
; (2) Precess the B1950 coordinates of Eps Ind (RA = 21h 59m,33.053s,
; DEC = (-56d, 59', 33.053") to equinox B1975.
;
; IDL> ra = ten(21, 59, 33.053)*15
; IDL> dec = ten(-56, 59, 33.053)
; IDL> precess, ra, dec ,1950, 1975, /fk4
;
; PROCEDURE:
; Algorithm from Computational Spherical Astronomy by Taff (1983),
; p. 24. (FK4). FK5 constants from "Astronomical Almanac Explanatory
; Supplement 1992, page 104 Table 3.211.1.
;
; PROCEDURE CALLED:
; Function PREMAT - computes precession matrix
;
; REVISION HISTORY
; Written, Wayne Landsman, STI Corporation August 1986
; Correct negative output RA values February 1989
; Added /PRINT keyword W. Landsman November, 1991
; Provided FK5 (J2000.0) I. Freedman January 1994
; Precession Matrix computation now in PREMAT W. Landsman June 1994
; Added /RADIAN keyword W. Landsman June 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Correct negative output RA values when /RADIAN used March 1999
; Work for arrays, not just vectors W. Landsman September 2003
; NAME:
; PUTAST
; PURPOSE:
; Put WCS astrometry parameters into a given FITS header.
;
; CALLING SEQUENCE:
; putast, hdr ;Prompt for all values
; or
; putast, hdr, astr, [EQUINOX =, CD_TYPE = ]
; or
; putast, hdr, cd,[ crpix, crval, ctype], [ EQUINOX =, CD_TYPE = ]
;
; INPUTS:
; HDR - FITS header, string array. HDR will be updated to contain
; the supplied astrometry.
; ASTR - IDL structure containing values of the astrometry parameters
; CDELT, CRPIX, CRVAL, CTYPE, LONGPOLE, PROJP1, and PROJP2
; See EXTAST.PRO for more info about the structure definition
; or
; CD - 2 x 2 array containing the astrometry parameters CD1_1 CD1_2
; CD2_1 CD2_2
; in units of DEGREES/PIXEL
; CRPIX - 2 element vector giving X and Y coord of reference pixel
; BE SURE THE COORDINATES IN CRPIX ARE GIVEN IN FITS STANDARD
; (e.g. first pixel in image is [1,1] ) AND NOT IDL STANDARD
; (first pixel in image is [0,0]
; CRVAL - 2 element vector giving R.A. and DEC of reference pixel
; in degrees
; CTYPE - 2 element string vector giving projection types for the two axes.
; For example, to specify a tangent projection one should set
; ctype = ['RA---TAN','DEC--TAN']
;
; OUTPUTS:
; HDR - FITS header now contains the updated astrometry parameters
; A brief HISTORY record is also added.
;
; OPTIONAL KEYWORD INPUTS:
; EQUINOX - numeric scalar giving the year of equinox of the reference
; coordinates. Default (if EQUINOX keyword is not already
; present in header) is 2000.
;
; CD_TYPE - Integer scalar, either 0, 1 or 2 specifying how the CD matrix
; is to be written into the header
; (0) write PCn_m values along with CDELT values
; (1) convert to rotation and write as a CROTA2 value (+ CDELT)
; (2) as CDn_m values (IRAF standard)
;
; All three forms are valid representations according to Greisen &
; Calabretta (2002, A&A, 395, 1061), also available at
; http://www.aoc.nrao.edu/~egreisen/) although form (0) is preferred.
; Form (1) is the former AIPS standard and is now deprecated.
; If CD_TYPE is not supplied, PUTAST will try to determine the
; type of astrometry already in the header. If there is no
; astrometry in the header then the default is CD_TYPE = 2
; NOTES:
; The recommended use of this procedure is to supply an astrometry
; structure.
;
; PROMPTS:
; If only a header is supplied, the user will be prompted for a plate
; scale, the X and Y coordinates of a reference pixel, the RA and
; DEC of the reference pixel, the equinox of the RA and Dec and a
; rotation angle.
;
; PROCEDURES USED:
; GETOPT(), GET_COORDS, GET_EQUINOX, SXADDPAR, SXPAR(), TAG_EXIST(),
; ZPARCHECK
; REVISION HISTORY:
; Written by W. Landsman 9-3-87
; Major rewrite, use new astrometry structure March, 1994
; Use both CD and CDELT to get plate scale for CD_TYPE=1 September 1995
; Use lower case for FITS keyword Comments W.L. March 1997
; Fixed for CD_TYPE=1 and CDELT = [1.0,1.0] W.L September 1997
; Default value of CD_TYPE is now 2, Use GET_COORDS to read coordinates
; to correct -0 problem W.L. September 1997
; Update CROTA1 if it already exists W.L. October 1997
; Convert rotation to degrees for CD_TYPE = 1 W. L. June 1998
; Convert to IDL V5.0 W.L. June 1998
; Accept CD_TYPE = 0 keyword input W.L October 1998
; Remove reference to obsolete !ERR W.L. February 2000
; No longer support CD001001 format, write default tangent CTYPE value
; consistent conversion between CROTA and CD matrix W.L. October 2000
; Use GET_EQUINOX to get equinox value W.L. January 2001
; Update CTYPE keyword if previous value is 'LINEAR' W.L. July 2001
; Use SIZE(/TNAME) instead of DATATYPE() W.L. November 2001
; Allow direct specification of CTYPE W.L. June 2002
; Don't assume celestial coordinates W. Landsman April 2003
; Make default CD_TYPE = 2 W. Landsman September 2003
; NAME:
; RDFLOAT
; PURPOSE:
; Quickly read a numeric ASCII data file into IDL floating/double vectors.
; EXPLANATION:
; Columns of data may be separated by tabs or spaces. This
; program is fast but is restricted to data files where all columns can
; be read as floating point (or all double precision).
;
; Use READCOL if greater flexibility is desired. Use READFMT to read a
; fixed-format ASCII file. Use FORPRINT to print columns of data.
;
; CALLING SEQUENCE:
; RDFLOAT, name, v1, [ v2, v3, v4, v5, ... v19]
; COLUMNS, /DOUBLE, SKIPLINE = , NUMLINE = ]
;
; INPUTS:
; NAME - Name of ASCII data file, scalar string. In VMS, an extension of
; .DAT is assumed, if not supplied.
;
; OPTIONAL INPUT KEYWORDS:
; COLUMNS - Numeric scalar or vector specifying which columns in the file
; to read. For example, if COLUMNS = [3,7,11] then the first
; output variable (v1) would contain column 3, the second would
; contain column 7 and the third would contain column 11. If
; the number of elements in the COLUMNS vector is less than the
; number of output parameters, then consecutive columns are
; implied. For example, if 3 output parameters are supplied
; (v1,v2,v3) and COLUMNS = 3, then columns 3,4, and 5 will be
; read.
; SKIPLINE - Integer scalar specifying number of lines to skip at the top
; of file before reading. Default is to start at the first line.
; NUMLINE - Integer scalar specifying number of lines in the file to read.
; Default is to read the entire file
; /DOUBLE - If this keyword is set, then all variables are read in as
; double precision.
; /SILENT - Set this keyword to suppress any informative messages
;
; OUTPUTS:
; V1,V2,V3,...V19 - IDL vectors to contain columns of data.
; Up to 19 columns may be read. All output vectors are of type
; float, unless the /DOUBLE keyword is set,
;
; EXAMPLES:
; Each row in a file 'position.dat' contains a star number and 6 columns
; of data giving an RA and Dec in sexigesimal format. Read into IDL
; variables.
;
; IDL> rdfloat,'position.dat',ID,hr,min,sec,deg,dmin,dsec
;
; All output vectors will be floating point. To only read the
; declination vectors (Deg,dmin,dsec)
;
; IDL> rdfloat,'position.dat',deg,dmin,dsec,col=4
;
; RESTRICTIONS:
; (1) All rows in the file must be formatted identically (except for
; those skipped by SKIPLINE). RDFLOAT reads the first line of
; the data (after SKIPLINE) to determine the number of columns of
; data.
; (2) Cannot be used to read strings
; PROCEDURES USED:
; NUMLINES()
; REVISION HISTORY:
; Written W. Landsman September 1995
; Call NUMLINES() function February 1996
; Read up to 19 columns August 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Allow to skip more than 32767 lines W. Landsman June 2001
; Added /SILENT keyword W. Landsman March 2002
; Added COLUMNS keyword, use STRSPLIT W. Landsman May 2002
; Use SKIP_LUN if V5.6 or later W. Landsamn Nov 2002
; NAME:
; RDPLOT
;
; PURPOSE:
; Like CURSOR but with a full-screen cursor and continuous readout option
; EXPLANATION:
; This program is designed to essentially mimic the IDL CURSOR command,
; but with the additional options of continuously printing out the data
; values of the cursor's position, and using a full-screen cursor rather
; than a small cross cursor. The Full screen cursor uses OPLOT and
; X-windows graphics masking to emulate the cursor.
; One difference is that IF the PRINT keyword is set but the DOWN, WAIT,
; or CHANGE keywords are not set, then the leftmost mouse button will
; print a "newline" line-feed, but not exit.
;
; CALLING SEQUENCE:
; RDPLOT, [X, Y, WaitFlag], [/DATA, /DEVICE, /NORMAL,
; /NOWAIT, /WAIT, /DOWN, /CHANGE, ERR=,
; PRINT=, XTITLE=, YTITLE=, XVALUES=, YVALUES=,
; /FULLCURSOR, /NOCLIP, LINESTYLE=, THICK=, COLOR=, /CROSS]
;
; REQUIRED INPUTS:
; None.
;
; OPTIONAL INPUTS:
; WAITFLAG = Uses the same table as the intrinsic CURSOR command, But note
; that unlike the CURSOR command, there is no UP keyword.
; WaitFlag=0 sets the NOWAIT keyword
; WaitFlag=1 sets the WAIT keyword {default}
; WaitFlag=2 sets the CHANGE keyword
; WaitFlag=3 sets the DOWN keyword
;
; OPTIONAL OUTPUTS:
; X - a named variable to receive the final cursor X position, scalar
; Y - a named variable to receive the final cursor Y position, scalar
; OPTIONAL KEYWORD INPUT PARAMETERS:
; /DATA = Data coordinates are displayed and returned.
; /DEVICE = device coordinates are displayed and returned.
; /NORMAL = normal coordinates are displayed and returned.
; Default is to use DATA coordinates if available (see notes).
; /NOWAIT = if non-zero the routine will immediately return the cursor's
; present position.
; WAIT = if non-zero will wait for a mouse key click before returning. If
; cursor key is already down, then procedure immediately exits.
; DOWN = equivalent to WAIT *except* that if the mouse key is already down
; when the procedure is called, the procedure will wait until the mouse
; key is clicked down again.
; CHANGE = returns when the mouse is moved OR a key is clicked up or down.
; PRINT = if non-zero will continuously print out (at the terminal) the data
; values of the cursor's position. If PRINT>1, program will printout a
; brief header describing the mouse button functions. However, note that
; the button functions are overridden if any of the DOWN, WAIT, mouse
; or CHANGE values are non-zero.
; XTITLE = label used to describe the values of the abscissa if PRINT>0.
; YTITLE = label used to describe the values of the ordinate if PRINT>0.
; XVALUES = a vector corresponding to the values to be printed when the
; PRINT keyword is set. This allows the user the option of printing
; out other values rather than the default X coordinate position of
; the cursor. E.g., if XVALUES is a string vector of dates such as
; ['May 1', 'May 2', ...], then those dates will be printed rather than
; the X value of the cursor's position: if X=1 then 'May 2' would be
; printed, etc. This requires that the values of the X coordinate read
; by the cursor must be positive (can't access negative elements).
; If XVALUES=-1, then NO values for X will be printed.
; YVALUES = analagous to the XVALUES keyword.
; FULLCURSOR = if non-zero default cursor is blanked out and full-screen
; (or full plot window, depending on the value of NOCLIP) lines are
; drawn; their intersecton is centered on the cursor position.
; NOCLIP = if non-zero will make a full-screen cursor, otherwise it will
; default to the value in !P.NOCLIP.
; LINESTYLE = style of line that makes the full-screen cursor.
; THICK = thickness of the line that makes the full-screen cursor.
; COLOR = color of the full-screen cursor.
; CROSS = if non-zero will show the regular cross AND full screen cursors.
;
; OPTIONAL KEYWORD OUTPUT PARAMETER:
; ERR = returns the most recent value of the !mouse.button value.
;
; NOTES:
; Note that this procedure does not allow the "UP" keyword/flag...which
; doesn't seem to work too well in the origianl CURSOR version anyway.
;
; If a data coordinate system has not been established, then RDPLOT will
; create one identical to the device coordinate system. Note that this
; kluge is required even if the user specified /NORMAL coordinates, since
; RDPLOT makes use of the OPLOT procedure. This new data coordinate system
; is effectively "erased" (!X.CRange and !Y.CRange are both set to zero)
; upon exit of the routine so as to not change the plot status from the
; user's point of view.
;
; Only tested on X-windows systems. If this program is interrupted, the
; graphics function might be left in a non-standard state; in that case,
; run the program RESET_RDPLOT to return the standard graphics functions,
; or type the command: DEVICE, /CURSOR_CROSS, SET_GRAPHICS=3, BYPASS=0
;
; BUGS:
; It is assumed that the current background of the plot is correctly
; defined by the value in !P.Background. Otherwise, the color of the
; long cursor probably will not be correct. Sometimes the color doesn't
; work anyway, and I'm not sure why.
;
; There may be some cases (e.g., when THICK>1 and NOCLIP=0) when the
; full-screen cursor is not correctly erased, leaving "ghost images" on the
; plot. It just seems that the screen updates get slow or the positions
; ambiguous with a thick line and the cursor off the plot.
;
; PROCEDURE:
; Basically is a bells-n-whistles version of the CURSOR procedure. All
; the details are covered in the above discussion of the keywords.
;
; EXAMPLE (a silly, but informative one):
; Months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', $
; 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
; plot, indgen(12), xrange=[-5, 15]
; rdplot, /FULL, /PRINT, XTITLE='Month: ', YTITLE='Y-value per month = ', $
; xvalues=Months
;
; MODIFICATION HISTORY:
; Written (originally named CURFULL) by J.Wm.Parker 1993 Nov 22
; Created data coordinates if not already present, W. Landsman Nov. 93
; Added continuous printout of data values, COLOR and FULLCURSOR keywords
; (so that default is that it acts just like the cursor command).
; Changed name from CURFULL to RDPLOT. J.Wm.Parker 1994 Apr 20
; Modified (with some translation table assistance from the IDL support
; group) to correctly plot the crosshair with the desired IDL
; color using the device's translation table to determine the XOR
; function and using the BYPASS function. Added the RESET_RDPLOT
; procedure to cleanup crashes that might occur while running
; RDPLOT. Other minor changes/bug fixes. J.Wm.Parker 1994 May 21
; Modified DOWN, WAIT, CHANGE functions to behave more similar to the
; generic CURSOR procedure. J.Wm.Parker 1995 April 24
; Added XVALUES, YVALUES keywords and cleanup. J.Wm.Parker 1995 April 24
; Convert to IDL V5.0, W. Landsman July 1998
; Change !D.NCOLORS to !D.TABLE_SIZE for 24 bit displays W. Landsman May 2000
; Skip translation table for TrueColor visuals W. Landsman March 2001
; NAME:
; READCOL
; PURPOSE:
; Read a free-format ASCII file with columns of data into IDL vectors
; EXPLANATION:
; Lines of data not meeting the specified format (e.g. comments) are
; ignored. Columns may be separated by commas, tabs or spaces.
;
; Use READFMT to read a fixed-format ASCII file. Use RDFLOAT for
; much faster I/O (but less flexibility). Use FORPRINT to write
; columns of data (inverse of READCOL).
;
; CALLING SEQUENCE:
; READCOL, name, v1, [ v2, v3, v4, v5, ... v25 , COMMENT=
; DELIMITER= ,FORMAT = , /DEBUG , /SILENT , SKIPLINE = , NUMLINE = ]
;
; INPUTS:
; NAME - Name of ASCII data file, scalar string. In VMS, an extension of
; .DAT is assumed, if not supplied.
;
; OPTIONAL INPUT KEYWORDS:
; FORMAT - scalar string containing a letter specifying an IDL type
; for each column of data to be read. Allowed letters are
; A - string data, B - byte, D - double precision, F- floating
; point, I - integer, L - longword, Z - longword hexadecimal,
; and X - skip a column.
;
; Columns without a specified format are assumed to be floating
; point. Examples of valid values of FMT are
;
; 'A,B,I' ;First column to read as a character string, then
; 1 column of byte data, 1 column integer data
; 'L,L,L,L' ;Four columns will be read as longword arrays.
; ' ' ;All columns are floating point
;
; If a FORMAT keyword string is not supplied, then all columns are
; assumed to be floating point.
;
; /SILENT - Normally, READCOL will display each line that it skips over.
; If SILENT is set and non-zero then these messages will be
; suppressed.
; /DEBUG - If this keyword is non-zero, then additional information is
; printed as READCOL attempts to read and interpret the file.
; COMMENT - single character specifying comment signal. Any line
; beginning with this character will be skipped. Default is
; no comment lines.
; DELIMITER - single character specifying delimiter used to separate
; columns. Default is either a comma, tab, or a blank.
; SKIPLINE - Scalar specifying number of lines to skip at the top of file
; before reading. Default is to start at the first line.
; NUMLINE - Scalar specifying number of lines in the file to read.
; Default is to read the entire file
;
; OUTPUTS:
; V1,V2,V3,...V25 - IDL vectors to contain columns of data.
; Up to 25 columns may be read. The type of the output vectors
; are as specified by FORMAT.
;
; EXAMPLES:
; Each row in a file position.dat contains a star name and 6 columns
; of data giving an RA and Dec in sexigesimal format. Read into IDL
; variables. (NOTE: The star names must not include the delimiter
; as a part of the name, no spaces or commas as default.)
;
; IDL> FMT = 'A,I,I,F,I,I,F'
; IDL> READCOL,'position.dat',F=FMT,name,hr,min,sec,deg,dmin,dsec
;
; The HR,MIN,DEG, and DMIN variables will be integer vectors.
;
; Alternatively, all except the first column could be specified as
; floating point.
;
; IDL> READCOL,'position.dat',F='A',name,hr,min,sec,deg,dmin,dsec
;
; To read just the variables HR,MIN,SEC
; IDL> READCOL,'position.dat',F='X,I,I,F',HR,MIN,SEC
;
; RESTRICTIONS:
; This procedure is designed for generality and not for speed.
; If a large ASCII file is to be read repeatedly, it may be worth
; writing a specialized reader.
;
; Columns to be read as strings must not contain the delimiter character
; (i.e. commas or spaces by default). Either change the default
; delimiter with the DELIMITER keyword, or use READFMT to read such files.
;
; Numeric values are converted to specified format. For example,
; the value 0.13 read with an 'I' format will be converted to 0.
;
; PROCEDURES CALLED
; GETTOK(), NUMLINES(), REPCHR(), STRNUMBER()
;
; REVISION HISTORY:
; Written W. Landsman November, 1988
; Modified J. Bloch June, 1991
; (Fixed problem with over allocation of logical units.)
; Added SKIPLINE and NUMLINE keywords W. Landsman March 92
; Read a maximum of 25 cols. Joan Isensee, Hughes STX Corp., 15-SEP-93.
; Call NUMLINES() function W. Landsman Feb. 1996
; Added DELIMITER keyword W. Landsman Nov. 1999
; Fix indexing typos (i for k) that mysteriously appeared W. L. Mar. 2000
; Hexadecimal support added. MRG, RITSS, 15 March 2000.
; Default is comma or space delimiters as advertised W.L. July 2001
; Faster algorithm, use STRSPLIT if V5.3 or later W.L. May 2002
; Restore default recognition of tab delimiter W.L. June 2002
; Use SKIP_LUN if V5.6 or later W.L. Nov. 2002
; NAME:
; READFITS
; PURPOSE:
; Read a FITS file into IDL data and header variables.
; EXPLANATION:
; See http://idlastro.gsfc.nasa.gov/fitsio.html for other ways of
; reading FITS files with IDL.
;
; CALLING SEQUENCE:
; Result = READFITS( Filename/Fileunit,[ Header, heap, /NOSCALE, EXTEN_NO=,
; NSLICE=, /SILENT, STARTROW =, NUMROW = , HBUFFER=,
; /COMPRESS, /No_Unsigned ] )
;
; INPUTS:
; Filename = Scalar string containing the name of the FITS file
; (including extension) to be read. If the filename has
; a *.gz extension, it will be treated as a gzip compressed
; file. If it has a .Z extension, it will be treated as a
; Unix compressed file.
; OR
; Fileunit - A scalar integer specifying the unit of an already opened
; FITS file. The unit will remain open after exiting
; READFITS(). There are two possible reasons for choosing
; to specify a unit number rather than a file name:
; (1) For a FITS file with many extensions, one can move to the
; desired extensions with FXPOSIT() and then use READFITS(). This
; is more efficient that repeatedly starting at the beginning of
; the file.
; (2) For reading a FITS file across a Web http: address after opening
; the unit with the SOCKET procedure (IDL V5.4 or later,
; Unix and Windows only)
; OUTPUTS:
; Result = FITS data array constructed from designated record.
; If the specified file was not found, then Result = -1
;
; OPTIONAL OUTPUT:
; Header = String array containing the header from the FITS file.
; If you don't need the header, then the speed may be improved by
; not supplying this parameter.
; heap = For extensions, the optional heap area following the main
; data array (e.g. for variable length binary extensions).
;
; OPTIONAL INPUT KEYWORDS:
; /CHECKSUM - If set, then READFITS() will call FITS_TEST_CHECKSUM to
; verify the data integrity if CHECKSUM keywords are present
; in the FITS header. Cannot be used with the NSLICE, NUMROW
; or STARTROW keywords, since verifying the checksum requires
; that all the data be read. See FITS_TEST_CHECKSUM() for more
; information.
;
; /COMPRESS - Signal that the file is gzip compressed. By default,
; READFITS will assume that if the file name extension ends in
; '.gz' then the file is gzip compressed. The /COMPRESS keyword
; is required only if the the gzip compressed file name does not
; end in '.gz'. Keyword only valid for IDL V5.3 or later.
;
; EXTEN_NO - positive scalar integer specifying the FITS extension to
; read. For example, specify EXTEN = 1 or /EXTEN to read the
; first FITS extension.
;
; /NOSCALE - If present and non-zero, then the ouput data will not be
; scaled using the optional BSCALE and BZERO keywords in the
; FITS header. Default is to scale.
;
; /NO_UNSIGNED -By default, if the header indicates an unsigned integer
; (BITPIX = 16, BZERO=2^15, BSCALE=1) then FITS_READ will output
; an IDL unsigned integer data type (UINT). But if /NO_UNSIGNED
; is set, then the data is converted to type LONG.
;
; NSLICE - Non-negative integer scalar specifying which N-1 dimensional
; slice of a N-dimensional array to read. For example, if the
; primary image of a file 'wfpc.fits' contains a 800 x 800 x 4
; array, then
;
; IDL> im = readfits('wfpc.fits',h, nslice=2)
; is equivalent to
; IDL> im = readfits('wfpc.fits',h)
; IDL> im = im[*,*,2]
; but the use of the NSLICE keyword is much more efficient.
;
; NUMROW - Scalar non-negative integer specifying the number of rows
; of the image or table extension to read. Useful when one
; does not want to read the entire image or table.
;
; POINT_LUN - Position (in bytes) in the FITS file at which to start
; reading. Useful if READFITS is called by another procedure
; which needs to directly read a FITS extension. Should
; always be a multiple of 2880, and not be used with EXTEN_NO
; keyword.
;
; /SILENT - Normally, READFITS will display the size the array at the
; terminal. The SILENT keyword will suppress this
;
; STARTROW - Non-negative integer scalar specifying the row
; of the image or extension table at which to begin reading.
; Useful when one does not want to read the entire table.
;
; HBUFFER - Number of lines in the header, set this to slightly larger
; than the expected number of lines in the FITS header, to
; improve performance when reading very large FITS headers.
; Should be a multiple of 36 -- otherwise it will be modified
; to the next higher multiple of 36. Default is 180
;
; EXAMPLE:
; Read a FITS file test.fits into an IDL image array, IM and FITS
; header array, H. Do not scale the data with BSCALE and BZERO.
;
; IDL> im = READFITS( 'test.fits', h, /NOSCALE)
;
; If the file contains a FITS extension, it could be read with
;
; IDL> tab = READFITS( 'test.fits', htab, /EXTEN )
;
; The function TBGET() can be used for further processing of a binary
; table, and FTGET() for an ASCII table.
; To read only rows 100-149 of the FITS extension,
;
; IDL> tab = READFITS( 'test.fits', htab, /EXTEN,
; STARTR=100, NUMR = 50 )
;
; To read in a file that has been compressed:
;
; IDL> tab = READFITS('test.fits.gz',h)
;
; ERROR HANDLING:
; If an error is encountered reading the FITS file, then
; (1) the system variable !ERROR is set (via the MESSAGE facility)
; (2) the error message is displayed (unless /SILENT is set),
; and the message is also stored in !ERR_STRING
; (3) READFITS returns with a value of -1
; RESTRICTIONS:
; (1) Cannot handle random group FITS
;
; NOTES:
; (1) If data is stored as integer (BITPIX = 16 or 32), and BSCALE
; and/or BZERO keywords are present, then the output array is scaled to
; floating point (unless /NOSCALE is present) using the values of BSCALE
; and BZERO. In the header, the values of BSCALE and BZERO are then
; reset to 1. and 0., while the original values are written into the
; new keywords O_BSCALE and O_BZERO. If the BLANK keyword was
; present, then any input integer values equal to BLANK in the input
; integer image are unchanged by BSCALE or BZERO
;
; (2) The use of the NSLICE keyword is incompatible with the NUMROW
; or STARTROW keywords.
;
; (3) READFITS() underwent a substantial rewrite in October 1998 to
; eliminate recursive calls, and improve efficiency when reading
; extensions.
; 1. The NUMROW and STARTROW keywords can now be used when reading
; a primary image (extension = 0).
; 2. There is no error check for moving past the end of file when
; reading the data array.
;
; (4) On some Unix shells, one may get a "Broken pipe" message if reading
; a UNIX compressed file, and not reading to the end of the file (i.e.
; the decompression has not gone to completion). This is an
; informative message only, and should not affect the output of
; READFITS().
;
; PROCEDURES USED:
; Functions: SXPAR()
; Procedures: SXADDPAR, SXDELPAR
;
; MODIFICATION HISTORY:
; Original Version written in 1988, W.B. Landsman Raytheon STX
; Revision History prior to October 1998 removed
; Major rewrite to eliminate recursive calls when reading extensions
; W.B. Landsman Raytheon STX October 1998
; Add /binary modifier needed for Windows W. Landsman April 1999
; Read unsigned datatypes, added /no_unsigned W. Landsman December 1999
; Output BZERO = 0 for unsigned data types W. Landsman January 2000
; Open with /swap_if_little_endian if since V5.1 W. Landsman February 2000
; Fixed logic error when using NSLICE keyword W. Landsman March 2000
; Fixed byte swapping problem for compressed files on little endian
; machines introduced in Feb 2000 W. Landsman April 2000
; Fix error handling so FREE_LUN is called in case of READU error
; W. Landsman N. Rich, Aug. 2000
; Assume since V5.1, remove NANvalue keyword W. Landsman July 2001
; Support the unofficial 64bit integer format W. Landsman September 2001
; Option to read a unit number rather than file name W.L October 2001
; Fix byte swapping problem for compressed files again (sigh...)
; W. Landsman March 2002
; Make sure gzip defined when using a unit number W. Landsman Oct. 2002
; Assume V5.2 or later W. Landsman Jan. 2003
; Don't read entire header unless needed W. Landsman Jan. 2003
; Added HBUFFER keyword W. Landsman Feb. 2003
; Added CHECKSUM keyword W. Landsman May 2003
; NAME:
; READFMT
; PURPOSE:
; Quickly read a fixed format ASCII data file into IDL variables.
; EXPLANATION:
; Lines of data not meeting the specified format (e.g. comments) are
; ignored.
;
; To read a free format ASCII data file use the procedures
; READCOL or RDFLOAT. To print (formatted or free) columns of data
; use the procedure FORPRINT.
;
; CALLING SEQUENCE:
; READFMT, name, fmt, v1,[ v2, v3, v4, ..., v25 ,
; /SILENT, /DEBUG, SKIPLINE= , NUMLINE =]
;
; INPUTS:
; NAME - Name of ASCII data file. An extension of .DAT is assumed,
; if not supplied.
; FMT - scalar string containing a valid FORTRAN read format.
; Must include a field length specification. Cannot include
; internal parenthesis. A format field must be included for
; each output vector. Multiple format fields are allowed, but
; the repetition factor must be less than 100, (.i.e. 19X is
; allowed but 117X is illegal)
;
; Examples of valid FMT values are
; FMT = 'A7,3X,2I4' or FMT = '1H ,5I7,2A7'
; Examples of INVALID FMT values are
; FMT = 'A7,B3' ;'B' is not a valid FORTRAN format
; FMT = 'A7,2(I3,F5.1)' ;Internal parenthesis not allowed
; FMT = 'A7,F,I' ;Field length not included
;
; OUTPUTS:
; V1,V2,V3,V4... - IDL vectors to contain columns of data.
; Up to 25 output vectors may be read. The type of the output
; vectors are specified by FMT.
;
; OPTIONAL KEYWORD INPUTS:
; /SILENT - If this keyword is set and non-zero, then certain terminal
; output is suppressed while reading the file
; /DEBUG - Set this keyword to display additional information while
; reading the file.
; SKIPLINE - Scalar specifying number of lines to skip at the top of
; file before reading. Default is to start at first line
; NUMLINE - Scalar specifying number of lines in the file to read.
; Default is to read the entire file
;
; EXAMPLES:
; Each row in a fixed-format file POSITION.DAT contains a 5 character
; star name and 6 columns of data giving an RA and Dec in sexigesimal
; format. A possible format for such data might be
;
; IDL> FMT = 'A5,2I3,F5.1,2x,3I3'
; and the file could be quickly read with
;
; IDL> READFMT,'POSITION', fmt, name, hr, min, sec, deg, dmin, dsec
;
; NAME will be a string vector,SEC will be a floating point vector, and
; the other vectors will be of integer type.
;
; RESTRICTIONS:
; This procedure is designed for generality and not for speed.
; If a large ASCII file is to be read repeatedly, it may be worth
; writing a specialized reader.
;
; NOTES:
; When reading a field with an integer format I, the output vector is
; byte - if n = 1
; integer*2 - if 1 < n < 5
; integer*4 - in all other cases
; Octal ('O') and hexadecimal ('Z') formats are read into longwords
;
; PROCEDURE CALLS:
; GETTOK(), NUMLINES(), REMCHAR, ZPARCHECK
;
; REVISION HISTORY:
; Written W. Landsman November, 1988
; Added SKIPLINE and NUMLINE keywords March 92
; Allow up to 25 columns to be read June 92
; Call NUMLINES() function Feb 1996
; Converted to IDL V5.0 W. Landsman September 1997
; Recognize 'O' and 'Z' formats W. Landsman September 1997
; NAME:
; SUNPOS
; PURPOSE:
; To compute the RA and Dec of the Sun at a given date.
;
; CALLING SEQUENCE:
; SUNPOS, jd, ra, dec, [elong, obliquity, /RADIAN ]
; INPUTS:
; jd - The Julian date of the day (and time), scalar or vector
; usually double precision
; OUTPUTS:
; ra - The right ascension of the sun at that date in DEGREES
; double precision, same number of elements as jd
; dec - The declination of the sun at that date in DEGREES
;
; OPTIONAL OUTPUTS:
; elong - Ecliptic longitude of the sun at that date in DEGREES.
; obliquity - the obliquity of the ecliptic, in DEGREES
;
; OPTIONAL INPUT KEYWORD:
; /RADIAN - If this keyword is set and non-zero, then all output variables
; are given in Radians rather than Degrees
;
; NOTES:
; The accuracy in the 20th century should be within 1"; however this
; has not been extensively tested.
;
; The returned RA and Dec are in the given date's equinox.
;
; Procedure was extensively revised in May 1996, and the new calling
; sequence is incompatible with the old one.
; METHOD:
; Uses a truncated version of Newcomb's Sun. Adapted from the IDL
; routine SUN_POS by CD Pike, which was adapted from a FORTRAN routine
; by B. Emerson (RGO).
; EXAMPLE:
; (1) Find the apparent RA and Dec of the Sun on May 1, 1982
;
; IDL> jdcnv, 1982, 5, 1,0 ,jd ;Find Julian date jd = 2445090.5
; IDL> sunpos, jd, ra, dec
; IDL> print,adstring(ra,dec,2)
; 02 31 32.61 +14 54 34.9
;
; The Astronomical Almanac gives 02 31 32.58 +14 54 34.9 so the error
; in SUNPOS for this case is < 0.5".
;
; (2) Find the apparent RA and Dec of the Sun for every day in 1997
;
; IDL> jdcnv, 1997,1,1,0, jd ;Julian date on Jan 1, 1997
; IDL> sunpos, jd+ dindgen(365), ra, dec ;RA and Dec for each day
;
; MODIFICATION HISTORY:
; Written by Michael R. Greason, STX, 28 October 1988.
; Accept vector arguments, W. Landsman April,1989
; Eliminated negative right ascensions. MRG, Hughes STX, 6 May 1992.
; Rewritten using the 1993 Almanac. Keywords added. MRG, HSTX,
; 10 February 1994.
; Major rewrite, improved accuracy, always return values in degrees
; W. Landsman May, 1996
; Added /RADIAN keyword, W. Landsman August, 1997
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; SXADDPAR
; PURPOSE:
; Add or modify a parameter in a FITS header array.
;
; CALLING SEQUENCE:
; SXADDPAR, Header, Name, Value, [ Comment, Location, /SaveComment,
; BEFORE =, AFTER = , FORMAT= , /PDU]
;
; INPUTS:
; Header = String array containing FITS or STSDAS header. The
; length of each element must be 80 characters. If not
; defined, then SXADDPAR will create an empty FITS header array.
;
; Name = Name of parameter. If Name is already in the header the value
; and possibly comment fields are modified. Otherwise a new
; record is added to the header. If name is equal to 'COMMENT'
; or 'HISTORY' or a blank string then the value will be added to
; the record without replacement. For these cases, the comment
; parameter is ignored.
;
; Value = Value for parameter. The value expression must be of the
; correct type, e.g. integer, floating or string. String values
; of 'T' or 'F' are considered logical values.
;
; OPTIONAL INPUT PARAMETERS:
; Comment = String field. The '/' is added by this routine. Added
; starting in position 31. If not supplied, or set equal to
; '', or /SAVECOMMENT is set, then the previous comment field is
; retained (when found)
;
; Location = Keyword string name. The parameter will be placed before the
; location of this keyword. This parameter is identical to
; the BEFORE keyword and is kept only for consistency with
; earlier versions of SXADDPAR.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
; BEFORE = Keyword string name. The parameter will be placed before the
; location of this keyword. For example, if BEFORE='HISTORY'
; then the parameter will be placed before the first history
; location. This applies only when adding a new keyword;
; keywords already in the header are kept in the same position.
;
; AFTER = Same as BEFORE, but the parameter will be placed after the
; location of this keyword. This keyword takes precedence over
; BEFORE.
;
; FORMAT = Specifies FORTRAN-like format for parameter, e.g. "F7.3". A
; scalar string should be used. For complex numbers the format
; should be defined so that it can be applied separately to the
; real and imaginary parts. If not supplied then the default is
; 'G19.12' for double precision, and 'G14.7' for floating point.
;
; /PDU = specifies keyword is to be added to the primary data unit
; header. If it already exists, it's current value is updated in
; the current position and it is not moved.
; /SAVECOMMENT = if set, then any existing comment is retained, i.e. the
; COMMENT parameter only has effect if the keyword did not
; previously exist in the header.
; OUTPUTS:
; Header = updated FITS header array.
;
; EXAMPLE:
; Add a keyword 'TELESCOP' with the value 'KPNO-4m' and comment 'Name
; of Telescope' to an existing FITS header h.
;
; IDL> sxaddpar, h, 'TELESCOPE','KPNO-4m','Name of Telescope'
; NOTES:
; The functions SXADDPAR() and FXADDPAR() are nearly identical, with the
; major difference being that FXADDPAR forces required FITS keywords
; BITPIX, NAXISi, EXTEND, PCOUNT, GCOUNT to appear in the required order
; in the header, and FXADDPAR supports the OGIP LongString convention.
; There is no particular reason for having two nearly identical
; procedures, but both are too widely used to drop either one.
;
; All HISTORY records are inserted in order at the end of the header.
;
; All COMMENT records are also inserted in order at the end of the header
; header, but before the HISTORY records. The BEFORE and AFTER keywords
; can override this.
;
; All records with no keyword (blank) are inserted in order at the end of
; the header, but before the COMMENT and HISTORY records. The BEFORE and
; AFTER keywords can override this.
; RESTRICTIONS:
; Warning -- Parameters and names are not checked
; against valid FITS parameter names, values and types.
;
; MODIFICATION HISTORY:
; DMS, RSI, July, 1983.
; D. Lindler Oct. 86 Added longer string value capability
; Converted to NEWIDL D. Lindler April 90
; Added Format keyword, J. Isensee, July, 1990
; Added keywords BEFORE and AFTER. K. Venkatakrishna, May '92
; Pad string values to at least 8 characters W. Landsman April 94
; Aug 95: added /PDU option and changed routine to update last occurence
; of an existing keyword (the one SXPAR reads) instead of the
; first occurence.
; Comment for string data can start after column 32 W. Landsman June 97
; Make sure closing quote supplied with string value W. Landsman June 98
; Converted to IDL V5.0 W. Landsman June 98
; Increase precision of default formatting of double precision floating
; point values. C. Gehman, JPL September 1998
; Mar 2000, D. Lindler, Modified to use capital E instead of lower case
; e for exponential formats.
; Apr 2000, Make user-supplied format upper-case W. Landsman
; Oct 2001, Treat COMMENT or blank string like HISTORY keyword W. Landsman
; Jan 2002, Allow BEFORE, AFTER to apply to COMMENT keywords W. Landsman
; June 2003, Added SAVECOMMENT keyword W. Landsman
;
; NAME:
; SXPAR
; PURPOSE:
; Obtain the value of a parameter in a FITS header
;
; CALLING SEQUENCE:
; result = SXPAR( Hdr, Name, [ Abort, COUNT=, COMMENT =, /NoCONTINUE ])
;
; INPUTS:
; Hdr = FITS header array, (e.g. as returned by READFITS)
; string array, each element should have a length of 80 characters
;
; Name = String name of the parameter to return. If Name is of the
; form 'keyword*' then an array is returned containing values of
; keywordN where N is an integer. The value of keywordN will be
; placed in RESULT(N-1). The data type of RESULT will be the
; type of the first valid match of keywordN found.
;
; OPTIONAL INPUTS:
; ABORT - string specifying that SXPAR should do a RETALL
; if a parameter is not found. ABORT should contain
; a string to be printed if the keyword parameter is not found.
; If not supplied, SXPAR will return quietly with COUNT = 0
; (and !ERR = -1) if a keyword is not found.
;
; OPTIONAL INPUT KEYWORDS:
; /NOCONTINUE = If set, then continuation lines will not be read, even
; if present in the header
; /SILENT - Set this keyword to suppress warning messages about duplicate
; keywords in the FITS header.
;
; OPTIONAL OUTPUT KEYWORDS:
; COUNT - Optional keyword to return a value equal to the number of
; parameters found by SXPAR, integer scalar
;
; COMMENT - Array of comments associated with the returned values
;
; OUTPUTS:
; Function value = value of parameter in header.
; If parameter is double precision, floating, long or string,
; the result is of that type. Apostrophes are stripped
; from strings. If the parameter is logical, 1b is
; returned for T, and 0b is returned for F.
; If Name was of form 'keyword*' then a vector of values
; are returned.
;
; SIDE EFFECTS:
; !ERR is set to -1 if parameter not found, 0 for a scalar
; value returned. If a vector is returned it is set to the
; number of keyword matches found. The use of !ERR is deprecated, and
; instead the COUNT keyword is preferred
;
; If a keyword (except HISTORY or COMMENT) occurs more than once in a
; header, a warning is given, and the *last* occurence is used.
;
; EXAMPLES:
; Given a FITS header, h, return the values of all the NAXISi values
; into a vector. Then place the history records into a string vector.
;
; IDL> naxisi = sxpar( h ,'NAXIS*') ; Extract NAXISi value
; IDL> history = sxpar( h, 'HISTORY' ) ; Extract HISTORY records
;
; PROCEDURE:
; The first 8 chacters of each element of Hdr are searched for a
; match to Name. The value from the last 20 characters is returned.
; An error occurs if there is no parameter with the given name.
;
; If a numeric value has no decimal point it is returned as type
; LONG. If it contains more than 8 numerals, or contains the
; characters 'D' or 'E', then it is returned as type DOUBLE. Otherwise
; it is returned as type FLOAT. Very large integer values, outside
; the range of valid LONG, are returned as DOUBLE.
;
; If the value is too long for one line, it may be continued on to the
; the next input card, using the OGIP CONTINUE convention. For more info,
; http://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/ofwg_recomm/r13.html
;
; Complex numbers are recognized as two numbers separated by one or more
; space characters.
;
; If a numeric value has no decimal point (or E or D) it is returned as
; type LONG. If it contains more than 8 numerals, or contains the
; character 'D', then it is returned as type DOUBLE. Otherwise it is
; returned as type FLOAT. If an integer is too large to be stored as
; type LONG, then it is returned as DOUBLE.
;
; NOTES:
; The functions SXPAR() and FXPAR() are nearly identical, although
; FXPAR() has slightly more sophisticated parsing. There is no
; particular reason for having two nearly identical procedures, but
; both are too widely used to drop either one.
;
; PROCEDURES CALLED:
; GETTOK(), VALID_NUM()
; MODIFICATION HISTORY:
; DMS, May, 1983, STPAR Written.
; D. Lindler Jan 90 added ABORT input parameter
; J. Isensee Jul,90 added COUNT keyword
; W. Thompson, Feb. 1992, added support for FITS complex values.
; W. Thompson, May 1992, corrected problem with HISTORY/COMMENT/blank
; keywords, and complex value error correction.
; W. Landsman, November 1994, fix case where NAME is an empty string
; W. Landsman, March 1995, Added COMMENT keyword, ability to read
; values longer than 20 character
; W. Landsman, July 1995, Removed /NOZERO from MAKE_ARRAY call
; T. Beck May 1998, Return logical as type BYTE
; W. Landsman May 1998, Make sure integer values are within range of LONG
; Converted to IDL V5.0, May 1998
; W. Landsman Feb 1998, Recognize CONTINUE convention
; W. Landsman Oct 1999, Recognize numbers such as 1E-10 as floating point
; W. Landsman Jan 2000, Only accept integer N values when name = keywordN
; W. Landsman Dec 2001, Optional /SILENT keyword to suppress warnings
; W. Landsman/D. Finkbeiner Mar 2002 Make sure extracted vectors
; of mixed data type are returned with the highest type.
; NAME:
; TDB2TDT
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
; UPDATED VERSIONs can be found on my WEB PAGE:
; http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
; Relativistic clock corrections due to Earth motion in solar system
;
; MAJOR TOPICS:
; Planetary Orbits
;
; CALLING SEQUENCE:
; corr = TDB2TDT(JD, TBASE=, DERIV=deriv)
;
; DESCRIPTION:
;
; The function TDB2TDT computes relativistic corrections that must
; be applied when performing high precision absolute timing in the
; solar system.
;
; According to general relativity, moving clocks, and clocks at
; different gravitational potentials, will run at different rates
; with respect to each other. A clock placed on the earth will run
; at a time-variable rate because of the non-constant influence of
; the sun and other planets. Thus, for the most demanding
; astrophysical timing applications -- high precision pulsar timing
; -- times in the accelerating earth observer's frame must be
; corrected to an inertial frame, such as the solar system
; barycenter (SSB). This correction is also convenient because the
; coordinate time at the SSB is the ephemeris time of the JPL
; Planetary Ephemeris.
;
; In general, the difference in the rate of Ti, the time kept by an
; arbitrary clock, and the rate of T, the ephemeris time, is given
; by the expression (Standish 1998):
;
; dTi/dT = 1 - (Ui + vi^2/2) / c^2
;
; where Ui is the potential of clock i, and vi is the velocity of
; clock i. However, when integrated, this expression depends on the
; position of an individual clock. A more convenient approximate
; expression is:
;
; T = Ti + (robs(Ti) . vearth(T))/c^2 + dtgeo(Ti) + TDB2TDT(Ti)
;
; where robs is the vector from the geocenter to the observer;
; vearth is the vector velocity of the earth; and dtgeo is a
; correction to convert from the observer's clock to geocentric TT
; time. TDB2TDT is the value computed by this function, the
; correction to convert from the geocenter to the solar system
; barycenter.
;
; As the above equation shows, while this function provides an
; important component of the correction, the user must also be
; responsible for (a) correcting their times to the geocenter (ie,
; by maintaining atomic clock corrections); (b) estimating the
; observatory position vector; and and (c) estimating earth's
; velocity vector (using JPLEPHINTERP).
;
; Users may note a circularity to the above equation, since
; vearth(T) is expressed in terms of the SSB coordinate time. This
; appears to be a chicken and egg problem since in order to get the
; earth's velocity, the ephemeris time is needed to begin with.
; However, to the precision of the above equation, < 25 ns, it is
; acceptable to replace vearth(T) with vearth(TT).
;
; The method of computation of TDB2TDT in this function is based on
; the analytical formulation by Fairhead, Bretagnon & Lestrade, 1988
; (so-called FBL model) and Fairhead & Bretagnon 1990, in terms of
; sinusoids of various amplitudes. TDB2TDT has a dominant periodic
; component of period 1 year and amplitude 1.7 ms. The set of 791
; coefficients used here were drawn from the Princeton pulsar timing
; program TEMPO version 11.005 (Taylor & Weisberg 1989).
;
; Because the TDB2TDT quantity is rather expensive to compute but
; slowly varying, users may wish to also retrieve the time
; derivative using the DERIV keyword, if they have many times to
; convert over a short baseline.
;
; Verification
;
; This implementation has been compared against a set of FBL test
; data found in the 1996 IERS Conventions, Chapter 11, provided by
; T. Fukushima. It has been verified that this routine reproduces
; the Fukushima numbers to the accuracy of the table, within
; 10^{-14} seconds.
;
; Fukushima (1995) has found that the 791-term Fairhead & Bretagnon
; analytical approximation use here has a maximum error of 23
; nanoseconds in the time range 1980-2000, compared to a numerical
; integration. In comparison the truncated 127-term approximation
; has an error of ~130 nanoseconds.
;
;
; PARAMETERS:
;
; JD - Geocentric time TT, scalar or vector, expressed in Julian
; days. The actual time used is (JD + TBASE). For maximum
; precision, TBASE should be used to express a fixed epoch in
; whole day numbers, and JD should express fractional offset
; days from that epoch.
;
;
; KEYWORD PARAMETERS:
;
; TBASE - scalar Julian day of a fixed epoch, which provides the
; origin for times passed in JD.
; Default: 0
;
; DERIV - upon return, contains the derivative of TDB2TDT in units
; of seconds per day. As many derivatives are returned as
; values passed in JD.
;
;
; RETURNS:
; The correction offset(s) in units of seconds, to be applied as
; noted above.
;
;
; EXAMPLE:
;
; Find the correction at ephemeris time 2451544.5 (JD):
; IDL> print, tdb2tdt(2451544.5d)
; -0.00011376314
; or 0.11 ms.
;
;
; REFERENCES:
;
; Princeton TEMPO Program
; http://pulsar.princeton.edu/tempo/
;
; FBL Test Data Set
; ftp://maia.usno.navy.mil/conventions/chapter11/fbl.results
;
; Fairhead, L. & Bretagnon, P. 1990, A&A, 229, 240
; (basis of this routine)
;
; Fairhead, L. Bretagnon, P. & Lestrade, J.-F. 1988, in *The Earth's
; Rotation and Reference Frames for Geodesy and Geodynamics*,
; ed. A. K. Babcock and G. A. Wilkins, (Dordrecht: Kluwer), p. 419
; (original "FBL" paper)
;
; Fukushima, T. 1995, A&A, 294, 895 (error analysis)
;
; Irwin, A. W. & Fukushima, T. 1999, A&A, 348, 642 (error analysis)
;
; Standish, E. M. 1998, A&A, 336, 381 (description of time scales)
;
; Taylor, J. H. & Weisberg, J. M. 1989, ApJ, 345, 434 (pulsar timing)
;
;
; SEE ALSO
; JPLEPHREAD, JPLEPHINTERP, JPLEPHTEST
;
; MODIFICATION HISTORY:
; Original logic from Fairhead & Bretagnon, 1990
; Drawn from TEMPO v. 11.005, copied 20 Jun 2001
; Documented and vectorized, 30 Jun 2001
;
;
; $Id: tdb2tdt.pro,v 1.4 2001/07/01 07:37:40 craigm Exp $
;
; NAME:
; TSC
;
; PURPOSE:
; Interpolate an irregularly sampled field using a Triangular Shaped Cloud
;
; EXPLANATION:
; This function interpolates an irregularly sampled field to a
; regular grid using Triangular Shaped Cloud (nearest grid point
; gets weight 0.75-dx^2, points before and after nearest grid
; points get weight 0.5*(1.5-dx)^2, where dx is the distance
; from the sample to the grid point in units of the cell size).
;
; CATEGORY:
; Mathematical functions, Interpolation
;
; CALLING SEQUENCE:
; Result = TSC, VALUE, POSX, NX[, POSY, NY, POSZ, NZ,
; AVERAGE = average, WRAPAROUND = wraparound,
; ISOLATED = isolated, NO_MESSAGE = no_message]
;
; INPUTS:
; VALUE: Array of sample weights (field values). For e.g. a
; temperature field this would be the temperature and the
; keyword AVERAGE should be set. For e.g. a density field
; this could be either the particle mass (AVERAGE should
; not be set) or the density (AVERAGE should be set).
; POSX: Array of X coordinates of field samples, unit indices: [0,NX>.
; NX: Desired number of grid points in X-direction.
;
; OPTIONAL INPUTS:
; POSY: Array of Y coordinates of field samples, unit indices: [0,NY>.
; NY: Desired number of grid points in Y-direction.
; POSZ: Array of Z coordinates of field samples, unit indices: [0,NZ>.
; NZ: Desired number of grid points in Z-direction.
;
; KEYWORD PARAMETERS:
; AVERAGE: Set this keyword if the nodes contain field samples
; (e.g. a temperature field). The value at each grid
; point will then be the weighted average of all the
; samples allocated to it. If this keyword is not
; set, the value at each grid point will be the
; weighted sum of all the nodes allocated to it
; (e.g. for a density field from a distribution of
; particles). (D=0).
; WRAPAROUND: Set this keyword if you want the first grid point
; to contain samples of both sides of the volume
; (see below).
; ISOLATED: Set this keyword if the data is isolated, i.e. not
; periodic. In that case total `mass' is not conserved.
; This keyword cannot be used in combination with the
; keyword WRAPAROUND.
; NO_MESSAGE: Suppress informational messages.
;
; Example of default allocation of nearest grid points: n0=4, *=gridpoint.
;
; 0 1 2 3 Index of gridpoints
; * * * * Grid points
; |---|---|---|---| Range allocated to gridpoints ([0.0,1.0> --> 0, etc.)
; 0 1 2 3 4 posx
;
; Example of ngp allocation for WRAPAROUND: n0=4, *=gridpoint.
;
; 0 1 2 3 Index of gridpoints
; * * * * Grid points
; |---|---|---|---|-- Range allocated to gridpoints ([0.5,1.5> --> 1, etc.)
; 0 1 2 3 4=0 posx
;
;
; OUTPUTS:
; Prints that a TSC interpolation is being performed of x
; samples to y grid points, unless NO_MESSAGE is set.
;
; RESTRICTIONS:
; Field data is assumed to be periodic with the sampled volume
; the basic cell, unless ISOLATED is set.
; All input arrays must have the same dimensions.
; Postition coordinates should be in `index units' of the
; desired grid: POSX=[0,NX>, etc.
; Keywords ISOLATED and WRAPAROUND cannot both be set.
;
; PROCEDURE:
; Nearest grid point is determined for each sample.
; TSC weights are computed for each sample.
; Samples are interpolated to the grid.
; Grid point values are computed (sum or average of samples).
;
; EXAMPLE:
; nx=20
; ny=10
; posx=randomu(s,1000)
; posy=randomu(s,1000)
; value=posx^2+posy^2
; field=tsc(value,posx*nx,nx,posy*ny,ny,/average)
; surface,field,/lego
;
; NOTES:
; Use csc.pro or ngp.pro for lower order interpolation schemes. A
; standard reference for these interpolation methods is: R.W. Hockney
; and J.W. Eastwood, Computer Simulations Using Particles (New York:
; McGraw-Hill, 1981).
;
; MODIFICATION HISTORY:
; Written by Joop Schaye, Feb 1999.
; Check for overflow for large dimensions P. Riley/W. Landsman Dec. 1999
; NAME:
; TVCIRCLE
; PURPOSE:
; Draw circle(s) of specified radius at specified position(s)
; EXPLANATION:
; If a position is not specified, and device has a cursor, then a circle
; is drawn at the current cursor position.
;
; CALLING SEQUENCE:
; TVCIRCLE, rad, x, y, color, [ /DATA, /FILL, _EXTRA = ]
;
; INPUTS:
; RAD - radius of circle(s) to be drawn, positive numeric scalar
;
; OPTIONAL INPUT:
; X - x position for circle center, vector or scalar
; Y - y position for circle center, vector or scalar
; If X and Y are not specified, and the device has a cursor,
; then program will draw a circle at the current cursor position
; COLOR - intensity value(s) (0 - !D.N_COLORS) used to draw the circle(s)
; If COLOR is a scalar then all circles are drawn with the same
; color value. Otherwise, the Nth circle is drawn with the
; Nth value of color. Default = !P.COLOR.
;
; OPTONAL KEYWORD INPUTS:
; /DATA - if this keyword is set and non-zero, then the circle width and
; X,Y position center are interpreted as being in DATA
; coordinates. Note that data coordinates must be previously
; defined (with a PLOT or CONTOUR call). TVCIRCLE will
; internally convert to device coordinates before drawing the
; circle, in order to maintain optimal smoothness.
; /FILL - If set, fill the circle using POLYFILL
;
; Any keyword recognized by PLOTS (or POLYFILL if /FILL is set)
; is also recognized by TVCIRCLE. In particular, the color,
; linestyle, and thickness of the circles are controlled by the
; COLOR, LINESTYLE, and THICK keywords. If POLYFILL is set
; then available keywords are LINE_FILL and FILL_PATTERN.
; OUTPUTS:
; None
;
; RESTRICTIONS:
; (1) TVCIRCLE does not check whether it writes off of the edge of the
; display
; (2) Some round-off error may occur when non-integral values are
; supplied for both the radius and the center coordinates
; (3) TVCIRCLE does not accept /NORMAL coordinates, only data coordinates
; (if /DATA is set) or device coordinates (the default)
; (4) TVCIRCLE always draws a circle --- even if /DATA is set, and the
; X and Y data scales are unequal. (The X data scale is used to
; define the circle radius.) If this is not the behaviour
; you want, then use TVELLIPSE instead.
; EXAMPLE:
; (1) Draw circles of radius 9 pixels at the positions specified by
; X,Y vectors, using double thickness lines
;
; IDL> tvcircle, 9, x, y, THICK = 2
;
; Now fill in the circles using the LINE_FILL method
;
; IDL> tvcircle, 9, x, y, /FILL, /LINE_FILL
; METHOD:
; The method used is that of Michener's, modified to take into account
; the fact that IDL plots arrays faster than single points. See
; "Fundamental of Interactive Computer Graphics" by Foley and Van Dam"
; p. 445 for the algorithm.
;
; REVISON HISTORY:
; Original version written by B. Pfarr STX 10-88
; Major rewrite adapted from CIRCLE by Allyn Saroyan LNLL
; Wayne Landsman STX Sep. 91
; Added DATA keyword Wayne Landsman HSTX June 1993
; Added FILL keyword. R. S. Hill, HSTX, 4-Nov-1993
; Always convert to device coords, add _EXTRA keyword, allow vector
; colors. Wayne Landsman, HSTX, May 1995
; Allow one to set COLOR = 0, W. Landsman, HSTX, November 1995
; Check if data axes reversed. P. Mangifico, W. Landsman May 1996
; Converted to IDL V5.0 W. Landsman September 1997
; NAME:
; TVLASER
; PURPOSE:
; Prints screen or image array onto a Postscript file or printer.
; Information from FITS header is optionally used for labeling.
;
; CALLING SEQUENCE:
; TVLASER, [header, Image, BARPOS = ,CARROWS =, CLABELS = ,/COLORPS,
; COMMENTS = ,CSIZE = ,CTITLE = , DX = , DY =, /ENCAP, FILENAME =
; HEADER = ,/HELP, IMAGEOUT = ,/INTERP, /MAGNIFY, /NoCLOSE,
; /NoDELETE, /NO_PERS_INFO, /NoEIGHT, /NoPRINT, /NoRETAIN,
; /PORTRAIT, PRINTER = , /REVERSE, /SCALE, TITLE = , /TrueColor,
; XDIM=, XSTART=, YDIM=, YSTART=, BOTTOMDW=, NCOLORSDW= ]
;
; Note that the calling sequence was changed in May 1997
; OPTIONAL INPUTS:
; HEADER - FITS header string array. Object and astrometric info from
; the FITS header will be used for labeling, if available
; IMAGE - if an array is passed through this parameter, then this image
; will be used rather than reading off the current window. This
; allows easy use of large images. It is usually preferable
; to optimally byte scale IMAGE before supplying it to TVLASER
;
; OPTIONAL KEYWORD INPUT PARAMETERS:
; BARPOS - A four- or five-element vector giving the position and
; orientation of the color bar. The first four elements
; [X0,Y0,XSize,YSize] indicate the position and size of the color
; bar in INCHES, relative to origin of the displayed image.
; (X0,Y0) are the position of the lower left corner and
; (XSize,YSize) are the width and height. The fifth element is
; optional, and if present, the color bar will be printed
; horizontally rather than vertically. If BARPOS is set to
; anything but a four- or five-element vector, the bar is NOT
; printed. The default value is BARPOS = [-0.25, 0.0, 0.2, 2.0]
; BOTTOMDW - The lowest value to use in building the density
; wedge. Used with NCOLORSDW. Compatible with BOTTOM and
; NCOLORS keywords of XLOADCT.
; CARROWS - The color to print the North-East arrows. Default is dark.
; Three types of values can be passed:
; SCALAR: that value's color in the current color table
; 3-ELEMENT VECTOR: the color will be [R,G,B]
; STRING: A letter indicating the color. Valid names are:
; 'W' (white), 'D' (dark/black), 'R' (red), 'G' (green),
; 'B' (blue), 'T' (turquoise), 'V' (violet), 'Y' (yellow),
; If the keyword is set to a value of -1, the arrows are
; NOT printed.
; COLORPS - If present and non-zero, the idl.ps file is written using
; color postscript.
; COMMENTS - A string that will be included in the comment line below the
; image. For multi-line comments you can either use "!C" in the
; string as a carriage return {although the vertical spacing
; might be a little off} or, preferably, make the COMMENTS a
; string array with each line as a separate element.
; CLABELS - Color to print the labels, same format as for CARROWS.
; CSIZE - Color to print the size-scale bar and label, same format as for
; CARROWS.
; CTITLE - Color to print the title, same format as for CARROWS.
; DX,DY - offsets in INCHES added to the position of the figure on the
; paper. As is the case for the device keywords XOFFSET and
; YOFFSET, when in landscape mode DX and DY are the same
; *relative to the paper*, not relative to the plot (e.g., DX is
; the horizontal offset in portrait mode, but the *vertical*
; offset in landscape mode).
; ENCAP - If present and non-zero, the IDL.PS file is written in
; encapsulated postscript for import into LaTeX documents
; FILENAME - scalar string giving name of output postscript file.
; Default is idl.ps. Automatically sets /NODELETE
; HEADER = FITS header. This is an alternative to supplying the FITS
; header in the first parameter.
; HELP - print out the sytax for this procedure.
; INTERP - If present and non-zero, current color table will be
; interpolated to fill the full range of the PostScript color
; table (256 colors). Otherwise, the current color table will be
; directly copied. You probably will want to use this if you
; are using IMAGE keyword and a shared color table.
; MAGNIFY - The net magnification of the entire figure. At this point,
; the figure is not automatically centered on the paper if the
; value of MAGNIFY is not equal to 1, but the DX and DY keywords
; can be used to shift location. For example, to fit a full plot
; on the printable area (8.5x8.5 inches) of the Tek PhaserIISD
; color printer use: MAGNIFY=0.8, DX=0.5, DY=0.5.;
; NCOLORSDW - The number of values to include in the density
; wedge. Used with BOTTOMDW. Compatible with
; BOTTOM/NCOLORS keywords of XLOADCT.
; NoCLOSE - If present and non-zero, then the postscript file is not
; closed (or printed), the device is set to 'PS', and the data
; coordinate system is set to match the image size. This allows the
; user to add additional plotting commands before printing. For
; example, to include a 15 pixel circle around a source at
; coordinates (150,160), around an image, im, with FITS header
; array, h
;
; IDL> tvlaser,h,im,/NoClose ;Write image & annotation
; IDL> tvcircle,15,150,160,/data ;Draw circle
; IDL> device,/close ;Close postscript file & print
;
; NoDELETE - If present and non-zero, the postscript file is kept AND is
; also sent to the printer
; NoEIGHT - if set then only four bits sent to printer (saves space)
; NO_PERS_INFO - if present and non-zero, output notation will NOT
; include date/user block of information.
; NoPRINT - If present and non-zero, the output is sent to a file (default
; name 'idl.ps'), which is NOT deleted and is NOT sent to the
; printer.
; NoRETAIN - In order to avoid possible problems when using TVRD with
; an obscured window, TVLASER will first copy the current window
; to a temporary RETAIN=2 window. Set /NORETAIN to skip this
; step and improve performance
; PORTRAIT - if present and non-zero, the printer results will be in
; portrait format; otherwise, they will be in landscape format.
; If labels are requested, image will be in portrait mode,
; regardless
; PRINTER - scalar string giving the OS command to send a the postscript
; file to the printer. Under Unix, the default value of PRINTER
; is 'lpr ' while for other OS it is 'print '
; REVERSE - if present and non-zero, color table will be fliped, so black
; and white are reversed.
; SCALE - if present and non-zero, image will be bytscaled before being
; sent to postscript file.
; TITLE - if present and non-zero, the string entered here will be the
; title of the picture. Default is the OBJECT field in the
; header (if present).
; TRUECOLOR - if present and non-zero, the postscript file is created
; using the truecolor switch (i.e. true=3). The colorbar is
; not displayed in this mode.
; XDIM,YDIM - Number of pixels. Default is from !d.x_size and !d.y_size,
; or size of image if passed with IMAGE keyword.
; XSTART,YSTART - lower left corner (default of (0,0))
;
; OPTIONAL KEYWORD OUTPUT PARAMETER
; IMAGEOUT = the image byte array actually sent to the postscript file.
;
; SIDE EFFECTS:
; A postscript file is created in the current directory. User must have
; write privileges in the current directory. The file is named idl.ps
; unless the FILENAME keyword is given. The file is directed to the
; printer unless the /ENCAP, /NoCLOSE, or /NOPRINT keywords are given.
; After printing, the file is deleted unless the /NODELETE or FILENAME
; keywords are given.
; PROCEDURE:
; Read display or take IMAGE and then redisplay into a postscript file.
; If a header exists, printout header information. If header has
; astrometry, then print out orientation and scale information.
; PROCEDURES USED:
; ARROWS, EXTAST, FDECOMP, GETROT, PIXCOLOR, SXPAR(), XYAD, ZPARCHECK
;
;*EXAMPLE:
; 1) Send a true color image (xsize,ysize,3) to a printer (i.e. print23l),
; tvlaser,huv,cpic,/colorps,/truecolor,printer="print23l"
; % TVLASER: Now printing image: $print23l idl.ps
;
; MODIFICATION HISTORY:
; Major rewrite from UIT version W. Landsman Dec 94
; Massive rewrite. Added North-East arrows, pixel scale bar, color bar,
; and keywords DX, DY, MAGNIFY, INTERP, HELP, and COMMENTS.
; Created ablility to define colors for annotation and
; text. Repositioned text labels. J.Wm.Parker, HITC, 5/95
; Make Header and Image parameters instead of keywords. Add PRINTER
; keyword. Include alternate FITS keywords. W. Landsman May 97
; Copy to a RETAIN=2 window, work without FITS header W. Landsman June 97
; Cleaner output when no astrometry in header W. Landsman June 97
; Added /INFO to final MESSAGE W. Landsman July 1997
; 12/4/97 jkf/acc - added TrueColor optional keyword.
; Added /NoClose keyword, trim Equinox format W. Landsman 9-Jul-1998
; Don't display coordinate labels if no astrometry, more flexible
; formatting of exposure time W. Landsman 30-Aug-1998
; BottomDW and NColorsDW added. R. S. Hill, 1-Mar-1999
; Apply func tab to color bar if not colorps. RSH, 21 Mar 2000
; Fix problem with /NOCLOSE and unequal X,Y sizes W. Landsman Feb 2001
; Use TVRD(True=3) if /TRUECOLOR set W. Landsman November 2001
; NAME:
; UVBYBETA
; PURPOSE:
; Derive dereddened colors, metallicity, and Teff from Stromgren colors.
; EXPLANATION:
; Adapted from FORTRAN routine of same name published by T.T. Moon,
; Communications of University of London Observatory, No. 78. Parameters
; can either be input interactively (with /PROMPT keyword) or supplied
; directly.
;
; CALLING SEQUENCE:
; uvbybeta, /PROMPT ;Prompt for all parameters
; uvbybeta,by,m1,c1,Hbeta,n ;Supply inputs, print outputs
; uvbybeta, by, m1, c1, Hbeta, n, Te, Mv, Eby, delm0, radius,
; [ TEXTOUT=, Eby_in =, Name = ]
;
; INPUTS:
; by - Stromgren b-y color, scalar or vector
; m1 - Stromgren line-blanketing parameter, scalar or vector
; c1 - Stromgren Balmer discontinuity parameter, scalar or vector
; Hbeta - H-beta line strength index. Set Hbeta to 0 if it is not
; known, and UVBYBETA will estimate a value based on by, m1,and c1.
; Hbeta is not used for stars in group 8.
; n - Integer (1-8), scalar or vector, giving approximate stellar
; classification
;
; (1) B0 - A0, classes III - V, 2.59 < Hbeta < 2.88,-0.20 < c0 < 1.00
; (2) B0 - A0, class Ia , 2.52 < Hbeta < 2.59,-0.15 < c0 < 0.40
; (3) B0 - A0, class Ib , 2.56 < Hbeta < 2.61,-0.10 < c0 < 0.50
; (4) B0 - A0, class II , 2.58 < Hbeta < 2.63,-0.10 < c0 < 0.10
; (5) A0 - A3, classes III - V, 2.87 < Hbeta < 2.93,-0.01 < (b-y)o< 0.06
; (6) A3 - F0, classes III - V, 2.72 < Hbeta < 2.88, 0.05 < (b-y)o< 0.22
; (7) F1 - G2, classes III - V, 2.60 < Hbeta < 2.72, 0.22 < (b-y)o< 0.39
; (8) G2 - M2, classes IV _ V, 0.20 < m0 < 0.76, 0.39 < (b-y)o< 1.00
;
;
; OPTIONAL INPUT KEYWORD:
; Eby_in - numeric scalar specifying E(b-y) color to use. If not
; supplied, then E(b-y) will be estimated from the Stromgren colors
; NAME - scalar or vector string giving name(s) of star(s). Used only
; when writing to disk for identification purposes.
; /PROMPT - if set, then uvbybeta.pro will prompt for Stromgren indicies
; interactively
; TEXTOUT - Used to determine output device. If not present, the
; value of the !TEXTOUT system variable is used (see TEXTOPEN)
; textout=1 Terminal with /MORE (if a tty)
; textout=2 Terminal without /MORE
; textout=3 uvbybeta.prt (output file)
; textout=4 Laser Printer
; textout=5 User must open file
; textout=7 Append to existing uvbybeta.prt file
; textout = filename (default extension of .prt)
; /PRINT - if set, then force display output information to the device
; specified by !TEXTOUT. By default, UVBYBETA does not display
; information if output variables are supplied (and TEXTOUT is
; not set).
;
; OPTIONAL OUTPUTS:
; Te - approximate effective temperature
; MV - absolute visible magnitude
; Eby - Color excess E(b-y)
; delm0 - metallicity index, delta m0, (may not be calculable for early
; B stars).
; radius - Stellar radius (R/R(solar))
; EXAMPLE:
; Suppose 5 stars have the following Stromgren parameters
;
; by = [-0.001 ,0.403, 0.244, 0.216, 0.394 ]
; m1 = [0.105, -0.074, -0.053, 0.167, 0.186 ]
; c1 = [0.647, 0.215, 0.051, 0.785, 0.362]
; hbeta = [2.75, 2.552, 2.568, 2.743, 0 ]
; nn = [1,2,3,7,8] ;Processing group number
;
; Determine stellar parameters and write to a file uvbybeta.prt
; IDL> uvbybeta, by,m1,c1,hbeta, nn, t=3
; ==> E(b-y) = 0.050 0.414 0.283 0.023 -0.025
; Teff = 13060 14030 18420 7250 5760
; M_V = -0.27 -6.91 -5.94 2.23 3.94
; radius= 2.71 73.51 39.84 2.02 1.53
; SYSTEM VARIABLES:
; The non-standard system variables !TEXTOUT and !TEXTUNIT will be
; automatically defined if they are not already present.
;
; DEFSYSV,'!TEXTOUT',1
; DEFSYSV,'!TEXTUNIT',0
;
; NOTES:
; (1) **This procedure underwent a major revision in January 2002
; and the new calling sequence may not be compatible with the old** (NAME
; is now a keyword rather than a parameter.)
;
; (2) Napiwotzki et al. (1993, A&A, 268, 653) have written a FORTRAN
; program that updates some of the Moon (1985) calibrations. These
; updates are *not* included in this IDL procedure.
; PROCEDURES USED:
; DEREDD, TEXTOPEN, TEXTCLOSE
; REVISION HISTORY:
; W. Landsman IDL coding February, 1988
; Keyword textout added, J. Isensee, July, 1990
; Made some constants floating point. W. Landsman April, 1994
; Converted to IDL V5.0 W. Landsman September 1997
; Added Eby_in, /PROMPT keywords, make NAME a keyword and not a parameter
; W. Landsman January 2002
; NAME:
; WCSSPH2XY
; PURPOSE:
; Convert spherical coordinates to x and y (map) angular coordinates
; EXPLANATION:
; Convert spherical (longitude and latitude -- sky) coordinates to x
; and y (map) angular coordinates. This procedure is the inverse of
; WCSXY2SPH. See WCS_DEMO for example of use.
;
; This is a lower level procedure -- given a FITS header, the user will
; usually use ADXY which will then call WCSSPH2XY with the appropriate
; parameters.
; CATEGORY:
; Mapping and Auxiliary FITS Routine
;
; CALLING SEQUENCE:
; wcssph2xy, longitude, latitude, x, y, [ map_type , CTYPE = ,
; FACE =,PROJP1 = , PROJP2= , CRVAL = , CRXY = , LONGPOLE = ,
; LATPOLE = , NORTH_OFFSET =, SOUTH_OFFSET =, BADINDEX =]
;
; INPUT PARAMETERS:
; longitude - longitude of data, scalar or vector, in degrees
; latitude - latitude of data, same number of elements as longitude,
; in degrees
; map_type - optional positional parameter, numeric scalar (0-25)
; corresponding to a particular map projection. This is not a
; FITS standard, it is simply put in to allow function similar
; to that of less general map projection procedures (eg AITOFF).
; The following list gives the map projection types and their
; respective numbers.
;
; FITS Number Name Comments
; code code
; ---- ------ ----------------------- -----------------------------------
; DEF 0 Default = Cartesian
; AZP 1 Zenithal perspective projp1 required
; TAN 2 Gnomic AZP w/ projp1 = 0
; SIN 3 Orthographic AZP w/ projp1 = Infinity (>10^14)
; STG 4 Stereographic AZP w/ projp1 = 1
; ARC 5 Zenithal Equidistant
; ZPN 6 Zenithal polynomial prop1-projp9 required, useless
; ZEA 7 Zenithal equal area
; AIR 8 Airy projp1 required
; CYP 9 Cylindrical perspective projp1 and projp2 required
; CAR 10 Cartesian
; MER 11 Mercator
; CEA 12 Cylindrical equal area projp1 required
; COP 13 Conical perspective projp1 and projp2 required
; COD 14 Conical equidistant projp1 and projp2 required
; COE 15 Conical equal area projp1 and projp2 required
; COO 16 Conical orthomorphic projp1 and projp2 required
; BON 17 Bonne's equal area projp1 required
; PCO 18 Polyconic
; SFL 19 Sanson-Flamsteed
; PAR 20 Parabolic
; AIT 21 Hammer-Aitoff
; MOL 22 Mollweide
; CSC 23 Cobe Quadrilateralized convergence of inverse is poor
; Spherical Cube
; QSC 24 Quadrilateralized
; Spherical Cube
; TSC 25 Tangential Spherical Cube
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;
; CTYPE - One, two, or three element vector containing 8 character
; strings corresponding to the CTYPE1, CTYPE2, and CTYPE3
; FITS keywords:
;
; CTYPE[0] - first four characters specify standard system
; ('RA--','GLON' or 'ELON' for right ascension, Galactic
; longitude or ecliptic longitude respectively), second four
; letters specify the type of map projection (eg '-AIT' for
; Aitoff projection)
; CTYPE[1] - first four characters specify standard system
; ('DEC-','GLAT' or 'ELAT' for declination, galactic latitude
; or ecliptic latitude respectively; these must match
; the appropriate system of ctype1), second four letters of
; ctype2 must match second four letters of ctype1.
; CTYPE[2] - if present must be the 8 character string,'CUBEFACE',
; only used for spherical cube projections to identify an axis
; as containing the face on which each x and y pair of
; coordinates lie.
; PROJP1 - scalar with first projection parameter, this may
; or may not be necessary depending on the map projection used
; PROJP2 - scalar with second projection parameter, this may
; or may not be necessary depending on the map projection used
; CRVAL - 2 element vector containing standard system coordinates (the
; longitude and latitude) of the reference point
; CRXY - 2 element vector giving the x and y coordinates of the
; reference point, if this is not set the offset is [0,0]
; This is not a FITS standard -- it is similar to CRPIX but in
; angular X,Y coordinates (degrees) rather than pixel coordinates
; LATPOLE - native latitude of the standard system's North Pole
; LONGPOLE - native longitude of standard system's North Pole, default
; is 180 degrees for Zenithal systems
; NORTH_OFFSET - offset (radians) added to input points near north pole.
; SOUTH_OFFSET - offset (radians) added to input points near south pole.
; BADINDEX - vector, list of transformed points too close to poles.
;
;
; OUTPUT PARAMETERS:
;
; x - x coordinate of data, same number of elements as longitude, in
; degrees; if CRXY is set, then x will be returned offset by
; crxy(0). NOTE: x in all map projections increases to the
; left, not the right.
; y - y coordinate of data, same number of elements as longitude, in
; degrees; if CRXY is set, y will be returned offset by crxy[1]
; bad - vector returning index to transformed points close to pole.
;
; OPTIONAL OUTPUT KEYWORD PARAMETERS:
; FACE - a output variable used for spherical cube projections to
; designate the face of the cube on which the x and y
; coordinates lie. Will contain the same number of elements as
; X and Y. Must contain at least 1 arbitrary element on input
; If FACE is NOT defined on input, it is assumed that the
; spherical cube projection is laid out over the whole sky
; in the "sideways T" configuration.
; NOTES:
; The conventions followed here are described in more detail in
; "Representations of Celestial Coordinates in FITS" by Mark Calabretta
; and Eric Greisen (2002, A&A, 395, 1077; also see
; http://www.aoc.nrao.edu/~egreisen). The general
; scheme outlined in that article is to first use WCS_ROTATE to convert
; coordinates in one of three standard systems (celestial, galactic,
; or ecliptic) into a "native system" of latitude and longitude. The
; latitude and longitude are then converted into x and y coordinates
; which depend on the map projection which is performed. The rotation
; from standard to native coordinates can be skipped if one so desires.
; This procedure necessitates two basic sections. The first converts
; "standard" coordinates to "native" coordinates while the second converts
; "native" coordinates to x and y coordinates. The first section is
; simply a call to WCS_ROTATE, while the second contains the guts of
; the code in which all of the map projection is done. This procedure
; can be called in a form similar to AITOFF, EQPOLE, or QDCB by calling
; wcssph2xy with a fifth parameter specifying the map projection by
; number and by not using any of the keywords related to the map
; projection type (e.g. CTYPE).
;
; PROCEDURE:
;
; The first task of the procedure is to do general error-checking to
; make sure the procedure was called correctly and none of the
; parameters or keywords conflict. This is particularly important
; because the procedure can be called in two ways (either using
; FITS-type keywords or using a number corresponding to a map projection
; type). All variables are converted into double precision values and
; angular measurements are converted from degrees into radians.
; If necessary, longitude values are converted into the range -pi to pi.
; Any latitude points close to the of the poles are mapped to a specific
; latitude of from the pole so that the map transformations become
; completely invertible. The magnitude of this correction is given by
; the keywords NORTH_OFFSET and SOUTH_OFFSET and a list of affected
; points is optionally returned in the "badindex" output parameter.
; The next task of the procedure is to convert the "standard"
; coordinates to "native" coordinates by rotating the coordinate system.
; This rotation is performed by the procedure WCS_ROTATE and is governed
; by the keywords CRVAL and LONGPOLE. The final task of the WCSSPH2XY
; is to take "native" latitude and longitude coordinates and convert
; them into x and y coordinates. Any map specific error-checking is
; done at this time. All of the equations were obtained from
; "Representations of Celestial Coordinates in FITS" and cases needing
; special attention are handled appropriately (see the comments with
; individual map projections for more information on special cases).
;
; Note that a further transformation (using the CD matrix) is required
; to convert the (x,y) coordinates to pixel coordinates.
; COMMON BLOCKS:
;
; none
;
; PROCEDURES CALLED:
; WCS_ROTATE
;
; COPYRIGHT NOTICE:
;
; Copyright 1993, The Regents of the University of California. This
; software was produced under U.S. Government contract (W-7405-ENG-36)
; by Los Alamos National Laboratory, which is operated by the
; University of California for the U.S. Department of Energy.
; The U.S. Government is licensed to use, reproduce, and distribute
; this software. Neither the Government nor the University makes
; any warranty, express or implied, or assumes any liability or
; responsibility for the use of this software.
;
; AUTHOR:
;
; Rick Balsano
;
; MODIFICATIONS/REVISION LEVEL:
;
; 1.1 8/31/93
; 2.3 9/15/93 W. Landsman (HSTX) Update quad cube coords, vectorize
; keywords
; 2.4 12/29/93 I. Freedman (HSTX) Eliminated LU decomposition
; 2.5 1/5/93 I. Freedman (HSTX) Offset keywords / bad point index
; 2.6 Dec 94 Compute pole for transformations where the reference
; pixel is at the native origin W. Landsman (HSTX)
; 2.7 May 95 Change internal variable BETA for V4.0 compatibility
; 2.8 June 95 Change loop indices from integer to long
; 2.9 3/18/96 Change FACE usage for cube projections to match WCSLIB
; C/FORTRAN software library.
; Converted to IDL V5.0 W. Landsman September 1997
; 2.10 02/18/99 Fixed implementation of ARC algorithm
; 2.11 June 2003 Update conic projections, add LATPOLE keyword
; 2.12 Aug 2003, N.Rich - Fix pre-V5.5 bug from previous update
; 2.13 Sep 2003, W. Landsman CTYPE keywords need not be 8 characters
; NAME:
; WCSXY2SPH
;
; PURPOSE:
; Convert x and y (map) coordinates to spherical coordinates
; EXPLANATION:
; To convert x and y (map) coordinates to spherical (longitude and
; latitude or sky) coordinates. This procedure is the inverse of
; WCSSPH2XY.
;
; This is a lower level procedure -- given a FITS header, the user will
; usually use XYAD which will then call WCSXY2SPH with the appropriate
; parameters.
; CATEGORY:
; Mapping and Auxilary FITS Routine
;
; CALLING SEQUENCE:
;
; wcsxy2sph, x, y, longitude, latitude, [map_type], [ CTYPE = ,$
; FACE = ,PROJP1 = , PROJP2 = ,CRVAL =, CRXY =, LONGPOLE=, LATPOLE=]
;
; INPUT PARAMETERS:
;
; x - x coordinate of data, scalar or vector, in degrees, NOTE: x
; increases to to the left, not the right
; y - y coordinate of data, same number of elements as x, in degrees
; map_type - optional positional parameter, scalar corresponding to a
; particular map projection. This is not a FITS standard, it is
; simply put in to allow function similar to that of less general
; map projection procedures (eg AITOFF). The following list gives
; the map projection types and their respective numbers.
;
; FITS Number Name Comments
; code code
; ---- ------ ----------------------- -----------------------------------
; DEF 0 Default = Cartesian
; AZP 1 Zenithal perspective projp1 required
; TAN 2 Gnomic AZP w/ projp1 = 0
; SIN 3 Orthographic AZP w/ projp1 = Infinity (>10^14)
; STG 4 Stereographic AZP w/ projp1 = 1
; ARC 5 Zenithal Equidistant
; ZPN 6 Zenithal polynomial prop1-projp9 required, useless
; ZEA 7 Zenithal equal area
; AIR 8 Airy projp1 required
; CYP 9 Cylindrical perspective projp1 and projp2 required
; CAR 10 Cartesian
; MER 11 Mercator
; CEA 12 Cylindrical equal area projp1 required
; xy 13 Conical perspective projp1 and projp2 required
; COD 14 Conical equidistant projp1 and projp2 required
; COE 15 Conical equal area projp1 and projp2 required
; COO 16 Conical orthomorphic projp1 and projp2 required
; BON 17 Bonne's equal area projp1 required
; PCO 18 Polyconic
; SFL 19 Sanson-Flamsteed
; PAR 20 Parabolic
; AIT 21 Hammer-Aitoff
; MOL 22 Mollweide
; CSC 23 Cobe Quadrilateralized inverse converges poorly
; Spherical Cube
; QCS 24 Quadrilateralized
; Spherical Cube
; TSC 25 Tangential Spherical Cube
;
; OPTIONAL KEYWORD PARAMETERS:
;
; CTYPE - One, two, or three element vector containing 8 character
; strings corresponding to the CTYPE1, CTYPE2, and CTYPE3
; FITS keywords:
;
; CTYPE[0] - first four characters specify standard system
; ('RA--','GLON' or 'ELON' for right ascension, galactic
; longitude or ecliptic longitude respectively), second four
; letters specify the type of map projection (eg '-AIT' for
; Aitoff projection)
; CTYPE[1] - first four characters specify standard system
; ('DEC-','GLAT' or 'ELAT' for declination, galactic latitude
; or ecliptic latitude respectively; these must match
; the appropriate system of ctype1), second four letters of
; ctype2 must match second four letters of ctype1.
; CTYPE[2] - if present must be the 8 character string,'CUBEFACE',
; only used for spherical cube projections to identify an axis
; as containing the face on which each x and y pair of
; coordinates lie.
; FACE - a input variable used for spherical cube projections to
; designate the face of the cube on which the x and y
; coordinates lie. Must contain the same number of elements
; as X and Y.
; CRVAL - 2 element vector containing standard system coordinates (the
; longitude and latitude) of the reference point
; CRXY - 2 element vector giving the x and y coordinates of the
; reference point, if this is not set the offset of the x
; coordinate is assumed to be 0.
; LATPOLE - native latitude of the standard system's North Pole
; LONGPOLE - native longitude of standard system's North Pole, default
; is 180 degrees, numeric scalar
; PROJP1 - scalar with first projection parameter (PV2_1), this may
; or may not be necessary depending on the map projection used
; PROJP2 - scalar with second projection parameter (PV2_2), this may
; or may not be necessary depending on the map projection used
;
; OUTPUT PARAMETERS:
;
; longitude - longitude of data, same number of elements as x, in degrees
; latitude - latitude of data, same number of elements as x, in degrees
;
; NOTES:
; The conventions followed here are described in more detail in the paper
; "Representations of Celestial Coordinates in FITS" by Calabretta &
; Greisen (2002, A&A, 395, 1077, also see
; http://www.aoc.nrao.edu/~egreisen). The general scheme
; outlined in that article is to convert x and y coordinates into a
; "native" longitude and latitude and then rotate the system into one of
; three generally recognized systems (celestial, galactic or ecliptic).
;
; This procedure necessitates two basic sections. The first converts
; x and y coordinates to "native" coordinates while the second converts
; "native" to "standard" coordinates. The first section contains the
; guts of the code in which all of the map projection is done. The
; second step is performed by WCS_ROTATE and only involves rotation of
; coordinate systems. WCSXY2SPH can be called in a form similar to
; AITOFF, EQPOLE, or QDCB by calling wcsxy2sph with a fifth parameter
; specifying the map projection by number and by not using any of the
; keywords related to the map projection type (eg ctype1 and ctyp2).
;
; PROCEDURE:
; The first task of the procedure is to do general error-checking to
; make sure the procedure was called correctly and none of the
; parameters or keywords conflict. This is particularly important
; because the procedure can be called in two ways (either using
; FITS-type keywords or using a number corresponding a map projection
; type). All variables are converted into double precision values.
;
; The second task of the procedure is to take x and y coordinates and
; convert them into "native" latitude and longitude coordinates.
; Map-specific error-checking is done at this time. All of the
; equations were obtained from "Representations of Celestial
; Coordinates in FITS" and cases needing special attention are handled
; appropriately (see the comments with individual map projections for
; more information on special cases). WCS_ROTATE is then called to
; convert the "native" coordinates to "standard" coordinates by rotating
; the coordinate system. This rotation is governed by the keywords
; CRVAL, and LONGPOLE. The transformation is a straightforward
; application of euler angles. Finally, longitude values are converted
; into the range from 0 to 360 degrees.
;
; COMMON BLOCKS:
; none
; PROCEDURES CALLED:
; WCS_ROTATE
;
; COPYRIGHT NOTICE:
;
; Copyright 1991, The Regents of the University of California. This
; software was produced under U.S. Government contract (W-7405-ENG-36)
; by Los Alamos National Laboratory, which is operated by the
; University of California for the U.S. Department of Energy.
; The U.S. Government is licensed to use, reproduce, and distribute
; this software. Neither the Government nor the University makes
; any warranty, express or implied, or assumes any liability or
; responsibility for the use of this software.
;
; AUTHOR:
;
; Rick Balsano
;
; MODIFICATIONS/REVISION LEVEL:
;
; 1.1 8/31/93
; 1.2 9/12/93 W. Landsman Vectorized CRXY, CRVAL, CTYPE
; 1.3 29/12/93 I. Freedman Eliminated LU decomposition
; 1.4 22/09/94 W. Landsman If scalar input, then scalar output
; 1.5 02/03/05 W. Landsman Change variable name BETA for V4.0 compatibility
; 1.6 06/07/05 W. Landsman Change loop index from integer to long
; Converted to IDL V5.0 W. Landsman September 1997
; 1.7 02/18/99 W. Landsman Fixed implementation of ARC algorithm
; 1.8 June 2003 W. Landsman Update conic projections, add LATPOLE keyword
; 1.81 Sep 2003 W. Landsman Avoid divide by zero
; 1.82 Sep 2003 W. Landsman CTYPE keywords need not be 8 characters
; 1.83 Sep 2003 W. Landsman Preserve input array sizes