Viewing contents of file '../idllib/contrib/buie/astrom.pro'
;+
; NAME:
; astrom
; PURPOSE:
; Astrometry from a digital image.
; DESCRIPTION:
;
; This program is designed to permit doing astrometry and catalog driven
; photometry of digital images. It is implicitly assumed that these are
; true digital images, ie., that the images are strictly linear up to some
; signal level. The possiblities supported by this program are quite
; extensive. Read the PROCEDURE section below for more details.
;
; CATEGORY:
; Astrometry
; CALLING SEQUENCE:
; astrom,root,fileno
;
; INPUTS:
; root - Root of data file name (ie., 970309), must be a string.
; fileno - File number to load (0-999), integer.
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD INPUT PARAMETERS:
; BINFAC - Binning factor for displayed image. Default=2
;
; BORDER - Optional inset from each edge of array to be considered valid
; for measurements. The default is 20 pixels in from each
; edge. The value can either be scalar which is applied to
; all edges, or a 4 element vector that specifies the inset
; relative to [left,right,top,bottom]
;
; CATPATH - Location for finding USNO A1.0 catalog for use with refnet.
; default takes on the default of refnet.
;
; CENTER - Controls what is used for center of image for star extraction
; 0 - (default) with objects, uses object ephemeris as center
; no objects, interactive query for choice.
; 1 - Use the header coordinates, no corrections.
; 2 - Use the header coordinates after taking out last known
; position offset.
;
; CLEAN - Flag, if set requests filtering out
; cosmic ray strikes. This step doesn't take too long on
; a Sun Ultra 1/170 but may be prohibitive on slower machines.
;
; DRTHRESH - Threshold on the radial distance (in pixels) of the catalog
; to source matching step when working with a pre-computed
; source list. (Default=4.3)
;
; EDIT - Flag, if set allows interactive culling of bad astrometric
; measurements in the reference net.
;
; EXTLIST - If image is a multi-extension FITS image, this list will
; force the reduction of only the extension numbers listed.
; The default is to do all the extensions, one at a time.
;
; ETATERMS - Vector controlling the astrometric fit in the ETA direction (Dec).
; Must be a 10-element vector, default=[1,1,1,0,0,0,0,0,0,0]
; which is a pure linear fit. See ASTTERMS for a description
; of the terms available.
;
; GAIN - Gain (e-/ADU) of the image. Default=1.0
;
; KEYLIST - Name of a file containing a correspondence list. This list
; associates a set of standard names with the actual keyword
; names found in a FITS file header. If this keyword is
; omitted, a default list is used, as if a file with the
; following contents had been supplied:
; AIRMASS K AIRMASS
; DATE K DATE-OBS
; DATETMPL T DD-MM-YYYY
; EXPDELTA V 0.0
; EXPTIME K EXPTIME
; FILTER K FILTERS
; FILENAME K CCDFNAME
; OBJECT K OBJECT
; UT K UT
; RA K RA
; DEC K DEC
; The middle column is a flag. It may be K, for Keyword,
; T, for Template, or V, for Value. If it is V, the contents
; of the third field on that line should make sense for the
; name in the first field.
;
; MAGLIM - Limiting (faint) magnitude for catalog extraction (default=30.0)
;
; MINREFRESH - Flag, if set will minimize the amount of screen refreshes.
;
; MAXSTARS - Maximum number of reference stars to measure. (default=all)
;
; NEWCAT - Flag, if set forces the recreation of the STARFILE.
;
; NOSCALE - Flag, if set suppresses conversion from integer to floating point
; when the FITS file is read.
;
; NOOBJECTS - Flag, if set suppresses the final step of measuring astrometric
; unknowns.
;
; NOREMIND - Flag, if set suppresses the query for the reminder location.
;
; OBJNAME - Name of object to collect astrometry for. By default, the
; FITS header from the image is scanned for the OBJECT keyword
; and this becomes the OBJNAME after compressing multiple blanks,
; trimming leading and trailing blanks, and replacing single
; blanks by "_".
;
; OBJRAD - Radius (in pixels) of object aperture for astrometry and
; photometry. Default=10. The sky aperture is set between
; objrad+10 and objrad+30.
;
; PATH - String, this is the name of the directory where the data are
; stored. The actual data directory used is PATH+'/'+root.
; The default is '' (blank) and the file would be root.NNN
; which would permit putting a leading path on the root.
;
; PLASTFILE - This gives the name of the "plast" output file. The default
; is OBJNAME.pla. This file contains a list of asteroids that
; may be found on the image and is created by a separate
; program. To disable this feature, set PLASTFILE='none'.
;
; PHOTSTARS - Flag. If set, turns on a special mode that performs photometry
; on all good astrometric reference stars. The photometry is
; added to the root+'.log' file. Multiple reductions are
; weeded out and the photometry log file is left sorted by
; file. The log file is in the ALTLOG format (see basphote).
; The star name is automatically created from its coordinate.
; Ex: a star at 12:12:34.2 and +04:23:45 would be named
; NV1212342+042345. The position used for the name comes from
; the catalog, not astrometry from the image.
;
; PRETTY - String, if set, will cause the final image with overlays to
; be sent to a color postscript output file. This will only
; be used if you are measuring objects (id., NOOBJECTS not set).
; If not specified, this file will not be created.
;
; QUEUE - String, name of printer to send output to. If supplied, a
; hardcopy of the image is generated along with a list of
; stars identified. If blank or not specified, no printed output
; is generated.
;
; RESFILE - Filename where astrometric measurements are written to. The
; default file name is OBJNAME.ast. Only one line per image
; is allowed. Subsequent measurements made by astrom will
; override the measures for this image.
;
; ROAM - Flag, if set will provide a continuous running display of the
; RA and Dec of the cursor when you are in the object
; measurement loop.
;
; SAVECLEAN - Flag, if set will save the cleaned image to disk.
;
; SPOT - Array containing explicit x,y coordinates. Default=none.
; After the astrometric fit is complete, the RA,DEC of each
; x,y pair (2xN array) is computed, printed to the screen, and
; saved to an ancillary file, spot.dat, in the current directory.
;
; STARFILE - Filename where a list of astrometric reference stars is to
; be found. The default file name is OBJNAME.cat. If this file
; is not found, then this program will attempt to create the
; file by calling "refnet", a program that accesses the
; USNO A1.0 star catalog provided by David Monet of USNOFS.
;
; SUBEXP - This keyword controls reducing images with multiple exposures.
; This keyword should contain one or more strings that will
; serve to identify the multiple exposures. Ideally, this
; id string would be a single character, eg., 'a', 'b', etc.
; This program will loop over the string list for multiple
; reductions of the frame. The id string will be appended to
; the frame # where ever it is used. So, in the .ast file and
; fitcoeff.dat file the file name will be root.suffix_tag. In
; Refstars, the files are suffix_tag.ref. The default is to
; process one exposure per image and the _tag will not be added
; to any names.
;
; TRUSTCENTER - Flag, if set indicates that the 2-star solution and the
; previously known plate center are to be trusted. This
; removes the need to do the catalog star identifications.
; If the astrom.inf file is not found, or, if the plate
; center is not found in the centers.dat file, then this
; flag is ignored. This flag is also ignored if there
; is a valid Refstars file.
;
; TWOSTAR - Flag, if set suppresses the automatic correlation of the source
; list and catalog and switches to an interactive two star
; solution to get the initial image location. This should
; be done on the first frame of a night. It also seems to
; be necessary for large fields with non-linear distortions.
; This flag has no effect if there are no pre-existing
; source lists generated by findsrc.pro.
;
; XCENTER - Optional override of location of center of optical axis in
; pixel coordinates. The default is the center of the array.
; This location is considered the location of the tangent plane
; and the location of x=y=0 for the (x,y) <--> (xi,eta)
; transformation.
;
; XITERMS - Vector that controls the astrometric fit in the XI direction (RA)
; Must be a 10-element vector, default=[1,1,1,0,0,0,0,0,0,0]
; which is a pure linear fit. See ASTTERMS for a description
; of the terms available.
;
; WINDOW - window number to display image into.
;
; YCENTER - Optional override of location of center of optical axis in
; pixel coordinates. The default is the center of the array.
; This location is considered the location of the tangent plane
; and the location of x=y=0 for the (x,y) <--> (xi,eta)
; transformation.
;
; OUTPUTS:
; output is graphical and to a series of files.
; astrom.inf - Records the last 2-star astrometric solution.
; If the image being reduced is a multi-extension FITS
; file, this file will be named astromNN.inf where
; NN is the image extension of interest.
; centers.dat - Records the image center for all measured frames.
; root.log - Photometry (if PHOTSTARS set)
; objname.ast - Astrometry of object
; objname.cat - Star catalog extraction
; objname.pla - Output of PLAST (list of asteroids on image).
; root.stars - List of stars and coordinates for those where photometry
; was measured. Intended for inclusion in starcat.dat
; (see GETSTARS.PRO or LOADSTAR.PRO).
; Refstars/fileno.ref - Binary file containing positions, mag, fwhm for
; all catalog stars measured. This file will be
; be reused in later runs of ASTROM on this image
; as long as the object aperture radius is the same.
; fitcoeff.dat- List of fit coefficients for each of the xi,eta axes.
; position.dat- List of x,y positions for all objects measured.
;
; KEYWORD OUTPUT PARAMETERS:
; None.
;
; COMMON BLOCKS:
; None.
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; Input files must all be FITS and the file names must be of the form:
; root.NNN where "root" is some string and NNN is a 3-digit number.
;
; PROCEDURE:
;
; This program automates astrometric reductions of CCD images. Once the
; astrometric solution is determined for the image, you can then proceed
; to measure any source in the image to ascertain it's position. As you
; might guess from the above list, there are far too many options to this
; program. It is rare that you will use all the options, instead, some
; subset can be tweaked and tuned to _your_ data to make the process run
; as quickly as possible and with as little user interaction as possible.
; Under certain circumstances this program can run complete automatically
; but only if a great deal is already known about the images.
;
; To illustrate one use of the program consider doing astrometry of 1 or
; more objects on an image. Typically on the first invocation on a new
; image that you know little about you will not use any optional information.
; However, the first and most important optional input you can provide is
; the FITS keyword correspondence information through the KEYLIST option.
; If you have a decent header, then the program will have a good object name,
; time for the image, and (hopefully) a good coordinate for the image center.
;
; The object name is important because it is used to form file names for
; the output astrometry and the ancillary star catalog. If the images you
; are reducing have the same name for different sky locations then this
; program will get hopelessly confused.
;
; Next, ASTROM will look for a special file, astrom.inf, that contains clues
; about the image scale and orientation. If not found, you will be prompted
; for the needed information. This file can be edited after creation if you
; need to try other guesses (such as flipping image, trying different scales,
; trying different rotation angles, etc.). This step can be very frustrating
; if you don't have much information about your image. I quite often find
; it necessary to edit this file many times and re-run ASTROM on the same
; image until it makes sense.
;
; Next, ASTROM will look for a list of stars that should be on the frame and
; will serve as the astrometric reference network. If found, the file is
; read. If not found, then ASTROM will try to create this file using another
; external program (refnet). The set of stars requested will depend on the
; scale and orientation known at this point so if you get it wrong, delete the
; the star catalog and start over. The center for the star search comes
; from a number of places. This is what ASTROM tries, in order, (1)
; using the object name (first one if given an array), try to get an
; ephemeris position for the object for the time of the image, if this
; fails --> (2) ask for RA and DEC of image center. (3) If NOOBJECTS is set,
; (1) and (2) are bypassed and the RA,DEC from the header are used (after
; precessing to J2000, if needed). If you aren't doing objects, and the
; header value is bad, then you may need to insert your guess for the
; center directly into the "centers.dat" file. Sometimes this is the only
; way to proceed if the headers are really screwy.
;
; Next, ASTROM will try to get a list of
; known asteroids that _might_ be on the image. For this to work, you
; must have an external program that ASTROM will invoke to collect this
; information to a file. If the file already exists, then it will be
; read directly without the need for the external program. Note: this
; external program is non-trivial and cannot be easily exported away from
; Lowell Observatory. Fortunately, it can be disabled, but, if you can
; use it then you will see positional overlays of the expected locations
; of any asteroids on the frame along with their line of variation scaled
; by the orbit uncertainty (green line) and a 1-hour motion vector (yellow
; line).
;
; With this information in hand, ASTROM begins to work with the screen image.
; The image is displayed with a fairly hard linear stretch, -7 sigma to
; +16 sigma about the mean sky signal. The border is drawn and you are
; (possibly) asked to indicate a location on the image to be remembered. I
; find this useful in marking the location of a specific object that will
; be highlighted throughout the analysis. This is not used with NOOBJECTS
; set.
;
; If you allow it, ASTROM will proceed to remove cosmic rays from the image.
; This step is almost always useful but can sometimes take a long time to
; run depending on the size of the image and the speed of your computer.
; On a SUN Ultra 1/170E, a 4k x 2k image can take 10-20 seconds to clean.
; After cleaning the image is redisplayed with the same scaling. If you watch
; carefully you will see the cosmic ray strikes "blink off".
;
; Next, the asteroid overlay is generated and plotted. This involves
; generating ephemeris positions for all the asteroids for the time of this
; frame. Again, an external program is called to generate these positions
; and you find the relevant documentation in my IDL front-end program, EPHEM.
;
; Now, we get on to the steps of getting the astrometric grid in place.
; The image center, scale and orientation are used in calculating the locations
; of all catalog stars. Red diamonds, scaled in size by magnitude, are
; plotted on the image at these predicted locations. Your job at this
; point is to match up the overlay with the image. This can be either
; very easy, very hard, or anywhere in between. You may find need to
; tweak the contents of astrom.inf, change MAGLIM, or more depending on
; the situation. Sometimes the center is no good. One of two things will
; happen, (1) the wrong star overlay is plotted on the image, or (2) no
; stars are plotted. In the latter case, you will be prompted for a new
; image center. If you're lucky, this will allow you to proceed.
;
; Assuming all goes well, you now must establish the correspondence between
; the overlay and the image. At this step you must identify two stars
; in the image and the same two stars in the overlay. The prompts from ASTROM
; will indicate what information it next desires. Just remember that at
; any time ASTROM is looking for a mouse click to proceed, a right click will
; exit the program directly at that point. After you identify the first star,
; the overlay is replotted so that the overlay star sits over the image.
; After you identify the second star ASTROM has enough information to predict
; the location of all the other catalog stars and can proceed to measure
; their locations. Once all the catalog stars are measured, an astrometric
; function is fit to the positions. See, ASTTERMS for all the possible
; terms. The default is a linear solution which is usually pretty good for
; most CCD images. At this step you may need to fiddle with the aperture
; radius used to measure the stars. A radius (OBJRAD) that is too big
; or too small can lead to excessive scatter in the astrometry. I usually
; try to set OBJRAD=fwhm or just slightly under the fwhm. These positions
; the subsequent fit are saved to a couple of files so that if you come
; back to the image later, the fit is regenerated much quicker and you get
; directly to the next step, that of measuring unknown objects.
;
; Note that the fit to the star positions is done in a robust fashion. Stars
; with large residuals, unusual fwhm, signal too weak, or signal too strong,
; are avoided in the fit. These stars will plot as red circles in the final
; overlay while the good stars will overplot with yellow circles. You can
; control the "too strong" threshold with MAXPHOTSIG or you can do a purely
; interactive editing of the star positions used with /EDIT. Using /EDIT
; is required only if there just aren't that many stars (say 3-6) and there's
; no statistical basis for chucking out anything. If there are lots of stars,
; the automatic stuff works just fine.
;
; If you find that the header gives consistently good predictions of the
; catalog star locations, then use the TRUSTCENTER keyword to bypass the
; 2-star interactive step. Be careful, you need to have pretty good
; for this to work effectively. Also note that if you turn on the PHOTSTARS
; flag, ASTROM will automatically collect and save aperture photometry data
; on _all_ the catalog stars.
;
; At this point, you've arrived at the time you can measure new objects.
; You are prompted to click left on the object to measure. Nothing is saved
; until you click middle (to go to next object) or click right (to quit).
; So you can feel free to poke around in the image measuring lots of objects
; without saving all of it. When you click left it measures the location
; and computes RA,DEC, and, it puts up a small window on the object with the
; aperture location overplotted. In this window, there is a non-linear
; stretch so that you can see the wings of the PSF as well as the core.
; You also get a radial profile for additional diagnostic information.
; If you continue with a middle click, ASTROM will step through the OBJNAME
; array (if provided) or when it runs of out names it knows, it will begin
; prompting you for a name. This name is used to form the file name for
; all saved astrometry. Once you have processed an entire nights worth of
; data then you can use the ASTCOL program to collection it all and save it
; to another master storage location. The ASTCOL step is where the
; observatory code information is added to the astrometry.
;
; If you set /PRETTY, the very last thing done is to create a fancy postscript
; output image file. This image has all the overlays and asteroid locations
; and is occasionally useful.
;
; This description is by no means exhaustive. There are a very large number
; of options that must be tweaked to get the most out of the reduction. I
; ususally find it useful to create a master script that contains all the
; flags and options as I develop a means to reduce some data.
;
; Here's one example:
; astrom,fdir,fnum,maglim=17, $
; binfac=1,objrad=objrad,path=d,key='../site.key',objname=o, $
; /noremind,noobjects=noob,gain=2.5,plastfile='none',trustcenter=trust, $
; xi=[1,1,1,0,0,0,0,0,0,0],eta=[1,1,1,0,0,0,0,0,0,0],maxphotsig=23000.0
;
; Here I've set the limiting magnitude for the reference stars to 17.0,
; no cosmic ray cleaning, don't ask for reminder location, set gain
; of CCD, turn off looking for field asteroids, force linear fit, and set
; saturation level of CCD. The other options are either obvious or set
; to fixed values in the script handling this call. Note that by splitting
; the file name into root and suffix, you can generate a loop on the suffix
; and step through all the images you have.
;
; MODIFICATION HISTORY:
; 97/04/05, Written by Marc W. Buie, Lowell Observatory
; 97/06/13, MWB, added keylist fits header reading generalization
; 97/06/14, MWB, fixed line of variations plotting bug, also added saving
; information on reference stars used.
; 97/06/17, MWB, added controls over terms in Xi and Eta fits.
; 97/06/18, MWB, added saving fit coeffcients and object positions.
; 97/06/19, MWB, added NOREMIND keyword
; 97/06/20, MWB, added TRUSTCENTER keyword
; 97/09/09, MWB, added SUBEXP keyword
; 97/10/08, MWB, added X,YCENTER keywords.
; 97/10/16, MWB, added SPOT keyword.
; 97/10/21, MWB, rewrote initial 2-star fit.
; 97/11/24, MWB, added MAXSTARS keyword.
; 98/01/06, MWB, added support for pre-cleaned images.
; 98/03/13, MWB, some heavy rewriting. NOCLEAN is now CLEAN (default none)
; plus internal cleanups, some changes to ancillary plots,
; added support for external lists of image sources.
; 98/06/24, MWB, added TWOSTAR flag.
; 98/08/26, MWB, added ROAM keyword.
; 98/10/07, MWB, a few optimizations and some bug fixes and enhancements
; for multi-extension files.
; 98/10/08, MWB, added DRTHRESH keyword
; 98/11/04, MWB, added NOSCALE keyword
; 98/12/02, MWB, changed usage of PATH so that you can append 'root' or not.
; This now allows the image file name roots to be different
; from the name of the directory they are stored in.
;-
PRO astrom_dim,window,bim,binfac,border
setwin,window
sz=size(bim)
tv,bim
plot,[0],[1],/nodata,xmargin=[0,0],ymargin=[0,0], $
xr=[0,sz[1]*binfac-1],yr=[0,sz[2]*binfac-1],xstyle=5,ystyle=5,/noerase
oplot,[border[0],sz[1]*binfac-border[1],sz[1]*binfac-border[1],border[0],border[0]], $
[border[3],border[3],sz[2]*binfac-border[2],sz[2]*binfac-border[2],border[3]],color=4
END
PRO astrom_svctr,thisexp,info
rastr,info.racen,4,ras
decstr,info.deccen,3,decs
print,'Image@',ras,decs,'scale',info.pscale,'arcsec/pixel', $
info.rang*!radeg,'deg rotation', $
format='(a,1x,a,1x,a,1x,a,f6.3,1x,a,1x,f7.2,1x,a)'
cinfo=thisexp+' '+ras+' '+decs
repwrite,'centers.dat',thisexp,cinfo
END
pro astrom_regen,xpos,ypos,sra,sdec,z,zs1,zs2,info
xs1m = xpos[zs1]
xs2m = xpos[zs2]
ys1m = ypos[zs1]
ys2m = ypos[zs2]
info.raref = sra[zs1]
info.decref = sdec[zs1]
info.xcref = xs1m
info.ycref = ys1m
dx = info.xflip*(xs2m-xs1m)
dy = info.yflip*(ys2m-ys1m)
pixsep = sqrt(dx^2+dy^2)
astrd2sn,sra[zs2],sdec[zs2],info.raref,info.decref,xi2,eta2
xi2 = xi2*180.0d0/!dpi*3600.0d0
eta2 = eta2*180.0d0/!dpi*3600.0d0
angsep = sqrt(xi2^2+eta2^2)
info.pscale = angsep/pixsep
info.rang = atan(dy*xi2-dx*eta2,dx*xi2+dy*eta2)
xiterms0 = [0,1,1,0,0,0,0,0,0,0]
cxi0 = [cos(info.rang),sin(info.rang)]*info.pscale
etaterms0 = [0,1,1,0,0,0,0,0,0,0]
ceta0 = [-sin(info.rang),cos(info.rang)]*info.pscale
dx = info.xflip*(info.xc0 - xs1m)
dy = info.yflip*(info.yc0 - ys1m)
xi0 = asteval(dx,dy,cxi0,xiterms0)/3600.0d0*!dpi/180.0d0
eta0 = asteval(dx,dy,ceta0,etaterms0)/3600.0d0*!dpi/180.0d0
astsn2rd,xi0,eta0,info.raref,info.decref,ra,dec
info.ra0 = ra[0]
info.dec0 = dec[0]
end
PRO astrom_svref,reffile,o,c,x,y,f,m,e,mx,ra,dec,sm
openw,lref,reffile,/get_lun
writeu,lref,float(o),long(c) ; object aperture radius, number of objects.
writeu,lref,float(x) ; x location of source
writeu,lref,float(y) ; y location of source
writeu,lref,float(f) ; fwhm of source
writeu,lref,float(m) ; raw magnitude of source
writeu,lref,float(e) ; uncertainty on source magnitude
writeu,lref,float(mx) ; peak signal on source
writeu,lref,double(ra) ; J2000 ra of source in radians
writeu,lref,double(dec) ; J2000 dec of source in radians
writeu,lref,float(sm) ; Catalog magnitude of source
free_lun,lref
END
; x,y Pixel coordinates of catalog stars
; valid - Flag, true if successful, if not true, do not continue.
PRO astrom_twostar,dv,info,sra,sdec,x,y,ssz,valid
valid=0
retry:
; Force xi,eta reference point to be center of image
info.raref = info.racen
info.decref = info.deccen
info.xcref = info.xc
info.ycref = info.yc
astrd2xy,sra,sdec,info,x,y,XI=sxi,ETA=seta
astrom_dim,dv.win,dv.bim,dv.binfac,dv.border
for i=0,n_elements(x)-1 do $
oplot,[x[i]],[y[i]],psym=4,color=2,symsize=ssz[i]
; Select first anchor star from catalog list
print,' Click left on symbol for first anchor star (right=quit).'
cursor,xs1,ys1,3
if !mouse.button eq 4 then return
dist = sqrt((xs1-x)^2+(ys1-y)^2)
zs1 = where(dist eq min(dist))
zs1 = zs1[0]
oplot,[x[zs1]],[y[zs1]],psym=4,color=4,symsize=2
oldx = x[zs1]
oldy = y[zs1]
; Click on the star in the image
print,' Click left on the star in the image (right=quit).'
cursor,xs1,ys1,3
if !mouse.button eq 4 then return
basphote,dv.gain,dv.image,dv.exptime,xs1,ys1, $
dv.objrad,dv.objrad+10,dv.objrad+30, $
/nolog,/silent,xcen=xs1m,ycen=ys1m
; Update overlay based on first star.
info.raref = sra[zs1]
info.decref = sdec[zs1]
info.xcref = xs1m
info.ycref = ys1m
astrd2xy,sra,sdec,info,x,y,XI=sxi,ETA=seta
if sqrt((oldx-xs1m)^2+(oldy-ys1m)^2) gt 4.0 then begin
astrom_dim,dv.win,dv.bim,dv.binfac,dv.border
for i=0,n_elements(x)-1 do $
oplot,[x[i]],[y[i]],psym=4,color=2,symsize=ssz[i]
oplot,[xs1m],[ys1m],psym=4,color=1
endif
; Select second anchor star from catalog list
print,' Click left on symbol for second anchor star '+ $
'(middle=restart, right=quit).'
cursor,xs2,ys2,3
if !mouse.button eq 2 then goto,retry
if !mouse.button eq 4 then return
dist = sqrt((xs2-x)^2+(ys2-y)^2)
zs2 = where(dist eq min(dist))
zs2 = zs2[0]
oplot,[x[zs2]],[y[zs2]],psym=4,color=4,symsize=2
; Click on the star in the image
print,' Click left on the star in the image (right=quit).'
cursor,xs2,ys2,3
if !mouse.button eq 4 then return
basphote,dv.gain,dv.image,dv.exptime,xs2,ys2,dv.objrad,dv.objrad+10,dv.objrad+30, $
/nolog,/silent,xcen=xs2m,ycen=ys2m
oplot,[xs2m],[ys2m],psym=4,color=1
; Separation in pixels
dx = info.xflip*(xs2m-xs1m)
dy = info.yflip*(ys2m-ys1m)
pixsep = sqrt(dx^2+dy^2)
; Compute standard coordinates of two stars
; for simplicity, force ra0,dec0 to be on star 1, that
; means x1=y1=0 and xi1=eta1=0
astrd2sn,sra[zs2],sdec[zs2],info.raref,info.decref,xi2,eta2
xi2 = xi2*180.0d0/!dpi*3600.0d0
eta2 = eta2*180.0d0/!dpi*3600.0d0
angsep = sqrt(xi2^2+eta2^2)
; plate scale in arcsec/pixel and image rotation angle.
pscale = angsep/pixsep
rang = atan(dy*xi2-dx*eta2,dx*xi2+dy*eta2)
; cast into standard transformation
xiterms0 = [0,1,1,0,0,0,0,0,0,0]
cxi0 = [cos(rang),sin(rang)]*pscale
etaterms0 = [0,1,1,0,0,0,0,0,0,0]
ceta0 = [-sin(rang),cos(rang)]*pscale
; Save values
info.pscale = pscale
info.rang = rang
; Redetermine the ra and dec center
dx = info.xflip*(info.xc - xs1m)
dy = info.yflip*(info.yc - ys1m)
xi0 = asteval(dx,dy,cxi0,xiterms0)/3600.0d0*!dpi/180.0d0
eta0 = asteval(dx,dy,ceta0,etaterms0)/3600.0d0*!dpi/180.0d0
astsn2rd,xi0,eta0,info.raref,info.decref,ra,dec
info.racen = ra[0]
info.deccen = dec[0]
; Redetermine the ra and dec of the optical center
dx = info.xflip*(info.xc0 - xs1m)
dy = info.yflip*(info.yc0 - ys1m)
xi0 = asteval(dx,dy,cxi0,xiterms0)/3600.0d0*!dpi/180.0d0
eta0 = asteval(dx,dy,ceta0,etaterms0)/3600.0d0*!dpi/180.0d0
astsn2rd,xi0,eta0,info.raref,info.decref,ra,dec
info.ra0 = ra[0]
info.dec0 = dec[0]
astrd2xy,sra,sdec,info,x,y,XI=sxi,ETA=seta
valid=1
END
pro astrom,root,fileno, $
WINDOW=window,EDIT=edit,OBJRAD=objrad,BINFAC=binfac, $
STARFILE=starfile,RESFILE=resfile,PATH=in_path,OBJNAME=in_objname, $
CLEAN=clean,PLASTFILE=plastfile,NEWCAT=newcat,PHOTSTARS=photstars, $
QUEUE=queue,NOOBJECTS=noobjects,MAXPHOTSIG=maxphotsig,MAGLIM=maglim, $
KEYLIST=keylist, PRETTY=pretty,BORDER=border,ETATERMS=etaterms, $
XITERMS=xiterms,NOREMIND=noremind,TRUSTCENTER=in_trustcenter,GAIN=gain, $
SUBEXP=subexp,XCENTER=in_xcenter,YCENTER=in_ycenter,SPOT=spot, $
MINREFRESH=minrefresh,MAXSTARS=maxstars,CATPATH=catpath, $
SAVECLEAN=saveclean,EXTLIST=extlist,CENTER=center,TWOSTAR=twostar, $
ROAM=roam,DRTHRESH=drthresh,NOSCALE=noscale
if badpar(root,7,0,CALLER='ASTROM: (root) ') then return
if badpar(fileno,[2,3],0,CALLER='ASTROM: (fileno) ') then return
if badpar(maglim,[0,2,3,4,5],0,CALLER='ASTROM: (MAGLIM) ', $
default=30.0) then return
if badpar(drthresh,[0,2,3,4,5],0,CALLER='ASTROM: (DRTHRESH) ', $
default=4.3) then return
if badpar(binfac,[0,1,2,3],0,CALLER='ASTROM: (BINFAC) ', $
default=2) then return
if badpar(edit,[0,1,2,3],0,CALLER='ASTROM: (EDIT) ',default=0) then return
if badpar(newcat,[0,1,2,3],0,CALLER='ASTROM: (NEWCAT) ', $
default=0) then return
if badpar(clean,[0,1,2,3],0,CALLER='ASTROM: (CLEAN) ', $
default=0) then return
if badpar(noobjects,[0,1,2,3],0,CALLER='ASTROM: (NOOBJECTS) ', $
default=0) then return
if badpar(in_objname,[0,7],[0,1],CALLER='ASTROM: (OBJNAME) ', $
default='[[DEFAULT]]') then return
if badpar(objrad,[0,1,2,3,4,5],0,CALLER='ASTROM: (OBJRAD) ', $
default=10) then return
if badpar(in_path,[0,7],0,CALLER='ASTROM: (PATH) ', $
default='[[DEFAULT]]') then return
if badpar(maxphotsig,[0,2,3,4,5],0,CALLER='ASTROM: (MAXPHOTSIG) ', $
default=13000.0) then return
if badpar(plastfile,[0,7],0,CALLER='ASTROM: (PLASTFILE) ', $
default='[[DEFAULT]]') then return
if badpar(photstars,[0,1,2,3],0,CALLER='ASTROM: (PHOTSTARS) ', $
default=0) then return
if badpar(pretty,[0,7],0,CALLER='ASTROM: (PRETTY) ', $
default='[[DEFAULT]]') then return
if badpar(roam,[0,1,2,3],0,CALLER='ASTROM: (ROAM) ', $
default=0) then return
if badpar(queue,[0,7],0,CALLER='ASTROM: (QUEUE) ', $
default='') then return
if badpar(resfile,[0,7],0,CALLER='ASTROM: (RESFILE) ', $
default='[[DEFAULT]]') then return
if badpar(starfile,[0,7],0,CALLER='ASTROM: (STARFILE) ', $
default='[[DEFAULT]]') then return
if badpar(keylist,[0,7],0,CALLER='ASTROM: (KEYLIST) ', $
default='[[DEFAULT]]') then return
if badpar(window,[0,1,2,3],0,CALLER='ASTROM: (WINDOW) ', $
default=0) then return
if badpar(border,[0,1,2,3],[0,1],CALLER='ASTROM: (BORDER) ', $
default=20,rank=brank) then return
if badpar(xiterms,[0,1,2,3],1,CALLER='ASTROM: (XITERMS) ', $
default=[1,1,1,0,0,0,0,0,0,0],npts=nxiterms) then return
if badpar(etaterms,[0,1,2,3],1,CALLER='ASTROM: (ETATERMS) ', $
default=[1,1,1,0,0,0,0,0,0,0],npts=netaterms) then return
if badpar(saveclean,[0,1,2,3],0,CALLER='ASTROM: (SAVECLEAN) ', $
default=0) then return
if badpar(noremind,[0,1,2,3],0,CALLER='ASTROM: (NOREMIND) ', $
default=0) then return
if badpar(in_trustcenter,[0,1,2,3],0,CALLER='ASTROM: (TRUSTCENTER) ', $
default=0) then return
if badpar(gain,[0,2,3,4,5],0,CALLER='ASTROM: (GAIN) ', $
default=1.0) then return
if badpar(debug,[0,1,2,3],0,CALLER='ASTROM: (DEBUG) ',default=0) then return
if badpar(subexp,[0,7],[0,1],CALLER='ASTROM: (SUBEXP) ', $
default='[[DEFAULT]]') then return
if badpar(spot,[0,2,3,4,5],[1,2],CALLER='ASTROM: (SPOT) ', $
default=[-1.0e9,-1.0e9]) then return
if badpar(minrefresh,[0,1,2,3],0,CALLER='ASTROM: (MINREFRESH) ', $
default=0) then return
if badpar(maxstars,[0,1,2,3],0,CALLER='ASTROM: (MAXSTARS) ', $
default=0) then return
if badpar(extlist,[0,1,2,3],[0,1],CALLER='ASTROM: (EXTLIST) ', $
default=-1) then return
if badpar(center,[0,1,2,3],0,CALLER='ASTROM: (CENTER) ', $
default=0) then return
if badpar(twostar,[0,1,2,3],0,CALLER='ASTROM: (TWOSTAR) ', $
default=0) then return
if badpar(noscale,[0,1,2,3],0,CALLER='ASTROM: (NOSCALE) ', $
default=0) then return
IF nxiterms ne 10 THEN BEGIN
print,'ASTROM: Error! the length of XITERMS must be 10.'
return
ENDIF
IF netaterms ne 10 THEN BEGIN
print,'ASTROM: Error! the length of ETATERMS must be 10.'
return
ENDIF
;===========================================
; This section is only run once upon starting.
maglim=float(maglim) ; force type to float
do_objects = not noobjects
fullrefresh = not minrefresh
ans=''
blanks=' '
bel=STRING( 7B )
; Get header correspondence list.
loadkeys,keylist,hdrlist
; First, detangle the input information and get pointed at the data file
; and directory.
IF in_path eq '[[DEFAULT]]' THEN BEGIN
path = ''
ENDIF ELSE BEGIN
if exists(addslash(in_path)+root+'/') then $
path = addslash(in_path)+root+'/' $
else $
path = addslash(in_path)
ENDELSE
IF not exists('Refstars') THEN spawn,'mkdir Refstars'
IF brank eq 0 THEN BEGIN
border=replicate(border,4)
ENDIF ELSE BEGIN
IF n_elements(border ne 4) THEN BEGIN
print,'ASTROM: Illegal border, must be scalar or 4-element vector.'
print,' [left,right,top,bottom] in pixels from edge.'
return
ENDIF
IF min(border) lt 0 THEN BEGIN
print,'ASTROM: Illegal border, all values must be greater than 0.'
return
ENDIF
ENDELSE
suffix=string(fileno,format='(i3.3)')
; if precleaned version exists we'll use that (regardless of flag).
imfile=root+'c.'+suffix
if not exists(path+imfile) then begin
imfile=root+'.'+suffix
preclean=0
IF not exists(path+imfile) THEN BEGIN
print,'ASTROM: Image ',path+imfile,' not found.'
print,' Aborting.'
return
ENDIF
endif else begin
preclean=1
endelse
; Check header of image to see if it is a multi-extension image.
hdr=headfits(path+imfile)
numext=sxpar(hdr,'NEXTEND')
IF numext eq 0 THEN BEGIN
extlist=0
ENDIF ELSE BEGIN
IF extlist[0] eq -1 THEN BEGIN
extlist=indgen(numext)+1
ENDIF ELSE BEGIN
IF max(extlist) gt numext THEN BEGIN
print,'ASTROM: Input extension list is incompatible with the number of extensions'
print,'in the file. This file had ',numext,' extensions and the largest item in'
print,'your list is ',max(extlist)
print,'Aborting.'
return
ENDIF ELSE IF min(extlist) le 0 THEN BEGIN
print,'ASTROM: Input extension list is invalid. You have one or more values less'
print,'than or equal to zero.'
print,'Aborting.'
return
ENDIF
ENDELSE
ENDELSE
numext=n_elements(extlist)
FOR ext=0,numext-1 DO BEGIN
IF extlist[ext] eq 0 THEN BEGIN
astinf = 'astrom.inf'
exttag = ''
extstr = ''
ENDIF ELSE BEGIN
extstr = strb36(extlist[ext])
astinf = 'astrom'+extstr+'.inf'
exttag = 'x'+extstr
if plastfile ne 'none' then plastfile = '[[DEFAULT]]'
starfile = '[[DEFAULT]]'
ENDELSE
; Read in image, and some header values
print,'Loading ',path+imfile,' ',exttag
raw=0.
raw=readfits(path+imfile,hdr,exten_no=extlist[ext],/silent,noscale=noscale)
parsekey,hdr,hdrlist,hdrinfo
jdstr,hdrinfo.jd,0,jds
sz=size(raw)
IF (sz[1]/binfac)*binfac ne sz[1] or (sz[2]/binfac)*binfac ne sz[2] THEN $
raw=raw[0:(sz[1]/binfac)*binfac-1,0:(sz[2]/binfac)*binfac-1]
IF in_objname[0] eq '[[DEFAULT]]' THEN $
objname=nobname(strtrim(strcompress(hdrinfo.object),2)) $
ELSE $
objname=in_objname
objname=strlowcase(objname)
fnobj = objname[0]+exttag+'.obj'
exptime=hdrinfo.exptime
if exptime eq 0. then exptime=1.0
if badpar(in_xcenter,[0,2,3,4,5],[0,1],CALLER='ASTROM: (XCENTER) ', $
default=sz[1]/2.0) then return
if badpar(in_ycenter,[0,2,3,4,5],[0,1],CALLER='ASTROM: (YCENTER) ', $
default=sz[2]/2.0) then return
if extlist[ext] eq 0 then begin
xcenter = in_xcenter[0]
ycenter = in_ycenter[0]
endif else begin
if n_elements(in_xcenter) eq 1 then $
xcenter=in_xcenter $
else if extlist[ext] le n_elements(in_xcenter) then $
xcenter = in_xcenter[extlist[ext]-1] $
else begin
print,'ASTROM: xcenter does not have enough entries, no value present for'
print,' extension',strcompress(extlist[ext]),' unable to continue.'
endelse
if n_elements(in_ycenter) eq 1 then $
ycenter=in_ycenter $
else if extlist[ext] le n_elements(in_ycenter) then $
ycenter = in_ycenter[extlist[ext]-1] $
else begin
print,'ASTROM: ycenter does not have enough entries, no value present for'
print,' extension',strcompress(extlist[ext]),' unable to continue.'
endelse
endelse
IF hdrinfo.ra ge 0. THEN BEGIN
IF hdrinfo.epoch NE 2000.0 THEN BEGIN
precess,hdrinfo.ra,hdrinfo.dec,hdrinfo.epoch,2000.0,/radian
print,' *** Header coordinates precessed to J2000 from ', $
strtrim(string(hdrinfo.epoch,format='(f10.2)'),2)
hdrinfo.epoch=2000.0
ENDIF
ENDIF
IF starfile eq '[[DEFAULT]]' THEN BEGIN
starfile = objname[0] + exttag + '.cat'
ENDIF
IF plastfile eq '[[DEFAULT]]' THEN BEGIN
plastfile = objname[0] + exttag + '.pla'
ENDIF
; See if PLAST file exists and set a flag for future use.
doplast = exists(plastfile) and plastfile ne 'none'
; Read in astrometry information saved from last run. If not found,
; ask for information.
if exists(astinf) then begin
version=''
openr,lun,astinf,/get_lun
readf,lun,version,format='(a)'
if version ne 'ASTROM v1.0' then begin
print,'Illegal astrom.inf file, version tag is wrong.'
return
endif
readf,lun,raoff,decoff
readf,lun,pscale0
readf,lun,rang0
xflip=1
yflip=1
readf,lun,xflip,yflip,format='(i2,1x,i2)'
free_lun,lun
endif else begin
raoff = 0.0
decoff = 0.0
read,pscale0,prompt='Initial plate scale guess (arcsec/pixel) '
read,rang0,prompt='Rotation angle (degrees) '
read,ans,prompt='Is the X-axis flipped with respect to the sky? '
if strmid(ans,0,1) eq 'y' then xflip=-1 else xflip=1
read,ans,prompt='Is the Y-axis flipped with respect to the sky? '
if strmid(ans,0,1) eq 'y' then yflip=-1 else yflip=1
openw,lun,astinf,/get_lun
printf,lun,'ASTROM v1.0'
printf,lun,raoff,decoff
printf,lun,pscale0
printf,lun,rang0
printf,lun,xflip,yflip,format='(i2,1x,i2)'
free_lun,lun
in_trustcenter=0
endelse
; If starcat or plast file is needed, get a center to use for the
; extraction.
IF not exists(starfile) or newcat or (not doplast and plastfile ne 'none') THEN BEGIN
; First try to get the location of the object, this assumes that the object
; was the intended target and the pointing was set to this object. This
; protects against a bad header value. If object name does not yield a
; valid position then we'll have to resort to other means to get the center.
IF do_objects and center eq 0 THEN BEGIN
ephem,hdrinfo.jd,688,2+50,'e'+objname[0],eph
objra=eph[0]
objdec=eph[1]
rastr,objra,1,ras
decstr,objdec,0,decs
print,'Object location: ',ras,' ',decs
ENDIF ELSE BEGIN
objra = -99.9
objdec = -99.0
ENDELSE
IF objra lt 0.0 and center eq 0 THEN BEGIN
rastr,hdrinfo.ra,1,ras
decstr,hdrinfo.dec,0,decs
print,' Header indicates ra,dec = ',ras,' ',decs,' @ ',jds
rastr,hdrinfo.ra+raoff*cos(hdrinfo.dec),1,ras
decstr,hdrinfo.dec+decoff,0,decs
print,' with offset -- ra,dec = ',ras,' ',decs,' @ ',jds
read,ans,prompt='Use raw header (R), with offset (O), or manual (default)'
ans=strlowcase(strmid(ans,0,1))
IF ans eq 'r' THEN BEGIN
objra = hdrinfo.ra
objdec = hdrinfo.dec
ENDIF ELSE IF ans eq 'o' THEN BEGIN
objra = hdrinfo.ra+raoff*cos(hdrinfo.dec)
objdec = hdrinfo.dec+decoff
ENDIF ELSE BEGIN
read,ans,prompt='RA of image center (J2000) ? ',format='(a)'
objra=raparse(strcompress(ans,/remove_all))
read,ans,prompt='Dec of image center (J2000) ? ',format='(a)'
objdec=decparse(strcompress(ans,/remove_all))
ENDELSE
ENDIF ELSE IF center eq 0 THEN BEGIN
print,'Using nominal object position for star catalog extraction.'
ENDIF ELSE IF center eq 1 THEN BEGIN
print,'Using raw header for star catalog extraction.'
objra = hdrinfo.ra
objdec = hdrinfo.dec
ENDIF ELSE BEGIN
print,'Using offset from header for star catalog extraction.'
objra = hdrinfo.ra+raoff*cos(hdrinfo.dec)
objdec = hdrinfo.dec+decoff
ENDELSE
psize=fix(sqrt(sz[1]^2+sz[2]^2)*pscale0*1.05+0.5)
ENDIF
IF not exists(starfile) or newcat THEN BEGIN
killrefstars=1
print,' --> running refnet....'
refnet,objra,objdec,psize,psize,maglim,maglim,starfile,CATPATH=catpath
ENDIF ELSE BEGIN
killrefstars=0
ENDELSE
IF not exists(starfile) THEN BEGIN
print,'ASTROM: The star catalog file '+starfile+ $
' does not exist. Aborting'
return
ENDIF
IF not doplast and plastfile ne 'none' THEN BEGIN
print,' --> running plast....'
plast,hdrinfo.jd,objra,objdec,psize,psize,plastfile,title=imfile
doplast = exists(plastfile)
ENDIF
;Setup display image and coordinate system on display.
dispim = 0.
dispim = rebin(raw,sz[1]/binfac,sz[2]/binfac,sample=0)
setwin,window,xsize=sz[1]/binfac,ysize=sz[2]/binfac
skysclim,dispim,lowval,hival,meanval,sigma
lowval = max([meanval-3*sigma,min(dispim)])
hival = min([meanval+8*sigma,max(dispim)])
savec=5
bim=0B
bim = bytscl(dispim,min=lowval,max=hival,top=!d.n_colors-1-savec)+savec
loadct,0,bottom=savec,/silent
tvlct,red,gre,blu,/get
; 0 black
; 1 green
; 2 reddish
; 3 yellow
; 4 red
red[0]=0
gre[0]=0
blu[0]=0
red[1]=0
gre[1]=255
blu[1]=0
red[2]=200
gre[2]=50
blu[2]=80
red[3]=255
gre[3]=255
blu[3]=0
red[4]=255
gre[4]=0
blu[4]=0
tvlct,red,gre,blu
if fullrefresh then astrom_dim,window,bim,binfac,border
if not clean or preclean then begin
image=raw
endif else begin
print,' --> cleaning cosmic rays from image....'
acre,raw,image,5,4
if saveclean then begin
savename=root+'c.'+suffix
writefits,path+savename,image,hdr
endif
endelse
dispim = rebin(image,sz[1]/binfac,sz[2]/binfac,sample=0)
bim = bytscl(dispim,min=lowval,max=hival,top=!d.n_colors-1-savec)+savec
if fullrefresh then astrom_dim,window,bim,binfac,border
; Setup a structure with commonly used information
dv = { $
nx: sz[1], $
ny: sz[2], $
binfac: binfac, $
win: window, $
bim: bim, $
border: border, $
objrad: objrad, $
gain: gain, $
image: image, $
exptime: exptime }
; Get plast information and collect the information for plotting
IF doplast THEN BEGIN
openr,lun,plastfile,/get_lun
line=''
FOR i=1,4 DO readf,lun,line,format='(a1)'
numobj=0
astname0=''
WHILE not eof(lun) DO BEGIN
readf,lun,astname0,mag0,format='(a23,f5.2)'
readf,lun,line,format='(a1)'
IF numobj eq 0 THEN BEGIN
pastname=astname0
pastmag=mag0
numobj=1
ENDIF ELSE BEGIN
pastname=[pastname,astname0]
pastmag=[pastmag,mag0]
numobj=numobj+1
ENDELSE
ENDWHILE
IF numobj eq 0 THEN doplast=0
free_lun,lun
ENDIF
IF doplast THEN BEGIN
astcode=strarr(numobj)
FOR i=0,numobj-1 DO BEGIN
tmpstr=str_sep(strtrim(pastname[i],2),' ')
rawname=tmpstr[0]
IF strpos(rawname,'P/') ne -1 THEN BEGIN
astcode[i]='c'+rawname
ENDIF ELSE IF strpos(rawname,'C/') ne -1 THEN BEGIN
astcode[i]='c'+rawname
ENDIF ELSE BEGIN
astcode[i]='e'+rawname
ENDELSE
ENDFOR
ephem,replicate(hdrinfo.jd,numobj),688,23+50,astcode,eph
astra=eph[0,*]
astdec=eph[1,*]
astvra=eph[18,*]
astvdec=eph[19,*]
asterr=eph[20,*]
ephem,replicate(hdrinfo.jd+1.0/24.0d0,numobj),688,2+50,astcode,eph
astrap1=eph[0,*]
astdecp1=eph[1,*]
astrate=sqrt((astrap1-astra)^2+(astdecp1-astdec)^2)*!radeg*3600.0
astra=astra[*]
astdec=astdec[*]
astvra=astvra[*]
astvdec=astvdec[*]
asterr=asterr[*]
astrap1=astrap1[*]
astdecp1=astdecp1[*]
c1=sqrt((astvra*cos(astdec))^2+astvdec^2)
astvra=astvra*cos(astdec)/c1
astvdec=astvdec/c1
asterr=asterr/3600.0d0*!dpi/180.0d0
ENDIF
;===========================================
; This section is run for each sub-exposure.
FOR ie=0,n_elements(subexp)-1 DO BEGIN
IF subexp[ie] eq '[[DEFAULT]]' THEN BEGIN
suff = suffix+exttag
ENDIF ELSE BEGIN
suff = suffix+exttag+'_'+subexp[ie]
IF ie ne 0 THEN print,''
print,'Sub-exposure ',subexp[ie]
ENDELSE
trustcenter=in_trustcenter
reffile='Refstars/'+suff+'.ref'
fnsrc = imfile+'.src'+exttag
thisexp=root+'.'+suff
IF exists(reffile) and killrefstars THEN spawn,'rm -f '+reffile
; Look for image center from previous reduction.
thisfile=-1
last=''
found_center=0
IF exists('centers.dat') THEN BEGIN
openr,lun,'centers.dat',/get_lun
line=''
done=0
WHILE not eof(lun) and not done DO BEGIN
readf,lun,line,format='(a)'
cline=strcompress(strtrim(line,2))
if cline ne '' then begin
words=str_sep(cline,' ')
fnf=str_sep(words[0],'.')
if n_elements(words) lt 3 or n_elements(fnf) lt 2 then begin
print,'ASTROM: Fatal error occurred reading the centers.dat file. The last line'
print,'read is:'
print,'[',line,']'
print,'This should have three blank separated fields, the first would be a file'
print,'name of the form, root.suffix, and the next two fields are the RA and Dec'
print,'of the center. You must fix this file before astrom can continue.'
print,''
free_lun,lun
return
endif
IF fnf[1] gt suff THEN BEGIN
done=1
ENDIF ELSE BEGIN
IF thisexp eq words[0] THEN BEGIN
cra=raparse(words[1])
cdec=decparse(words[2])
found_center=1
done=1
ENDIF ELSE last=line
ENDELSE
endif
ENDWHILE
free_lun,lun
ENDIF
headercoord=0
IF not found_center THEN BEGIN
hra=hdrinfo.ra
hdec=hdrinfo.dec
IF hra ge 0. THEN BEGIN
headercoord=1
ENDIF ELSE BEGIN
read,ans,prompt='RA of image center (J2000) ? ',format='(a)'
hra=raparse(ans)
read,ans,prompt='Dec of image center (J2000) ? ',format='(a)'
hdec=decparse(ans)
headercoord=0
trustcenter=0
cra=hra
cdec=hdec
ENDELSE
rastr,hra,1,hras
decstr,hdec,0,hdecs
ENDIF
; display image
if fullrefresh then astrom_dim,window,bim,binfac,border
first=1
IF headercoord THEN BEGIN
cra=hra+raoff*cos(hdec)
cdec=hdec+decoff
print,' *** Using header coordinates for image center.'
ENDIF
info = { $
racen: cra, $ ; RA of center of image
deccen: cdec, $ ; Dec center of image
pscale: pscale0, $ ; Scale of image, arcsec/pixel
rang: rang0/!radeg, $ ; Rotation angle of image (radians).
xflip: xflip, $ ; 1 = no flip, -1 = X-axis flipped
yflip: yflip, $ ; 1 = no flip, -1 = Y-axis flipped
xc: sz[1]/2.0, $ ; X center of image
yc: sz[2]/2.0, $ ; Y center of image
raref: cra, $ ; RA of xi,eta reference point.
decref: cdec, $ ; Dec of xi,eta reference point.
xcref: sz[1]/2.0, $ ; X reference location
ycref: sz[2]/2.0, $ ; Y reference location
ra0: cra, $ ; RA of optical axis (might not be cra)
dec0: cdec, $ ; Dec of optical axis (might not be cdec)
xc0: xcenter, $ ; X location of optical axis (might not be xc)
yc0: ycenter } ; Y location of optical axis (might not be yc)
astxy2rd,info.xc0,info.yc0,info,ra0,dec0
info.ra0 = ra0
info.dec0 = dec0
; Mark reminder spot
if first and do_objects and not noremind then begin
print,' Click a spot to remember'
cursor,xloc,yloc,3
if !mouse.button eq 4 then return
first=0
endif
if do_objects and not noremind then begin
oplot,[-20,20,20,-20,-20]+xloc,[-20,-20,20,20,-20]+yloc,color=4
empty ; flush graphics buffer output
endif
IF doplast THEN BEGIN
if numobj gt 0 then $
print,'Asteroids found in plast support file that may be on image:'
FOR i=0,numobj-1 DO BEGIN
rastr,astra[i],1,str1
decstr,astdec[i],0,str2
print,i,pastname[i]+blanks,str1,str2,asterr[i]*!radeg*3600.0,'"', $
astrate[i],'"/hr V=',pastmag[i], $
format='(i3,1x,a12,1x,a,1x,a,1x,f7.1,a,2x,f4.1,a,f4.1)'
ENDFOR
ENDIF
; At this point, we need to establish the coordinate transformation
; from pixels to ra,dec. There are three ways to get there. All
; three methods require have a list of (x,y) positions for a set of
; catalog stars. (1) interactively identify the catalog <-> image
; matchup, measure the catalog stars. (2) read an external list of
; sources and match to the source catalog. (3) read the correspondence
; list that was saved from the last run. Once the list is in place,
; then we proceed to fit for the transformation. In order, first #3
; then #2, then method #1 is tried.
; Try to load the star measurements from the last run.
IF exists(reffile) THEN BEGIN
openr,lref,reffile,/get_lun
sv_objrad=0.0
nstars=0L
readu,lref,sv_objrad,nstars
IF abs(sv_objrad-objrad) gt 0.1 THEN BEGIN
print,'The object aperture radius in the saved star data is ', $
sv_objrad
print,' You have requested an object radius of ',objrad
print,' To continue, the saved star data will be deleted and '+ $
'will be regenerated'
print,' with the new radius. If this is not what you want to'+ $
' do, then you should'
print,' quit and run again with OBJRAD=',sv_objrad
print,''
read,ans,prompt='Do you want to delete and continue (y/n) ?'
if strmid(ans,0,1) ne 'y' then return
free_lun,lref
spawn,'rm '+reffile
ENDIF ELSE BEGIN
xpos = fltarr(nstars)
ypos = fltarr(nstars)
fwhm = fltarr(nstars)
mag = fltarr(nstars)
err = fltarr(nstars)
sig_max = fltarr(nstars)
sra = dblarr(nstars)
sdec = dblarr(nstars)
smag = fltarr(nstars)
bad = intarr(nstars)
readu,lref,xpos,ypos,fwhm,mag,err,sig_max,sra,sdec,smag
z = where(bad ne 1 and sig_max le maxphotsig)
free_lun,lref
ENDELSE
ENDIF
; If no old stars, we'll need to read in the raw catalog list.
IF not exists(reffile) THEN BEGIN
; Read in list of stars (this is the REFNET output format).
print,' --> reading star catalog....'
readcol,starfile,hr,m1,s1,dgas,m2,s2,smag, $
format='d,d,d,a,d,d,f',/silent
signas = strmid(dgas,0,1)
dg = fix(strmid(dgas,1,2))
hmstorad,hr,m1,s1,sra
sign = replicate(1.0,n_elements(dg))
z=where(signas eq '-',count)
if count ne 0 then sign[z] = -1.0
dmstorad,sign,abs(dg),m2,s2,sdec
if maglim eq 30.0 then $
ssz = (20.0-smag+1)/2 + 0.5 $
else $
ssz = (maglim-smag+1)/2 + 0.5
z=where(ssz gt 4,count)
if count ne 0 then ssz[z]=3.0
z=where(ssz lt 0.5,count)
if count ne 0 then ssz[z]=0.5
bad=intarr(n_elements(sra))
z=where(smag lt maglim,count)
IF count eq 0 THEN BEGIN
print,'All catalog stars are fainter than the limit,',maglim
return
ENDIF
sra=sra[z]
sdec=sdec[z]
smag=smag[z]
ssz=ssz[z]
bad=bad[z]
; Force xi,eta reference point to be center of image
info.raref = info.racen
info.decref = info.deccen
info.xcref = info.xc
info.ycref = info.yc
astrd2xy,sra,sdec,info,x,y,XI=sxi,ETA=seta
ENDIF
; If no old stars, and there is a source list, use it.
IF not exists(reffile) and exists(fnsrc) THEN BEGIN
print,' --> reading source list....'
list=readfits(fnsrc,hdrsrc,/silent)
gap=sxpar(hdrsrc,'GAP')
IF gap ne objrad THEN $
print,'Warning, GAP <> OBJRAD. gap=',gap,' objrad=',objrad
nlist=n_elements(list)/5
xpos=reform(list[*,0],nlist)
ypos=reform(list[*,1],nlist)
fwhm=reform(list[*,2],nlist)
mag =reform(list[*,3],nlist)
err =reform(list[*,4],nlist)
sig_max=reform(raw[fix(xpos+0.5),fix(ypos+0.5)],nlist)
zz=where(sig_max lt maxphotsig,nlist0)
IF nlist0 ne nlist THEN BEGIN
xpos=xpos[zz]
ypos=ypos[zz]
fwhm=fwhm[zz]
mag =mag[zz]
err = err[zz]
sig_max=sig_max[zz]
nlist=nlist0
ENDIF
nlist0 = min([1.3*n_elements(x),n_elements(xpos)])
IF nlist0 ne nlist THEN BEGIN
zz = sort(mag)
zz = zz[0:nlist0-1]
xpos=xpos[zz]
ypos=ypos[zz]
fwhm=fwhm[zz]
mag =mag[zz]
err = err[zz]
sig_max=sig_max[zz]
nlist=nlist0
ENDIF
; This block handles automatically cross-referencing the catalog and source
; list. The two-star solution must be pretty good, especially in rotation
; for this to work.
if not twostar then begin
astrom_dim,window,bim,binfac,border
oplot,xpos,ypos,psym=5,color=1,symsize=2.0
print,' --> determine source/catalog position offset....'
; dxdy=replicate(2,sz[1],sz[2])
dxdy=intarr(sz[1],sz[2])
FOR i=0,n_elements(x)-1 DO BEGIN
dx=fix(xpos-x[i]+0.5+sz[1]/2.0)
dy=fix(ypos-y[i]+0.5+sz[2]/2.0)
zd=where(dx ge 0 and dx lt sz[1] and $
dy ge 0 and dy lt sz[2],countzd)
IF countzd gt 0 THEN BEGIN
dxdy[dx[zd],dy[zd]]=dxdy[dx[zd],dy[zd]]+1
ENDIF
ENDFOR
dxdy2=rebin(dxdy,sz[1]/2,sz[2]/2)
zd=where(dxdy2 eq max(dxdy2))
xoff=((zd[0] mod (sz[1]/2))-(sz[1]/2)/2.0)*2.0
yoff=(zd[0]/(sz[1]/2)-(sz[2]/2)/2.0)*2.0
print,xoff,yoff
; zd=where(dxdy eq max(dxdy))
; xoff=(zd[0] mod sz[1])-sz[1]/2.0
; yoff=zd[0]/sz[1]-sz[2]/2.0
print,'First pass offset ',xoff,yoff,sz[1],sz[2]
fndrad=12.0
basphote,1,dxdy,1.,xoff+sz[1]/2.0,yoff+sz[2]/2.0,fndrad, $
fndrad+10,fndrad+40,/nolog,/silent,xcen=nxc,ycen=nyc
nxc = nxc - sz[1]/2.0
nyc = nyc - sz[2]/2.0
print,'Improved? offset ',nxc,nyc
endif else begin
nxc = 0
nyc = 0
xoff = 0
yoff = 0
endelse
;!! The cross-correlation doesn't quite get the two-star solution
;!! right, so that regions of the reference stars (like in the corner)
;!! are missed in large numbers. This is probably only important for
;!! images with appreciable non-linearities. (eg., LONEOS).
; If the centroid of the correlation spot doesn't agree with
; the location of the max, then it is likely the rotation angle
; isn't close enough. It's probably safer to do the interactive
; determination at this stage.
IF twostar or abs(nxc-xoff) gt drthresh/2.0 or $
abs(nyc-yoff) gt drthresh/2.0 THEN BEGIN
if not twostar then $
print,bel,'Warning! auto correlation was not good, try manual.'
astrom_twostar,dv,info,sra,sdec,x,y,ssz,valid
if not valid then return
; Redetermine ra and dec offset for "good" header coordinates.
IF headercoord THEN BEGIN
raoff=(info.racen-hra)/cos(info.deccen)
decoff=info.deccen-hdec
ENDIF
ENDIF ELSE BEGIN
x = x + xoff
y = y + yoff
ENDELSE
print,' --> cross-link source list and catalog....'
sidx=replicate(-1L,n_elements(x))
astrom_dim,window,bim,binfac,border
for i=0,n_elements(x)-1 do $
oplot,[x[i]],[y[i]],psym=4,color=2,symsize=ssz[i]
oplot,xpos,ypos,psym=5,color=1,symsize=0.8
FOR i=0,n_elements(x)-1 DO BEGIN
dr = (xpos-x[i])^2 + (ypos-y[i])^2
zt = where(dr eq min(dr) and dr lt drthresh^2)
sidx[i] = zt[0]
IF zt[0] ne -1 THEN $
oplot,[xpos[zt[0]]],[ypos[zt[0]]],psym=4,color=3,symsize=2.2 $
else begin
zt = where(dr eq min(dr) and dr lt (drthresh^2)*4)
if zt[0] ne -1 then begin
oplot,[xpos[zt[0]]],[ypos[zt[0]]],psym=4,color=3,symsize=5
endif
endelse
ENDFOR
zy = where(sidx ne -1,nstars)
if zy[0] eq -1 then begin
print,'Fatal error! no stars lined up. quitting on this image.'
return
endif
x = x[zy]
y = y[zy]
ssz = ssz[zy]
xpos = xpos[sidx[zy]]
ypos = ypos[sidx[zy]]
fwhm = fwhm[sidx[zy]]
mag = mag[sidx[zy]]
err = err[sidx[zy]]
sig_max = sig_max[sidx[zy]]
sra = sra[zy]
sdec = sdec[zy]
smag = smag[zy]
bad = intarr(nstars)
oplot,xpos,ypos,psym=4,color=3,symsize=3.2
z = where(bad ne 1)
; Compute a coverage statistic
xq = fix(3*xpos/float(dv.nx))
yq = fix(3*ypos/float(dv.ny))
qu = xq + yq*3
qhist = histogram(qu,min=0,max=8)
qmed = median(qhist)
zq=where(qhist lt qmed/2 or qhist eq 0,countzq)
if countzq ne 0 then $
print,bel,'Warning! poor star coverage on', $
strcompress(countzq),' enneanants: ',imfile,' ',exttag
if countzq gt 2 then begin
if countzq le 4 then begin
read,ans,prompt='Do you want to continue (default=n)? '
if strmid(ans,0,1) ne 'y' then return
endif else begin
print,'Very bad coverage, quitting on this image.'
return
endelse
endif
; Save the catalog star measurements to a binary file.
astrom_svref,reffile,gap,nstars,xpos,ypos,fwhm, $
mag,err,sig_max,sra,sdec,smag
ENDIF
; need to regenerate some of the info structure information
; for the sn2xy transformation if using source lists.
IF exists(reffile) or exists(fnsrc) THEN BEGIN
; first select two stars to setup a good two-star solution
idx = sort(xpos[z]^2+ypos[z]^2)
loptr=fix(0.1*n_elements(idx))
hiptr=fix(0.9*n_elements(idx))
;print,loptr,hiptr,n_elements(idx)
il = max([(loptr-5),0]) - loptr
ih = min([(hiptr+5),n_elements(idx)-1]) - hiptr
;print,il,ih
;return
angs=replicate(1000.0,ih-il+1)
for ii=il,ih do begin
zs1 = z[idx[loptr+ii]]
zs2 = z[idx[hiptr+ii]]
if zs1 ge 0 and zs1 lt n_elements(xpos) and $
zs2 ge 0 and zs2 lt n_elements(xpos) then begin
astrom_regen,xpos,ypos,sra,sdec,z,zs1,zs2,info
angs[ii-il]=info.rang
endif
; astrom_svctr,thisexp,info
endfor
gangs=angs[where(angs lt 500.0)]
zit=where(angs eq median(gangs))
; print,'ASTROM: best at ',zit[0]
zs1 = z[idx[fix(0.1*n_elements(idx))+zit[0]+il]]
zs2 = z[idx[fix(0.9*n_elements(idx))+zit[0]+il]]
astrom_regen,xpos,ypos,sra,sdec,z,zs1,zs2,info
astrom_svctr,thisexp,info
ENDIF
IF not exists(reffile) and not exists(fnsrc) THEN BEGIN
retry:
; Force xi,eta reference point to be center of image
info.raref = info.racen
info.decref = info.deccen
info.xcref = info.xc
info.ycref = info.yc
astrd2xy,sra,sdec,info,x,y,XI=sxi,ETA=seta
z = where(x gt border[0] and x lt sz[1]-border[1] and $
y gt border[3] and y lt sz[2]-border[2],count)
rastr,info.racen,1,ras1
decstr,info.deccen,0,decs1
print,' Center used: ',ras1,' ',decs1
if count eq 0 and last ne '' then begin
print,'??? No stars found, reloading last known center.'
print,last
words=str_sep(last,' ')
cra=raparse(words[1])
cdec=decparse(words[2])
last=''
info.racen=cra
info.deccen=cdec
headercoord=0
goto,retry
endif else if count eq 0 then begin
rastr,hdrinfo.ra,1,ras
decstr,hdrinfo.dec,0,decs
print,' Header indicates ra,dec = ',ras,' ',decs,' @ ',jds
print,'No stars on frame.'
read,ans,prompt='Do you want to (re)enter a corrected center? '
if strmid(ans,0,1) eq 'n' then return
read,ans,prompt='RA of image center? ',format='(a)'
hra=raparse(ans)
read,ans,prompt='Dec of image center? ',format='(a)'
hdec=decparse(ans)
info.racen=hra
info.deccen=hdec
headercoord=0
trustcenter=0
goto,retry
endif
; Do the catalog/field correlation for non-trusted center
IF not trustcenter THEN BEGIN
astrom_twostar,dv,info,sra,sdec,x,y,ssz,valid
IF not valid THEN return
; Redetermine ra and dec offset for "good" header coordinates.
IF headercoord THEN BEGIN
raoff=(info.racen-hra)/cos(info.deccen)
decoff=info.deccen-hdec
ENDIF
ENDIF ; end of catalog/field correlation for non-trusted center
astrom_svctr,thisexp,info
if fullrefresh then begin
astrom_dim,window,bim,binfac,border
if do_objects and not noremind then $
oplot,[-20,20,20,-20,-20]+xloc,[-20,-20,20,20,-20]+yloc,color=4
endif
astrd2xy,sra,sdec,info,x,y
for i=0,n_elements(x)-1 do $
oplot,[x[i]],[y[i]],psym=4,color=2,symsize=ssz[i]
; Select those stars on the frame
nstars=n_elements(sra)
xpos=fltarr(nstars)
ypos=xpos
perr=ypos+9999.0
mag=perr
err=perr
fwhm=perr
sig_max=fltarr(nstars)
z = where(x le border[0] or x ge sz[1]-border[1] or $
y le border[3] or y ge sz[2]-border[2],count)
IF count ne 0 THEN bad[z]=1
z = where(x gt border[0] and x lt sz[1]-border[1] and $
y gt border[3] and y lt sz[2]-border[2] and bad eq 0,count)
if count ne 0 then begin
IF maxstars gt 0 THEN BEGIN
zz=z
numtodo=maxstars
ENDIF ELSE numtodo=count
IF numtodo gt count THEN numtodo=count
for i=0,numtodo-1 do begin
IF numtodo eq count THEN BEGIN
iz=z[i]
ENDIF ELSE BEGIN
i1=0
i4=n_elements(zz)-1
is=fix(randomu(seed)*(i4+1))
i2=is-1
i3=is+1
iz=zz[is]
IF i eq 0 THEN BEGIN
z=iz
ENDIF ELSE BEGIN
z=[z,iz]
ENDELSE
IF i1 gt i2 THEN BEGIN
zz=zz[i3:i4]
ENDIF ELSE IF i3 gt i4 THEN BEGIN
zz=zz[i1:i2]
ENDIF ELSE BEGIN
zz=[zz[i1:i2],zz[i3:i4]]
ENDELSE
ENDELSE
basphote,gain,image,exptime,x[iz],y[iz], $
objrad,objrad+10,objrad+30,/nolog,/silent, $
xcen=xm,ycen=ym,fwhm=fwhm0,mag=mag0,err=err0, $
max=max0,boxmrad=max([5,objrad])
xpos[iz]=xm
ypos[iz]=ym
perr[iz]= sqrt((xm-x[iz])^2 + (ym-y[iz])^2)*info.pscale
fwhm[iz]=fwhm0
mag[iz]=mag0
err[iz]=err0
sig_max[iz]=max0
rastr,sra[iz],1,ras
decstr,sdec[iz],0,decs
setusym,-1
oplot,[xpos[iz]],[ypos[iz]],psym=8,color=3,symsize=2.5
setusym,1
empty ; flush graphics buffer output
endfor
IF numtodo ne count THEN count=numtodo
oplot,x[z],y[z],psym=4,color=2
endif else begin
print,'No stars on frame, unable to continue.'
return
endelse
; Save the catalog star measurements to a binary file.
astrom_svref,reffile,objrad,count,xpos[z],ypos[z],fwhm[z], $
mag[z],err[z],sig_max[z],sra[z],sdec[z],smag[z]
ENDIF ; end block for measuring star network
IF doplast THEN BEGIN
astrd2xy,astra,astdec,info,astx,asty
astrd2xy,astra+astvra*asterr/cos(astdec),astdec+astvdec*asterr, $
info,astvxp,astvyp
astrd2xy,astra-astvra*asterr/cos(astdec),astdec-astvdec*asterr, $
info,astvxm,astvym
astrd2xy,astrap1,astdecp1,info,astxp1,astyp1
astxvar=[[astvxm],[astx],[astvxp]]
astyvar=[[astvym],[asty],[astvyp]]
ENDIF
; Convert from ra,dec to xi,eta (Smart, p283)
astrd2sn,sra,sdec,info.ra0,info.dec0,xi,eta
xi = xi*180.0d0/!dpi*3600.0d0
eta = eta*180.0d0/!dpi*3600.0d0
robomean,fwhm[z],3.0,0.5,meanfwhm,avgdev,fwhmstddev
print,'Seeing --> (FWHM) ',meanfwhm,' +/- ', $
fwhmstddev,meanfwhm*info.pscale, $
format='(a,f4.1,a,f3.1," pixels ",f4.1," arcsec")'
; Exclude objects with discrepant sizes.
zz=where(abs(fwhm - meanfwhm) gt 3.0*fwhmstddev or fwhm le 0.7,countzz)
IF countzz ne 0 THEN bad[zz]=1
; Exclude excessivly faint objects.
zz=where(abs(mag gt 80.0) or err gt 0.1,countzz)
IF countzz ne 0 THEN bad[zz]=1
renormfac=sqrt(sz[1]^2+sz[2]^2)
dx = (xpos-info.xc0)/renormfac
dy = (ypos-info.yc0)/renormfac
astsolve,dx,dy,xi,eta,xiterms,etaterms,renormfac,bad,cxi,ceta, $
edit=edit,xflip=info.xflip ne 1,yflip=info.yflip ne 1
zg = where(bad eq 0,countgood)
dx = (xpos[zg]-info.xc0)/renormfac
dy = (ypos[zg]-info.yc0)/renormfac
xfit = asteval(dx,dy,cxi,xiterms)
efit = asteval(dx,dy,ceta,etaterms)
setwin,10,xsize=400,ysize=800
!p.multi=[0,1,7]
ytitle='Xi resid (arcsec)'
plot,dx*renormfac,xi[zg]-xfit,psym=7,symsize=0.5,xtitle='dx (pixels)',ytitle=ytitle
plot,dy*renormfac,xi[zg]-xfit,psym=7,symsize=0.5,xtitle='dy (pixels)',ytitle=ytitle
plot,sqrt(dx^2+dy^2)*renormfac,xi[zg]-xfit,psym=7,symsize=0.5,xtitle='r (pixels)',ytitle=ytitle
plot,(dx*renormfac)^2,xi[zg]-xfit,psym=7,symsize=0.5,xtitle='dx*dx (pixels)',ytitle=ytitle
plot,(dy*renormfac)^2,xi[zg]-xfit,psym=7,symsize=0.5,xtitle='dy*dy (pixels)',ytitle=ytitle
plot,mag[zg],xi[zg]-xfit,psym=7,symsize=0.5,xtitle='mag',ytitle=ytitle
plot,fwhm[zg],xi[zg]-xfit,psym=7,symsize=0.5,xtitle='fwhm (pixels)',ytitle=ytitle
setwin,11,xsize=400,ysize=800
plot,dx*renormfac,eta[zg]-efit,psym=7,symsize=0.5
plot,dy*renormfac,eta[zg]-efit,psym=7,symsize=0.5
plot,sqrt(dx^2+dy^2)*renormfac,eta[zg]-efit,psym=7,symsize=0.5
plot,(dx*renormfac)^2,eta[zg]-efit,psym=7,symsize=0.5
plot,(dy*renormfac)^2,eta[zg]-efit,psym=7,symsize=0.5
plot,mag[zg],eta[zg]-efit,psym=7,symsize=0.5
plot,fwhm[zg],eta[zg]-efit,psym=7,symsize=0.5
!p.multi=0
IF fullrefresh THEN BEGIN
setwin,14
plot,xpos[zg],ypos[zg],psym=8,symsize=0.5, $
xtitle='X (pixels)',ytitle='Y (pixels)', $
title='Xi (x-dir) and Eta (y-dir) error plot'
sfac=200.0/max([xi[zg]-xfit,eta[zg]-efit])
FOR i=0,n_elements(zg)-1 DO BEGIN
oplot,xpos[zg[i]]+sfac*[0.,xi[zg[i]]-xfit[i]], $
ypos[zg[i]]+sfac*[0.,eta[zg[i]]-efit[i]]
ENDFOR
ENDIF
;read,ans,prompt='continue -->'
astrom_dim,window,bim,binfac,border
IF do_objects and not noremind THEN $
oplot,[-20,20,20,-20,-20]+xloc,[-20,-20,20,20,-20]+yloc,color=4
IF doplast THEN BEGIN
oplot,astx,asty,psym=8,color=1,symsize=0.2
FOR i=0,numobj-1 DO BEGIN
xyouts,astx[i],asty[i],' '+strtrim(string(i),2),color=1
IF asterr[i]*!radeg*3600.0/info.pscale gt 5.0 THEN BEGIN
oplot,astxvar[i,*],astyvar[i,*],color=1
oplot,[astx[i],astxp1[i]],[asty[i],astyp1[i]],color=3
ENDIF ELSE BEGIN
oplot,[astx[i]],[asty[i]],psym=4,color=1,symsize=2.0
ENDELSE
ENDFOR
ENDIF
setusym,-1
FOR i=0,nstars-1 DO BEGIN
IF bad[i] eq 0 THEN BEGIN
oplot,[xpos[i]],[ypos[i]],psym=8,color=3,symsize=2.5
ENDIF ELSE BEGIN
oplot,[xpos[i]],[ypos[i]],psym=8,color=2,symsize=2.5
ENDELSE
ENDFOR
zok=where((sig_max le maxphotsig) and (bad eq 0) and $
(err lt 0.15),countok)
IF countok gt 0 THEN BEGIN
robomean,smag[zok]-mag[zok],3.0,0.5,magdiff,avgdev,stddev
zok1=where(sig_max le maxphotsig and bad eq 0 and $
abs(smag-mag-magdiff) lt 2.0*stddev,countok1)
meanerr,smag[zok1]-mag[zok1],err[zok1],magdiff,avgdev,stddev
print,'Mean catalog magnitude difference ',magdiff, $
stddev/sqrt(countok1-1),'from',countok1,'stars', $
format='(a,f7.3," +/- ",f6.3,1x,a,1x,i3,1x,a)'
fwhmrange=minmax(fwhm[zok])
zok=where((sig_max le maxphotsig) and (err lt 0.15) and $
fwhm ge fwhmrange[0] and fwhm le fwhmrange[1],countok)
ENDIF ELSE BEGIN
magdiff = 0.0
print,'Warning: no stars in valid photometric range, magdiff set to 0.'
print,'Range in peak counts: ',minmax(sig_max[where(bad eq 0)])
ENDELSE
; Save 2-star solution for next run.
openw,lun,astinf,/get_lun
printf,lun,'ASTROM v1.0'
printf,lun,raoff,decoff
printf,lun,info.pscale
printf,lun,info.rang*!radeg
printf,lun,info.xflip,info.yflip,format='(i2,1x,i2)'
free_lun,lun
; Save xi and eta coefficients to the master file.
astprmt ; promotes file type version of fitcoeff.dat
rastr,info.ra0,4,ras
decstr,info.dec0,3,decs
cinfo=string(info.xc0,info.yc0,format='(2(1x,e14.7))') + $
' '+ras+' '+decs+ $
string(xiterms,format='(10(1x,i1))') + $
string(cxi,format='(10(1x,e19.12))')
repwrite,'fitcoeff.dat',thisexp+' xi',thisexp+' xi'+cinfo, $
header='ASTFIT v1.0'
cinfo=string(info.xc0,info.yc0,format='(2(1x,e14.7))') + $
' '+ras+' '+decs+ $
string(etaterms,format='(10(1x,i1))') + $
string(ceta,format='(10(1x,e19.12))')
repwrite,'fitcoeff.dat',thisexp+' eta',thisexp+' eta'+cinfo, $
header='ASTFIT v1.0'
IF photstars and countok gt 0 THEN BEGIN
if fullrefresh then begin
astrom_dim,window,bim,binfac,border
setusym,-1
FOR i=0,countok-1 DO BEGIN
oplot,[xpos[zok[i]]],[ypos[zok[i]]],psym=8,color=3,symsize=2.5
xyouts,xpos[zok[i]],ypos[zok[i]],' '+strtrim(string(i),2),color=3
ENDFOR
setusym,1
endif
hardim,bytscl(image,min=lowval,max=hival,top=255),0,255, $
title=imfile+' '+objname[0],file=imfile+'.ps', $
autosize=1,queue=queue,/negative,width=18.0,/noprint
setusym,-1
FOR i=0,countok-1 DO BEGIN
oplot,[xpos[zok[i]]],[ypos[zok[i]]],psym=8,symsize=1.5
xyouts,xpos[zok[i]],ypos[zok[i]],' '+strtrim(string(i),2),/data, $
charsize=0.5
ENDFOR
setusym,1
device,/close
display
hardcopy = queue ne ''
IF hardcopy THEN BEGIN
spawn,'lpr -P'+queue+' idl.ps'
print,'output sent to printer ',queue
ENDIF
logfile=root+'.log'
starfile=root+'.stars'
openw,lunstar,starfile,/append,/get_lun
openw,lunlist,imfile+'.lst',/get_lun
printf,lunlist,'Image: ',imfile,' object: ',objname[0]
printf,lunlist,' '
print,'Doing automatic photometry on the good stars to ',logfile
; Construct object names
starname=strarr(countok)
FOR i=0,countok-1 DO BEGIN
rastr,sra[zok[i]],1,srastr
decstr,sdec[zok[i]],0,sdecstr
starname[i]='NVt'+strmid(srastr,0,2)+strmid(srastr,3,2)+ $
strmid(srastr,6,2)+strmid(srastr,9,1)+ $
strmid(sdecstr,0,3)+strmid(sdecstr,4,2)+ $
strmid(sdecstr,7,2)
printf,lunstar,starname[i],' ',srastr,' ',sdecstr,' 0.000 0.000 '
printf,lunlist,i,' ',starname[i],' ',srastr,' ',sdecstr
ENDFOR
free_lun,lunstar
free_lun,lunlist
IF hardcopy THEN BEGIN
spawn,'lpr -P'+queue+' '+imfile+'.lst'
ENDIF
spawn,'sort '+starfile+' | uniq > tmp.xxx; mv tmp.xxx '+starfile
filter=strtrim(hdrinfo.filter,2)
if exists(logfile) then $
logmanip,logfile,/JUSTCLEAN,DELETEFILE=imfile
objnum=10
basphote,gain,image,exptime,xpos[zok],ypos[zok], $
objrad,objrad+10,objrad+30,logfile,objnum,/ALTLOG,/SILENT, $
FNAME=imfile,FILTER=filter,JD=hdrinfo.jd,NAME=starname, $
PSCALE=info.pscale
logmanip,logfile,/JUSTCLEAN
ENDIF ; photstars block
IF spot[0] gt -0.9e9 THEN BEGIN
spot=double(spot)
nspots=n_elements(spot)/2
FOR sn=0,nspots-1 DO BEGIN
xi = asteval((spot[0,sn]-info.xc0)/renormfac, $
(spot[1,sn]-info.yc0)/renormfac, $
cxi,xiterms)/3600.0d0*!dpi/180.0d0
eta= asteval((spot[0,sn]-info.xc0)/renormfac, $
(spot[1,sn]-info.yc0)/renormfac, $
ceta,etaterms)/3600.0d0*!dpi/180.0d0
astsn2rd,xi,eta,info.ra0,info.dec0,ra,dec
rastr,ra,4,ras
decstr,dec,3,decs
infostr=string(thisexp,jds,spot[0,sn],spot[1,sn], $
format='(a,1x,a,1x,f6.1,1x,f6.1)')
print,infostr,' ',ras,' ',decs
repwrite,'spot.dat',infostr,infostr+' '+ras+' '+decs
ENDFOR
ENDIF
; Look for an object list file, these are object positions that
; have been located elsewhere. Re-measure the position, generate
; ra,dec and save.
objid=0
objtag=''
IF exists(fnobj) THEN BEGIN
rdoblist,fnobj,nobj,files,offvals,xyvals,flags,nfiles
zf=where(imfile eq files)
zf=zf[0]
IF zf ne -1 THEN BEGIN
words=str_sep(files[0],'.')
objtag=strmid(words[0],strlen(root)-2,2) + $
strmid(words[1],strlen(suffix)-3,3) + extstr
for objid=0,nobj-1 do begin
tobjname=objtag+strb36(objid,pad=2)
resfile = tobjname + '.ast'
IF flags[objid] eq 'y' THEN BEGIN
x = xyvals[zf*2,objid]
y = xyvals[zf*2+1,objid]
tag=strcompress(thisexp+' '+tobjname+' '+ $
string(objrad,format='(f10.1)'))
if x ge 0.0 and y ge 0.0 then begin
basphote,gain,image,exptime,x,y,objrad,objrad+10,objrad+30, $
/exact,/nolog,/silent,xcen=xm,ycen=ym,mag=mag0,fwhm=fwhm0
oplot,[xm],[ym],psym=5,color=4,symsize=1.5
xyouts,xm,ym,' '+strtrim(string(objid),2),color=4
astmag=(mag0+magdiff) < 99.9
xi = asteval((xm-info.xc0)/renormfac,(ym-info.yc0)/renormfac, $
cxi,xiterms)/3600.0d0*!dpi/180.0d0
eta= asteval((xm-info.xc0)/renormfac,(ym-info.yc0)/renormfac, $
ceta,etaterms)/3600.0d0*!dpi/180.0d0
astsn2rd,xi,eta,info.ra0,info.dec0,ra,dec
rastr,ra,4,ras
decstr,dec,3,decs
infostr=string(thisexp,hdrinfo.jd,ras,decs,astmag, $
format='(a,1x,f13.5,1x,a,1x,a,1x,f4.1)')
print,tobjname,infostr,' FWHM ',fwhm0*info.pscale,'"', $
format='(2x,a,1x,a,a,f4.1,a)'
repwrite,resfile,thisexp,infostr
resfile='[[DEFAULT]]'
infostr=strcompress(string(xm,ym,format='(f10.3,1x,f10.3)'))
repwrite,'position.dat',tag,tag+infostr
IF photstars THEN BEGIN
basphote,gain,image,exptime,xm,ym, $
objrad,objrad+10,objrad+30, $
logfile,objnum,/ALTLOG,/SILENT,FNAME=imfile, $
FILTER=filter,JD=hdrinfo.jd,NAME=tobjname, $
PSCALE=info.pscale,/EXACT
ENDIF
endif else begin
repwrite,resfile,thisexp,''
repwrite,'position.dat',tag,''
endelse
ENDIF ELSE IF flags[objid] eq '?' THEN BEGIN
print,'WARNING: object line ',objid,' has not yet been validated'
print,' in file ',fnobj,', skipping.'
if exists(resfile) then begin
print,'rm ',resfile
spawn,'rm '+resfile
endif
ENDIF else begin
if exists(resfile) then begin
print,'rm ',resfile
spawn,'rm '+resfile
endif
endelse
endfor
ENDIF ELSE BEGIN
print,'ERROR! There is no valid file match in object file ',fnobj
ENDELSE
ENDIF ELSE BEGIN
objres=findfile(objname[0]+'*.obj')
IF objres[0] ne '' THEN BEGIN
print,'WARNING!: There was no object list file found for this image.'
print,'Object name is ',objname[0]
ENDIF
ENDELSE
objnum=0
objcnt=0
if objtag eq '' then begin
cobjname=objname[0]
endif else begin
cobjname=objtag+strb36(objid,pad=2)
endelse
IF do_objects THEN BEGIN
print,bel
theta=findgen(361.0)/!radeg
xcirc=objrad*cos(theta)
ycirc=objrad*sin(theta)
; Get ready for astrometry on unknown
doneone=0
REPEAT BEGIN
setwin,window
if fullrefresh then tv,bim
plot,[0],[1],/nodata,xmargin=[0,0],ymargin=[0,0], $
xr=[0,sz[1]-1],yr=[0,sz[2]-1],xstyle=5,ystyle=5,/noerase
if fullrefresh then begin
oplot,[border[0],sz[1]-border[1],sz[1]-border[1], $
border[0],border[0]], $
[border[3],border[3],sz[2]-border[2], $
sz[2]-border[2],border[3]],color=4
IF not noremind THEN $
oplot,[-20,20,20,-20,-20]+xloc,[-20,-20,20,20,-20]+yloc,color=4
IF doplast THEN BEGIN
oplot,astx,asty,psym=8,color=1,symsize=0.2
FOR i=0,numobj-1 DO BEGIN
xyouts,astx[i],asty[i],' '+strtrim(string(i),2),color=1
IF asterr[i]*!radeg*3600.0/info.pscale gt 5.0 THEN BEGIN
oplot,astxvar[i,*],astyvar[i,*],color=1
oplot,[astx[i],astxp1[i]],[asty[i],astyp1[i]],color=3
ENDIF ELSE BEGIN
oplot,[astx[i]],[asty[i]],psym=4,color=1,symsize=2.0
ENDELSE
ENDFOR
ENDIF
setusym,-1
FOR i=0,nstars-1 DO BEGIN
IF bad[i] eq 0 THEN BEGIN
oplot,[xpos[i]],[ypos[i]],psym=8,color=3,symsize=2.5
ENDIF ELSE BEGIN
oplot,[xpos[i]],[ypos[i]],psym=8,color=2,symsize=2.5
ENDELSE
ENDFOR
endif
if doneone then oplot,[xm],[ym],psym=4,color=4
setusym,1
print,'Left click '+cobjname+ $
' to measure (middle-new object, right-exit)'
if roam then begin
cr = string("15b) ;"
fmt="($,a,a,2x,a)"
repeat begin
cursor,x,y,2
if !mouse.button eq 0 then begin
xi = asteval((x-info.xc0)/renormfac,(y-info.yc0)/renormfac, $
cxi,xiterms)/3600.0d0*!dpi/180.0d0
eta= asteval((x-info.xc0)/renormfac,(y-info.yc0)/renormfac, $
ceta,etaterms)/3600.0d0*!dpi/180.0d0
astsn2rd,xi,eta,info.ra0,info.dec0,ra,dec
rastr,ra,4,ras
decstr,dec,3,decs
print,format=fmt,cr,ras,decs
endif
endrep until !mouse.button ne 0
print,cr,' ',cr,format=fmt
endif else begin
cursor,x,y,3
endelse
IF !mouse.button eq 1 THEN BEGIN
basphote,gain,image,exptime,x,y,objrad,objrad+10,objrad+30, $
/nolog,/silent,xcen=xm,ycen=ym,mag=mag0,fwhm=fwhm0
astmag=(mag0+magdiff) < 99.9
xm0=fix(xm+0.5)
ym0=fix(ym+0.5)
dw=fix(ceil(objrad*3.5))
; Extract image sub-section around object, watch out for image
; edges.
x10=xm0-dw
x20=xm0+dw
y10=ym0-dw
y20=ym0+dw
x1=max([x10,0])
x2=min([x20,sz[1]-1])
y1=max([y10,0])
y2=min([y20,sz[2]-1])
swd=2*dw+1
robomean,image[x1:x2,y1:y2],3.0,0.5,meansky
subim=replicate(meansky,swd,swd)
subim[x1-x10,y1-y10] = image[x1:x2,y1:y2]
IF fwhm0 gt 1.5 THEN BEGIN
radp,subim,xm-x1,ym-y1, $
r,i,rfwhm,rcoefs,rfit,fzwid=fwhm0
ENDIF ELSE BEGIN
radp,subim,xm-(xm0-dw),ym-(ym0-dw), $
r,i,rfwhm,rcoefs,rfit,fzwid=meanfwhm
ENDELSE
idx=sort(r)
setwin,12
plot,r[idx],i[idx],psym=4
oplot,r[idx],rfit[idx],color=3
xyouts,0.6,0.8,'FWHM='+ $
strcompress(string(fwhm0,format='(f6.1)'),/remove)+ $
' pixels',/normal
xyouts,0.6,0.75,'FWHM='+ $
strcompress(string(fwhm0*info.pscale, $
format='(f6.1)'),/remove)+' arcsec',/normal
zbin=6
zsz=zbin*(2*dw+1)
setwin,13,xsize=zsz,ysize=zsz,/show
spotim=subim
robomean,spotim,3.0,0.5,meansky,dummy1,meanskysig
spotim=(spotim-min(meansky-3.0*meanskysig)+10)^0.1
tv,rebin( $
bytscl(spotim,min=min(spotim),max=max(spotim), $
top=!d.n_colors-1-savec)+savec, $
zsz,zsz,/sample)
plot,[0],[1],/nodata,xmargin=[0,0],ymargin=[0,0], $
xr=[xm0-dw-0.5,xm0+dw+0.5],yr=[ym0-dw-0.5,ym0+dw+0.5], $
xstyle=5,ystyle=5,/noerase
oplot,xcirc+xm,ycirc+ym,color=1
xi = asteval((xm-info.xc0)/renormfac,(ym-info.yc0)/renormfac, $
cxi,xiterms)/3600.0d0*!dpi/180.0d0
eta= asteval((xm-info.xc0)/renormfac,(ym-info.yc0)/renormfac, $
ceta,etaterms)/3600.0d0*!dpi/180.0d0
astsn2rd,xi,eta,info.ra0,info.dec0,ra,dec
rastr,ra,4,ras
decstr,dec,3,decs
infostr=string(thisexp,hdrinfo.jd,ras,decs,astmag, $
format='(a,1x,f13.5,1x,a,1x,a,1x,f4.1)')
print,infostr,' FWHM ',fwhm0*info.pscale,' arcsec', $
format='(3x,a,a,f4.1,a)'
doneone=1
ENDIF ELSE BEGIN ; new object or quit
; Save data if something was measured.
IF doneone THEN BEGIN
; Save final astrometry to the object file.
IF resfile eq '[[DEFAULT]]' THEN BEGIN
resfile = cobjname + '.ast'
ENDIF
repwrite,resfile,thisexp,infostr
resfile='[[DEFAULT]]'
; Save the x,y for the object to a single master file.
infostr=strcompress(string(xm,ym,format='(f10.3,1x,f10.3)'))
tag=strcompress(thisexp+' '+cobjname+' '+ $
string(objrad,format='(f10.1)'))
repwrite,'position.dat',tag,tag+infostr
doneone=0
IF photstars THEN BEGIN
basphote,gain,image,exptime,xm,ym, $
objrad,objrad+10,objrad+30, $
logfile,objnum,/ALTLOG,/SILENT,FNAME=imfile, $
FILTER=filter,JD=hdrinfo.jd,NAME=cobjname, $
PSCALE=info.pscale
ENDIF
ENDIF
; Ask for new name if middle click.
IF !mouse.button eq 2 THEN BEGIN
objcnt=objcnt+1
IF objcnt ge n_elements(objname) and objtag eq '' THEN BEGIN
read,cobjname,prompt='New object name? '
ENDIF ELSE IF objtag ne '' THEN BEGIN
objid=objid+1
cobjname=objtag+strb36(objid,pad=2)
ENDIF ELSE BEGIN
print,bel
cobjname=objname[objcnt]
ENDELSE
cobjname=nobname(strtrim(strcompress(cobjname),2))
ENDIF
ENDELSE
ENDREP UNTIL !mouse.button eq 4 ; loop on measuring object
ENDIF ; do_objects block
ENDFOR ; sub-exposure loop
;==================================
; Generate a final pretty postcript plot output file if requested
IF pretty ne '[[DEFAULT]]' THEN BEGIN
print,' --> creating final postscript plot file....'
hardim,bytscl(image,min=lowval,max=hival,top=255),0,255, $
title=imfile+' '+strupcase(cobjname),/color, $
autosize=1,queue='chani',width=18.0,/noprint,file=pretty
tvlct,red,gre,blu
oplot,[border[0],sz[1]-border[1],sz[1]-border[1],border[0],border[0]], $
[border[3],border[3],sz[2]-border[2],sz[2]-border[2],border[3]], $
color=4
IF not noremind THEN $
oplot,[-20,20,20,-20,-20]+xloc,[-20,-20,20,20,-20]+yloc,color=4
IF doplast THEN BEGIN
oplot,astx,asty,psym=8,color=1,symsize=0.2
FOR i=0,numobj-1 DO BEGIN
xyouts,astx[i],asty[i],' '+strtrim(string(i),2),color=1
IF asterr[i]*!radeg*3600.0/info.pscale gt 5.0 THEN BEGIN
oplot,astxvar[i,*],astyvar[i,*],color=1
oplot,[astx[i],astxp1[i]],[asty[i],astyp1[i]],color=3
ENDIF ELSE BEGIN
oplot,[astx[i]],[asty[i]],psym=4,color=1,symsize=2.0
ENDELSE
ENDFOR
ENDIF
setusym,-1
FOR i=0,nstars-1 DO BEGIN
IF bad[i] eq 0 THEN BEGIN
oplot,[xpos[i]],[ypos[i]],psym=8,color=1,symsize=1.5
ENDIF ELSE BEGIN
oplot,[xpos[i]],[ypos[i]],psym=8,color=2,symsize=1.5
ENDELSE
ENDFOR
device,/close
print,' --> postscript device now closed.'
display
ENDIF
ENDFOR ; End of extension for loop block
end