pro xray_contour, im, hdr, im2, hdr2, TYPE=type, PUTINFO=putinfo, $ XTITLE=xtitle, $ YTITLE=ytitle, NLEVELS = nlevels, MAX_VALUE=max_value, LEVELS=levels, $ YMINOR = yminor, SUBTITLE = subtitle, FOLLOW= follow, TITLE = title, $ XMINOR = xminor, CHAN = chan, TICKLEN = TICKLEN, L2 = l2, $ NL2 = nl2, NOPROJ = NOPROJ, GRID = GRID, VERBOSE = VERBOSE, OVER = OVER,$ COLOR1=COLOR1, color2=color2, asca=asca, equinox=equinox,stack=stack, $ ra=ra,dec=dec,cl=cl,$ raoff1 = raoff1,decoff1 = decoff1, raoff2 = raoff2, decoff2 = decoff2,$ C_COLOR1=C_COLOR1,C_COLOR2=C_COLOR2,fudge=fudge,nogrid=nogrid, $ noaxis=noaxis ;+ ; NAME: ; XRAY_CONTOUR2 ; PURPOSE: ; Contour plot labeled with astronomical coordinates. The type ; of coordinate display is controlled by the keyword TYPE ; Set TYPE=0 (default) to measure distances from the center of the image ; (XRAY_CONTOUR will decide whether the plotting units will be in ; arc seconds, arc minutes, or degrees depending on image size.) ; Set TYPE=1 for standard RA and Dec labeling ; ; CALLING SEQUENCE: ; XRAY_CONTOUR, im, hdr,[ TYPE=, PUTINFO=, XTITLE = , YTITLE = , NLEVELS= ; MAX_VALUE=, LEVELS=, TITLE =, SUBTITLE =, FOLLOW = , XMINOR =, ; YMINOR =, CHAN = ] ; ; 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 ; The other keywords XTITLE, YTITLE, NLEVELS, MAX_VALUE, LEVELS, ; TITLE, SUBTITLE, FOLLOW, XMINOR, and YMINOR have the same ; meaning as in the CONTOUR procedure. ; CHAN - On devices with multiple channels this specifies the output channel ; On the IVAS and DeAnza the default is the graphics overlay. ; ; 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 ; 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 ;- On_error,2 ;Return to caller if N_params() LT 2 then begin ;Sufficient parameters? print,'Syntax - imcontour, im, hdr, [ TYPE= ,PUTINFO= ]' return endif twoim = 1 if n_params() lt 4 then twoim = 0 fix=0.0 if keyword_set(fudge) then begin print,'Adding 0.5 to CRPIX values to compensate for difference between' print,'IDL and FITS convensions for coordinates of the center of a pixel' print,'See Legacy, No. 6, Aug 1995 pg 38 for details' fix = 0.5 endif check_fits, im, hdr, dimen, /NOTYPE ;Make sure header appropiate to image ; if !ERR EQ -1 then return ; Set defaults if keywords not set if not keyword_set(asca) then asca = 0 im_max = max( im, MIN=im_min ) ;Image MAX and MIN values if not keyword_set( LEVELS ) then $ ;Default is 6 equally spaced levels levels = im_min + (findgen(6)+1)*(float(im_max)-im_min)/7. if not keyword_set(NLEVELS) then $ ;Default is NLEVELS = 6 nlevels = N_elements(levels) noeras = keyword_set(over) verbose = keyword_set(verbose) if keyword_set(grid) then begin noeras = 1 nodata = 1 if n_elements(cl) eq 0 then cl = 0 endif else nodata = 0 if n_elements(cl) eq 0 then cl = 3 if twoim then begin check_fits, im2, hdr2, dimen2, /NOTYPE ;Make sure header appropiate to image if !ERR EQ -1 then return im_max2 = max( im2, MIN=im_min2 ) ;Image MAX and MIN values if not keyword_set( L2 ) then $ ;Default is 6 equally spaced levels l2 = im_min2 + (findgen(6)+1)*(float(im_max2)-im_min2)/7. if not keyword_set(NL2) then $ ;Default is NLEVELS = 6 nl2 = N_elements(l2) endif prj = not keyword_set(noproj) if prj then print,'Using tangental projection' if not keyword_set(MAX_VALUE) then max_value = im_max if not keyword_set( FOLLOW ) then follow = 1 else if follow eq 2 then follow=0 if not keyword_set( TYPE ) then type = 1 if not keyword_set( TICKLEN ) then TICKLEN= 0.02 if not keyword_set(PUTINFO) then putinfo = 0 if not keyword_set(XMINOR) then $ if !X.MINOR EQ 0 then xminor = 5 else xminor = !X.MINOR if not keyword_set(YMINOR) then $ if !Y.MINOR EQ 0 then yminor = 5 else yminor = !Y.MINOR if N_elements(chan) NE 1 then begin ;Device have a graphics overlay? case !D.NAME of 'DEANZA': chan = 4 'IVAS': chan = 3 ELSE: chan = 0 endcase endif if n_elements(raoff1) eq 0 then raoff1 = 0.0 if n_elements(decoff1) eq 0 then decoff1 = 0.0 if n_elements(raoff2) eq 0 then raoff2 = 0.0 if n_elements(decoff2) eq 0 then decoff2 = 0.0 extast, hdr, astr, noparams astr.crpix = astr.crpix + fix ; extast, hdr, cd, crpix, crval, noparams ;Extract astrometry from header ; if n_elements(cdsc) eq 0 then cdsc = 1.0 ; cd = cd*cdsc equ = sxpar(hdr,'epoch') if equ eq 0 then equ = sxpar(hdr,'equinox') if equ eq 0 then equ = 2000.0 print,'Equinox of image 1 = ',strn(equ) if keyword_set(equinox) then if equinox ne equ then begin oldra = astr.crval(0) & olddec = astr.crval(1) precess,oldra,olddec,equ,equinox astr.crval(0) = oldra astr.crval(1) = olddec print,"Precessed image 1's coords to ",strn(equinox) equ = equinox endif print,'raoff1, decoff1 = ',raoff1,decoff1 astr.crval(0) = astr.crval(0) - raoff1 astr.crval(1) = astr.crval(1) - decoff1 if asca eq 1 then astr.crpix = astr.crpix-1.0 if noparams LT 0 then $ ;Does astrometry exist? message,'FITS header does not contain astrometry' ; if prj then begin ; astr.cdelt = astr.cdelt/!RADEG & astr.crval = astr.crval/!RADEG ; endif if twoim then begin extast, hdr2, astr2, noparams2 astr2.crpix = astr2.crpix + fix ; extast, hdr2, cd2, crpix2, crval2, noparams2 ;Extract astrometry from header equ2 = sxpar(hdr2,'epoch') if equ2 eq 0 then equ2 = sxpar(hdr2,'equinox') if equ2 eq 0 then equ2 = 2000.0 print,'Equinox of image 2 = ',strn(equ2) if equ2 ne equ then begin oldra = astr2.crval(0) & olddec = astr2.crval(1) precess,oldra,olddec,equ2,equ astr2.crval(0) = oldra astr2.crval(1) = olddec print,"Precessed image 2's coords to image 1's equinox" endif astr2.crval(0) = astr2.crval(0) - raoff2 astr2.crval(1) = astr2.crval(1) - decoff2 if asca eq 2 then astr2.crpix = astr2.crpix-1.0 if noparams2 LT 0 then $ ;Does astrometry exist? message,'FITS header does not contain astrometry' ; if prj then begin ; astr2.cd = astr2.cd/!RADEG & astr2.crval = astr2.crval/!RADEG ; endif endif if verbose then print,'astr.crval(0) = ',strn(astr.crval(0)),$ ', astr.crval(1) = ',strn(astr.crval(1)) ; Adjust plotting window so that contour plot will have same dimensional ; ratio as the image xlength = !D.X_VSIZE & ylength = !D.Y_VSIZE xsize = fix( dimen(0) ) & ysize = fix( dimen(1) ) if twoim then begin xsize2 = fix( dimen2(0) ) & ysize2 = fix( dimen2(1) ) endif xsize1 = xsize-1 & ysize1 = ysize-1 xratio = xsize / float(ysize) yratio = ysize / float(xsize) if ( ylength*xratio LT xlength ) then begin xmax = 0.1+0.8*ylength*xratio/xlength ; xmax = 0.15 + 0.8*ylength*xratio/xlength ; pos = [ 0.15, 0.15, xmax, 0.95 ] pos = [0.1,0.1,xmax,0.9] if not putinfo then pos=[0.1,0.1,0.9,0.9] endif else begin xmax = 0.95 print,noeras if noeras then pos = [ 0.1, 0.1, 0.9, 0.1+0.8*xlength*yratio/ylength ] $ else begin if putinfo then pos = [0.15,0.15,xmax,0.15+0.8*xlength*yratio/ylength] $ else pos=[0.1,0.1,0.9,0.9] endelse ; pos = [0.0,0.0,0.8,0.8] print,pos ; stop endelse if !X.TICKS GT 0 then xtics = abs(!X.TICKS) else xtics = 8 if !Y.TICKS GT 0 then ytics = abs(!Y.TICKS) else ytics = 8 pixx = xsize/xtics ;Number of X pixels between tic marks pixy = ysize/ytics ;Number of Y pixels between tic marks getrot,astr,rot,cdelt ;Get the rotation and plate scale ; degx = pixx/cd(0,0) ;No. of degrees between tic marks ; degy = pixy/cd(1,1) xmid = xsize1/2. & ymid = ysize1/2. if prj then $ xy2ad,xmid,ymid,astr,ra_cen,dec_cen $ ;Get Ra and Dec of image center else begin dec_cen = (ymid-astr.crpix(1)+1)*astr.cdelt(1)+astr.crval(1) ra_cen = (xmid-astr.crpix(0)+1)*astr.cdelt(0)/cos(dec_cen/!radeg)+$ astr.crval(0) ; dec_cen = sxpar(hdr,'crval2') endelse ra_dec = adstring(ra_cen,dec_cen,1) ;Make a nice string ; Determine tic positions and labels for the different type of contour plots if type NE 0 then begin ;RA and Dec labeling xedge = [0,xsize1,0] ;X pixel values of the four corners yedge = [0,0,ysize1] ;Y pixel values of the four corners if prj then begin ; xy2ad, xedge, yedge, cd, crpix, crval, a, d ;RA and Dec of the corners xy2ad, xedge, yedge, astr, a, d ; a = a*!RADEG & d = d*!RADEG xy2ad,findgen(xsize),replicate(0,xsize),astr,ra,dummy ; xy2ad,findgen(xsize),replicate(0,xsize),cd,crpix,crval,ra,dummy xy2ad,replicate(0,ysize),findgen(ysize),astr,dummy,dec ; xy2ad,replicate(0,ysize),findgen(ysize),cd,crpix,crval,dummy,dec ; ra = ra*!radeg & dec = dec*!radeg endif else begin d = (yedge-astr.crpix(1)+1)*astr.cdelt(1)+astr.crval(1) a = (xedge-astr.crpix(0)+1)*astr.cdelt(0)/cos(d/!radeg)+astr.crval(0) dec = (findgen(ysize)-astr.crpix(1)+1)*astr.cdelt(1)+astr.crval(1) ra = (findgen(xsize)-astr.crpix(0)+1)*astr.cdelt(0)/cos(dec(0)/!radeg)+$ astr.crval(0) endelse if twoim then begin if prj then begin ; xy2ad,findgen(xsize2),findgen(ysize2),cd2,crpix2,crval2,ra2,dec2 ; xy2ad,findgen(xsize2),findgen(ysize2),astr2,ra2,dec2 xy2ad,findgen(xsize2),replicate(0,xsize2),astr2,ra2,dummy xy2ad,replicate(0,ysize2),findgen(ysize2),astr2,dummy,dec2 ; ra2 = ra2*!radeg & dec2 = dec2*!radeg endif else begin dec2 = (findgen(ysize2)-astr2.crpix(1)+1)*astr2.cdelt(1)+astr2.crval(1) ra2 = (findgen(xsize2)-astr2.crpix(0)+1)*astr2.cdelt(0)/$ cos(dec2/!radeg)+astr2.crval(0) endelse endif pixx = xsize/xtics ;Number of X pixels between tic marks pixy = ysize/ytics ;Number of Y pixels between tic marks tics, a(0), a(1), xsize, pixx, raincr,/RA ;Find an even increment for RA tics, d(0), d(2), ysize, pixy, decincr ;Find an even increment for Dec tic_one, a(0), pixx, raincr, botmin, xtic1, /RA ;Position of first RA tic tic_one, d(0), pixy, decincr,leftmin,ytic1 ;Position of first Dec tic nx = fix( (xsize1-xtic1-1)/pixx ) ;Number of X tic marks ny = fix( (ysize1-ytic1-1)/pixy ) ;Number of Y tic marks ra_grid = (botmin + findgen(nx+1)*raincr/4.);/!RADEG dec_grid = (leftmin + findgen(ny+1)*decincr/60.);/!RADEG ticlabels, botmin, nx+1, raincr, xlab, /RA, DELTA=1 ticlabels, leftmin, ny+1, decincr, ylab,DELTA=1 ;stop if prj then begin ; xpos = cons_ra(ra_grid/!radeg,0,cd,crpix,crval) ;Line of constant RA ; xpos = cons_ra(ra_grid/!radeg,0,astr) xpos = cons_ra(ra_grid,0,astr) ; ypos = cons_dec(dec_grid/!radeg,0,cd,crpix,crval) ;Line of constant Dec ; ypos = cons_dec(dec_grid/!radeg,0,astr) ypos = cons_dec(dec_grid,0,astr) endif else begin xpos = cos(dec_grid/!radeg)*(ra_grid-astr.crval(0))/astr.cdelt(0) + $ astr.crpix(0) ypos = (dec_grid-astr.crval(1))/astr.cdelt(1) + astr.crpix(1) endelse xunits = 'RIGHT ASCENSION' yunits = 'DECLINATION' endif else begin ;Label with distance from center ticpos, xsize1*astr.cdelt(0), xsize, pixx, incrx, xunits numx = fix(xsize/(2.*pixx)) ticpos, ysize1*astr.cdelt(0), ysize, pixy, incry, yunits numy = fix(ysize/(2.*pixy)) nx = 2*numx & ny = 2*numy xpos = xmid + (findgen(nx+1)-numx)*pixx ypos = ymid + (findgen(ny+1)-numy)*pixy xlab = string(indgen(nx+1)*incrx - incrx*numx,'(I3)') ylab = string(indgen(ny+1)*incry - incry*numy,'(I3)') endelse ; Get default values of XTITLE, YTITLE, TITLE and SUBTITLE if N_elements(xtitle) EQ 0 then $ if !X.TITLE eq '' then xtitle = xunits else xtitle = !X.TITLE if N_elements(ytitle) EQ 0 then $ if !Y.TITLE eq '' then ytitle = yunits else ytitle = !Y.TITLE if N_elements(title) EQ 0 then title = !P.TITLE if n_elements(color1) eq 0 then color1 = !color if n_elements(color2) eq 0 then color2 = !color if (not keyword_set( SUBTITLE) ) and (putinfo LT 1) then $ subtitle = 'CENTER: R.A. '+ strmid(ra_dec,1,13)+' DEC ' + $ strmid(ra_dec,13,13) if (not keyword_set( SUBTITLE) ) then subtitle = !P.SUBTITLE clab = levels clab(*) = 0 xtitle = xtitle + ' ('+string(equ,format='(f6.1)')+')' ytitle = ytitle + ' ('+string(equ,format='(f6.1)')+')' xcharsize=1 ycharsize=1 if n_elements(stack) gt 0 then begin xcharsize=0.5 ycharsize=0.5 if stack eq 1 then pos=[0.1,0.1,0.5,0.5] else begin xtitle='' xlab(*) = ' ' pos=[0.1,0.5,0.5,0.9] noeras=1 endelse endif print,'pos = ',pos if n_elements(c_color1) eq 0 then begin c_color1 = levels c_color1(*) = color1 endif xst = 1 & yst = 1 if keyword_set(noaxis) then begin xst = 4 & yst = 4 endif contour, im, ra, dec,$ XTICKS = nx, YTICKS = ny, POSITION=pos, $ XTICKV = ra_grid, YTICKV = dec_grid, $ XTITLE=xtitle, YTITLE=ytitle, COLOR = color1, $ ; XTICKV = xpos, YTICKV = ypos,XTITLE=xtitle, YTITLE=ytitle, $ XTICKNAME = xlab, YTICKNAME = ylab, TITLE = title, $ LEVELS = levels, NLEVELS = nlevels, MAX_VALUE=max_value, $ SUBTITLE = subtitle, FOLLOW = follow, XMINOR = xminor, $ YMINOR = yminor, CHAN = chan, TICKLEN = ticklen,c_labels=clab,$ XRANGE = [ra(0),ra(xsize1)],noeras=noeras,nodata=nodata,$ XCHARSIZE=xcharsize,YCHARSIZE=ycharsize,C_COLOR=c_color1, $ xstyle = xst, ystyle = yst ; if not prj then begin ; astr.cdelt = astr.cdelt/!radeg & astr.crval = astr.crval/!radeg ; endif xind = findgen(xsize) if not keyword_set(nogrid) then for i=0,ny do begin ; yind = cons_dec(dec_grid(i)/!radeg,xind,cd,crpix,crval) ; yind = cons_dec(dec_grid(i)/!radeg,xind,astr) yind = cons_dec(dec_grid(i),xind,astr) ; dec_yind = yind*astr.cdelt(1)*!radeg + dec(0) dec_yind = yind*astr.cdelt(1) + dec(0) ; ra_xind = xind*astr.cdelt(0)*!radeg/cos(dec(0)/!radeg) + ra(0) ra_xind = xind*astr.cdelt(0)/cos(dec(0)/!radeg) + ra(0) oplot,ra_xind,dec_yind ; oplot,ra(xind+0.5),dec(yind+0.5) ; print,'Grid ra,dec: ',adstring(ra(0),dec_grid(i)) ; print,'First point in grid line: ',adstring(ra_xind(0),dec_yind(0)) if verbose then print,'Range of grid line in DEC in arcsecs: ',$ ; strn(abs(max(yind)-min(yind))*astr.cdelt(1)*!radeg*3600) strn(abs(max(yind)-min(yind))*astr.cdelt(1)*3600) endfor yind = findgen(ysize) if not keyword_set(nogrid) then for i=0,nx do begin ; xind = cons_ra(ra_grid(i)/!radeg,yind,cd,crpix,crval) ; xind = cons_ra(ra_grid(i)/!radeg,yind,astr) xind = cons_ra(ra_grid(i),yind,astr) ; dec_yind = yind*astr.cdelt(1)*!radeg + dec(0) dec_yind = yind*astr.cdelt(1) + dec(0) ; ra_xind = xind*astr.cdelt(0)*!radeg/cos(dec(0)/!radeg) + ra(0) ra_xind = xind*astr.cdelt(0)/cos(dec(0)/!radeg) + ra(0) oplot,ra_xind,dec_yind ; oplot,ra(xind+0.5),dec(yind+0.5) ; print,'Grid ra,dec: ',adstring(ra_grid(i),dec(0)) ; print,'First point in grid line: ',adstring(ra_xind(0),dec_yind(0)) if verbose then print,'Range of grid line in RA in arcsecs: ',$ ; strn(abs(max(xind)-min(xind))*astr.cdelt(0)/cos(dec(0)/!radeg)*!radeg*3600) strn(abs(max(xind)-min(xind))*astr.cdelt(0)/cos(dec(0)/!radeg)*3600) endfor if twoim then begin clab = l2 clab(*) = 0 if n_elements(c_color2) eq 0 then begin c_color2 = levels c_color2(*) = color2 endif contour, im2, ra2, dec2,$ /OVER, LEVELS = l2, NLEVELS = nl2, FOLLOW = follow, c_line = cl,$ color=color2,c_labels=clab,c_color=c_color2 endif ; Write info about the contour plot if desired if putinfo GE 1 then begin xmax = xmax + 0.01 object = sxpar( hdr, 'OBJECT' ) if !ERR ne -1 then xyouts, xmax, 0.95, object, /NORM, CHAN = chan name = sxpar( hdr, 'IMAGE' ) if !ERR ne -1 then xyouts,xmax,0.90,name, /NORM, CHAN = chan xyouts, xmax, 0.85,'CENTER:',/NORM, CHAN = chan xyouts, xmax, 0.80, 'R.A. '+ strmid(ra_dec,1,13),/NORM, CHAN = chan xyouts, xmax, 0.75, 'DEC '+ strmid(ra_dec,13,13),/NORM, CHAN = chan xyouts, xmax, 0.70, 'IMAGE SIZE', /NORM, CHAN = chan xyouts, xmax, 0.65, 'X: ' + strtrim(xsize,2), /NORM, CHAN = chan xyouts, xmax, 0.60, 'Y: ' + strtrim(ysize,2), /NORM, CHAN = chan xyouts, xmax, 0.50, strmid(!STIME,0,17),/NORM, CHAN = chan xyouts, xmax, 0.40, 'CONTOUR LEVELS:',/NORM, CHAN = chan for i = 0,(nlevels < 6)-1 do $ xyouts,xmax,0.35-0.05*i,string(i,'(i2)') + ':' + $ string(levels(i)), /NORM, CHAN = chan endif return end