program xcdisp6 c c modification of xcdisp5, adding handling of different image sizes c c bug: some fields (eg 1447a0A) require different axffiles for direct c and spectral images. could enter two, w/default being the same. c include "cdlftn.inc" character*40 mname, axffile character*40 dname, maskname, directim character*40 mnamesave, dnamesave character key*1 character str*10, str4*4, str2*2 character*2 ch2arr(5) data ch2arr /'U', 'B', 'g', 'R', 'I'/ integer map, i, j, iap, ii, jj parameter (map = 200) real*4 cx, cy, dr, dl, rm, dy, col real*4 cenx(map), ceny(map), rmag(map) real*4 dright(map), dleft(map), dety(map) real*4 color(map) integer ippp(map), ireals(map), ip, irs integer iframe, ifb integer ifill, icolor integer m, ier, ix1, ix2, iy1, iy2 c parameter (m = 2048) parameter (m = 1024) integer mx, my, iyspec, ixim, iyim integer*2 ibspec, ibdir, ibdir4(4) real z1, z2, x, y, txsize, txsize2 integer nxpad, nypad, nxbox, nybox, nspec integer ispecshift, ixoff, iyoff integer xspecshift, yspecshift logical firstspec, done integer nlines, isize, ixOI, ixNa real*4 w55, wNa, disp, w, w0, z parameter (w55 = 5577., wNa = 5892., nlines = 7) real*4 wline(nlines) character*1 type(nlines) character*4 label(nlines) data wline /3727.3, 3933.7, 3968.5, 4861.33, $ 4958.91, 5006.84, 6562.82/ data type /'e', 'a', 'a', 'e', 'e', 'e', 'e'/ data label /'3727', 'Ca K', 'Ca H', ' Hb ', $ '4959', '5007', ' Ha '/ c c initialize cdl package and saotng display c call cfopen (0, ier) cc call cfopen ("unix:/tmp/.IMT%d", ier) cc call cfopen ("fifo:/dev/imt1i:/dev/imt1o", ier) if (ier .gt. 0) then write (*,*) 'open: Error return from CDL' goto 999 endif iframe = 1 cc ifb = 5 ifb = 3 call cfsetframe(iframe) call cfsetfbconfig(ifb) call cfgetfbconfig(ifb, nw, nh, nf) c c read input parameters c open(11, file='xcdisp6.in', status='old') read(11,'(a40)') maskname write(*,*) write(*,'(a7,a40)') 'mask = ', maskname mnamesave = maskname read(11,'(a40)') axffile read(11,'(a40)') directim write(*,'(a15,a40)') 'direct image = ', directim dnamesave = directim read(11,*) z1, z2 read(11,*) izt write(*,*) write(*,*) 'scale limits = ', z1, z2 if (izt .lt. 0 .or. izt .gt. 2) then write(*,*) 'warning: unknown izt = ', izt izt = 2 write(*,*) 'resetting to log10: izt = ', izt else write(*,*) 'izt = ', izt end if read(11,*) ibspec, ibdir write(*,*) 'subtract ', ibspec, ' from spectral images' write(*,*) 'subtract ', ibdir, ' from direct images' read(11,*) (ibdir4(j),j=1,4) write(*,*) 'subtract ', (ibdir4(j),j=1,4), $ ' from UBgI direct images' read(11,*) ispecshift read(11,*) nxpad, nypad read(11,*) nxbox, nybox write(*,*) write(*,*) 'ispecshift = ', ispecshift write(*,*) 'nxpad, nypad = ', nxpad, nypad write(*,*) 'nxbox, nybox = ', nxbox, nybox read(11,*) isize c c new lines c read(11,*) imxsize, imysize read(11,*) xspecshift, yspecshift read(11,*) nspec write(*,*) write(*,*) 'imxsize, imysize = ', imxsize, imysize write(*,*) 'xspecshift, yspecshift = ', xspecshift, yspecshift write(*,*) 'nspec = ', nspec close(11) c nspec = 3 c c read axffile c open(15, file=axffile, status='old') write(*,*) write(*,'(a16,a40)') 'reading axffile ', axffile nap = 0 read(15,*,iostat=ios) i, cx, cy, dr, dl, $ ip, rm, irs, dy, col do while (ios .eq. 0) nap = nap + 1 if (nap .gt. map) then write(*,*) 'error: too many apertures' write(*,*) 'nap > map: ', nap, map stop end if if (i .ne. nap) then write(*,*) 'error: aperture mismatch' write(*,*) 'i != nap: ', i, nap stop end if cenx(i) = cx ceny(i) = cy dright(i) = dr dleft(i) = dl ippp(i) = ip rmag(i) = rm ireals(i) = irs dety(i) = dy color(i) = col read(15,*,iostat=ios) i, cx, cy, dr, dl, $ ip, rm, irs, dy, col end do write(*,*) 'read ', nap, ' apertures' close(15) txsize = 1. iap = 0 firstspec = .true. ixOI = 0 ixNa = 0 disp = 0. c c wait for / act on image cursor input c done = .false. do while (.not. done) write(*,*) write(*,*) 'enter command' c write(*,*) 'ixOI = ', ixOI call cfreadcursor (0, x, y, key, ier) if (ier .gt. 0) then write (*,*) 'cfreadCursor: Error return from CDL' else c ------------------------------------------------ c choose aperture & display (dispap) c if (key .eq. 'a' .or. $ key .eq. 'n' .or. $ key .eq. 'p' .or. $ key .eq. 'r') then if (key .eq. 'a') then iap = 0 do while (iap .lt. 1 .or. iap .gt. nap) write(*,*) 'enter aperture: ' read(*,*) iap if (iap .lt. 1 .or. iap .gt. nap) then write(*,*) 'out of bounds; try again' end if end do else if (key .eq. 'n') then if (iap .lt. nap) then iap = iap + 1 else iap = nap write(*,*) 'defaulting to last aperture ', iap end if write(*,*) write(*,*) 'next aperture ', iap else if (key .eq. 'p') then if (iap .gt. 1) then iap = iap - 1 else iap = 1 write(*,*) 'defaulting to first aperture ', iap end if write(*,*) write(*,*) 'previous aperture ', iap else if (key .eq. 'r') then write(*,*) if (iap .lt. 1 .or. iap .gt. nap) then iap = 1 write(*,*) 'defaulting to aperture ', iap end if write(*,*) 'redrawing aperture ', iap end if c write(*,*) 'ixOI = ', ixOI write(*,'(''iap '', i3, '' ppp # '', i5, $ '' mag '', f6.2, '' obj '', i5, $ '' color '', f6.2)') $ iap, ippp(iap), rmag(iap), $ ireals(iap), color(iap) mname = maskname dname = directim c write(*,*) 'mnamesave = ', mnamesave c write(*,*) 'dnamesave = ', dnamesave c write(*,*) 'directim = ', directim c write(*,*) 'mname = ', mname c write(*,*) 'dname = ', dname call dispap(mname, dname, iap, $ cenx(iap), ceny(iap), dright(iap), $ dleft(iap), dety(iap), z1, z2, izt, $ ibspec, ibdir, ibdir4, $ nxpad, nypad, nxbox, nybox, $ nspec, firstspec, mx, my, iyspec, ispecshift, $ ixim, iyim, ifb, xspecshift, yspecshift) c write(*,*) 'ixOI = ', ixOI if (key .ne. 'r') then z = -999. end if ixoff = nw/2 - mx/2 iyoff = nh/2 - my/2 ifill = 0 icolor = c_red ix1 = mx/2 - 3 + ixoff ix2 = mx/2 + 3 + ixoff iy1 = iyspec - 3 + iyoff iy2 = iyspec + 3 + iyoff c write(*,*) 'marking box at: ', ix1, iy1, ix2, iy2 c call cfmarkbox(ix1, iy1, ix2, iy2, c $ ifill, icolor, ier) jj = iy1 c call cfmarkline(ix1, iy1, ix2, jj, icolor, ier) call cfmarkbox(ix1, iy1, ix2, jj, ifill, icolor, ier) jj = iy2 c call cfmarkline(ix1, iy2, ix2, jj, icolor, ier) call cfmarkbox(ix1, iy2, ix2, jj, ifill, icolor, ier) ix1 = mx/2 - 20 + ixoff iy1 = iyspec + 10 $ + nint(abs(dright(iap))) + iyoff write(str,'(i3,2x,i5)') iap, ireals(iap) call cfmarktext(ix1, iy1, str, txsize, $ 0., icolor, ier) icolor = c_yellow ix1 = ixim - 3 + ixoff ix2 = ixim + 3 + ixoff iy1 = iyim $ - nint(abs(dleft(iap))) + iyoff iy2 = iyim $ + nint(abs(dright(iap))) + iyoff call cfmarkbox(ix1, iy1, ix2, iy2, ifill, icolor, ier) icolor = c_red c iy1 = iyim - 3 + iyoff c iy2 = iyim + 3 + iyoff ix1 = ixim + ixoff iy1 = iyim + iyoff call cfmarkpoint(ix1, iy1, 0, 7, $ m_box, icolor, ier) c call cfmarkbox(ix1, iy1, ix2, iy2, c $ ifill, icolor, ier) txsize2 = 1.0 iy1 = iyim + iyoff do ifilt = 1, 3 ix1 = ixim + ixoff - (2*nybox+1) - 10*(3-ifilt) $ - (16-4*ifilt)*nybox/2 str2 = ch2arr(ifilt) c call cfmarktext(ix1, iy1, str2, txsize2, c $ 0., icolor, ier) end do do ifilt = 4, 5 ix1 = ixim + ixoff + 10*(ifilt-3) $ + (4*ifilt-14)*nybox/2 str2 = ch2arr(ifilt) c call cfmarktext(ix1, iy1, str2, txsize2, c $ 0., icolor, ier) end do ix1 = ixim + ixoff iy1 = iyim + iyoff c ------------------------------------------------ c print info on current aperture c else if (key .eq. 'i') then if (iap .ge. 1 .and. iap .le. nap) then write(*,'(''iap '', i3, '' ppp # '', i5, $ '' mag '', f6.2, '' obj '', i5, $ '' color '', f6.2)') $ iap, ippp(iap), rmag(iap), $ ireals(iap), color(iap) else write(*,*) 'warning: aperture ', iap, $ ' out of range' end if c else if (key .eq. 'c') then c write(*,*) c write(*,*) 'clearing all markers' c call cfclearoverlay(ier) c ------------------------------------------------ c other keys c else if (key .eq. 'd') then ix1 = nint(x) + ixoff iy1 = nint(y) + iyoff write(*,*) write(*,*) 'deleting marker at ', x, y call cfdeletemark(ix1, iy1, ier) else if (key .eq. 'q') then done = .true. write(*,*) write(*,*) 'done xcdisp' else if (key .eq. 's') then write(*,*) if (ixOI .ne. 0) then c write(*,*) 'deleting existing sky lines' ix1 = ixOI + ixoff iy1 = iyspec + iyoff c call cfdeletemark(ix1, iy1, ier) ix1 = ixNa + ixoff c call cfdeletemark(ix1, iy1, ier) end if ixOI = nint(x) write(*,*) 'marking 5577 sky at x = ', ixOI ix1 = ixOI + ixoff ix2 = ix1 iy1 = iyspec $ + nint(abs(dright(iap))) + iyoff iy2 = iyspec $ - nint(abs(dleft(iap))) + iyoff icolor = c_green call cfmarkbox(ix1, iy1, ix2, iy2, ifill, icolor, ier) write(*,*) 'now mark Na 5892' call cfreadcursor (0, x, y, key, ier) ixNa = nint(x) write(*,*) 'marking Na 5892 sky at x = ', ixNa ix1 = ixNa + ixoff ix2 = ix1 icolor = c_coral call cfmarkbox(ix1, iy1, ix2, iy2, ifill, icolor, ier) if (ixNa .le. ixOI) then write(*,*) 'warning: Na position <= 5577 position' write(*,*) 'Na, 5577: ', ixNa, ixOI disp = 0. else disp = (wNa - w55) / real(ixNa - ixOI) end if write(*,*) 'dispersion = ', disp, ' Angstroms/pixel' else if (key .eq. 't') then write(*,*) if (ixOI .eq. 0) then write(*,*) 'warning: no existing sky lines' else write(*,*) 'marking current sky lines' write(*,*) 'marking 5577 sky at x = ', ixOI ix1 = ixOI + ixoff ix2 = ix1 iy1 = iyspec $ + nint(abs(dright(iap))) + iyoff iy2 = iyspec $ - nint(abs(dleft(iap))) + iyoff icolor = c_green call cfmarkbox(ix1, iy1, ix2, iy2, $ ifill, icolor, ier) ix1 = ixNa + ixoff ix2 = ix1 write(*,*) 'marking Na 5892 sky at x = ', ixNa icolor = c_coral call cfmarkbox(ix1, iy1, ix2, iy2, $ ifill, icolor, ier) if (ixNa .le. ixOI) then write(*,*) $ 'warning: Na position <= 5577 position' write(*,*) 'Na, 5577: ', ixNa, ixOI disp = 0. else disp = (wNa - w55) / real(ixNa - ixOI) end if write(*,*) 'dispersion = ', disp, $ ' Angstroms/pixel' end if else if (key .eq. 'z' .or. $ key .eq. 'w') then write(*,*) if (disp .gt. 0.) then if (key .eq. 'z') then w = w55 + (x - real(ixOI)) * disp else if (key .eq. 'w') then write(*,*) 'enter observed wavelength: ' read(*,*) w x = real(ixOI) + (w - w55) / disp end if icolor = c_green ix1 = nint(x) + ixoff iy1 = iyspec + iyoff call cfmarkpoint(ix1, iy1, 0, isize, m_vline, $ icolor, ier) write(*,'(''marked pixel '', f7.1, $ '' wavelength '', f10.1)') x, w do i = 1, nlines write(*,'(x,i2,'':'',x,a4,$)') i, label(i) end do write(*,*) write(*,*) $ 'choose line or enter rest wavelength: ' read(*,*) w0 j = nint(w0) if (j .ge. 1 .and. $ j .le. nlines) then w0 = wline(j) write(*,'(''chose '', a4, $ '' rest wavelength = '', f7.2)') $ label(j), w0 end if z = w / w0 - 1. write(*,'(''redshift = '', f10.5)') z else write(*,*) 'warning: no dispersion defined' end if else if (key .eq. 'g') then write(*,*) if (disp .gt. 0.) then do i = 1, nlines write(*,'(x,i2,'':'',x,a4,$)') i, label(i) end do write(*,*) write(*,*) $ 'choose line or enter rest wavelength: ' read(*,*) w0 j = nint(w0) if (j .ge. 1 .and. $ j .le. nlines) then w0 = wline(j) write(*,'(''chose '', a4, $ '' rest wavelength = '', f7.2)') $ label(j), w0 end if write(*,*) 'enter guess redshift: ' read(*,*) z w = w0 * (1.+z) x = real(ixOI) + (w - w55) / disp icolor = c_green ix1 = nint(x) + ixoff iy1 = iyspec + iyoff call cfmarkpoint(ix1, iy1, 0, isize, m_vline, $ icolor, ier) write(*,'(''marked pixel '', f7.1, '' wavelength '', $ f10.1)') x, w else write(*,*) 'warning: no dispersion defined' end if else if (key .eq. 'l') then write(*,*) if (disp .gt. 0.) then if (z .eq. -999.) then write(*,*) 'enter guess redshift: ' read(*,*) z else write(*,*) $ 'using current guess redshift = ', z end if do i = 1, nlines w = wline(i) * (1.+z) x = real(ixOI) + (w - w55) / disp if (x .ge. 1. .and. x .le. real(mx)) then if (type(i) .eq. 'e') then icolor = c_blue else if (type(i) .eq. 'a') then icolor = c_red end if ix1 = nint(x) + ixoff iy1 = iyspec + iyoff call cfmarkpoint(ix1, iy1, 0, isize, m_vline, $ icolor, ier) c str4 = label(i) c call cfmarktext(ix1, iy1-isize/2-5, str4, c $ 0.75, 0., icolor, ier) end if end do else write(*,*) 'warning: no dispersion defined' end if else if (key .eq. '0') then write(*,*) if (disp .gt. 0.) then if (z .eq. -999.) then write(*,*) 'enter guess redshift: ' read(*,*) z else write(*,*) $ 'using current guess redshift = ', z end if w = w55 + (x - real(ixOI)) * disp w0 = w / (1.+z) icolor = c_green ix1 = nint(x) + ixoff iy1 = iyspec + iyoff call cfmarkpoint(ix1, iy1, 0, isize, m_vline, $ icolor, ier) write(*,'(''marked pixel '', f7.1, $ '' wavelength '', f10.1)') x, w write(*,'(''z = '', f10.5, $ '' rest wavelength = '', f10.1)') $ z, w0 else write(*,*) 'warning: no dispersion defined' end if else if (key .eq. 'h') then write(*,*) write(*,*) 'enter guess redshift: ' read(*,*) z write(*,'(''redshift = '', f10.5)') z else if (key .eq. 'b') then write(*,*) write(*,*) 'rereading xcdisp6.in' open(11, file='xcdisp6.in', status='old') read(11,'(a40)') maskname write(*,'(a7,a40)') 'mask = ', maskname mnamesave = maskname read(11,'(a40)') axffile read(11,'(a40)') directim write(*,'(a15,a40)') 'direct image = ', directim dnamesave = directim read(11,*) z1, z2 read(11,*) izt write(*,*) 'scale limits = ', z1, z2 if (izt .lt. 0 .or. izt .gt. 2) then write(*,*) 'warning: unknown izt = ', izt izt = 2 write(*,*) 'resetting to log10: izt = ', izt else write(*,*) 'izt = ', izt end if read(11,*) ibspec, ibdir write(*,*) 'subtract ', ibspec, $ ' from spectral images' write(*,*) 'subtract ', ibdir, $ ' from direct images' read(11,*) (ibdir4(j),j=1,4) write(*,*) 'subtract ', (ibdir4(j),j=1,4), $ ' from UBgI direct images' read(11,*) ispecshift read(11,*) nxpad, nypad read(11,*) nxbox, nybox write(*,*) 'ispecshift = ', ispecshift write(*,*) 'nxpad, nypad = ', nxpad, nypad write(*,*) 'nxbox, nybox = ', nxbox, nybox read(11,*) isize c c new lines c read(11,*) imxsize, imysize read(11,*) xspecshift, yspecshift read(11,*) nspec write(*,*) write(*,*) 'imxsize, imysize = ', imxsize, imysize write(*,*) 'xspecshift, yspecshift = ', xspecshift, yspecshift write(*,*) 'nspec = ', nspec close(11) else if (key .eq. '?') then call printhelp else write(*,*) write(*,*) 'unknown command: key, x, y = ', $ key, x, y end if end if c (ier check) end do c (cursor read loop) 999 continue C ----------------- C Clean up and exit C ----------------- call cfclose (ier) end c ************************************************************ subroutine printhelp( ) print '(" Command Summary")' print '(" ")' print '(" a - enter desired aperture")' print '(" n - next aperture")' print '(" p - previous aperture")' print '(" r - re-draw current aperture")' print '(" b - re-read input parameter file xcdisp6.in")' print '(" i - print info on current")' print '(" s - set 5577 and NaD positions")' print '(" t - show current 5577 and NaD positions")' print '(" l - mark lines given guess redshift")' print '(" 0 - get rest wavelength")' print '(" h - set guess redshift")' print '(" z - get redshift from marked position")' print '(" w - get redshift from observed wavelength")' print '(" g - mark line given guess redshift")' print '(" d - delete marker; this is buggy")' print '(" q - quit")' print '(" ? - print help")' print '(" ")' return end c ************************************************************ subroutine getimsect(image, pix, m, nx, ny, $ ix1, ix2, iy1, iy2, $ itx1, itx2, ity1, ity2, ier) character*40 image, errmsg integer im, axleng(7), naxis, ier, isfits integer m, nx, ny, ix1, ix2, iy1, iy2 integer itx1, itx2, ity1, ity2 integer*2 pix(m*m) c write(*,*) 'in getimsect' ier = 0 c abortive attempt to read in fits images c call cfisfits (image, isfits) c if (isfits .eq. 0) then c call cfisiraf (image, ier) call imopen(image, 1, im, ier) if (ier .ne. 0) then write(*,*) write(*,'(a40)') image write(*,*) 'getimsect error: cannot open ' call imemsg(ier,errmsg) write(*, '('' Error: '', a80)') errmsg return end if c write(*,*) 'opened ', image call imgsiz(im,axleng,naxis,dtype,ier) if (ier .ne. 0) then write(*,*) write(*,'(a40)') image write(*,*) $ 'getimsect error: cannot get size information' call imemsg(ier,errmsg) write(*, '('' Error: '', a80)') errmsg return end if c write(*,*) 'naxis = ', naxis if (naxis .ne. 2) then write(*,*) write(*,'(a40)') image write(*,*) 'getimsect error: not 2D image' ier = 1 return end if c write(*,*) 'axleng: ', (axleng(i),i=1,naxis) nx = axleng(1) ny = axleng(2) ix1 = max(ix1, itx1) c ix2 = min(ix2, itx2) ix2 = min(ix2, nx) iy1 = max(iy1, ity1) c iy2 = min(iy2, ity2) iy2 = min(iy2, ny) c image size is itx1:itx2,ity1:ity2 if (ix1 .ge. ix2 .or. iy1 .ge. iy2) then write(*,*) write(*,'(a40)') image write(*,*) $ 'getimsect error: image section no good' write(*,*) 'nx, ny = ', nx, ny write(*,*) 'itx1, itx2, ity1, ity2 = ', $ itx1, itx2, ity1, ity2 write(*,*) 'x: ', ix1, ix2 write(*,*) 'y: ', iy1, iy2 ier = 1 return end if mx = ix2-ix1+1 my = iy2-iy1+1 if (mx*my .gt. m*m) then write(*,*) write(*,'(a40)') image write(*,*) $ 'getimsect error: maximum image size m*m = ', $ m*m, ' exceeded by mx*my = ', mx*my write(*,*) 'mx, my = ', mx, my ier = 1 return end if call imgs2s(im, pix, ix1, ix2, iy1, iy2, ier) if (ier .ne. 0) then write(*,*) write(*,'(a40)') image write(*,*) $ 'getimsect error: cannot read in image section' write(*,*) 'x: ', ix1, ix2 write(*,*) 'y: ', iy1, iy2 call imemsg(ier,errmsg) write(*, '('' Error: '', a80)') errmsg return end if call imclos(im, ier) if (ier .ne. 0) then write(*,*) write(*,'(a40)') image write(*,*) $ 'getimsect error: cannot close' call imemsg(ier,errmsg) write(*, '('' Error: '', a80)') errmsg return end if c write(*,*) 'closed ', image return end c ************************************************************ subroutine displaypix(image, title, pix, m, nx, ny, $ ibitpix, iframe, ifb, $ izscale, z1, z2, izt, ier) character*40 image, title integer m, nx, ny, ibitpix, iframe, ifb integer izscale, izt, ier, nw, nh, lx, ly integer*2 pix(m*m) real a, b, c, d, tx, ty, z1, z2 call cfclearframe(ier) if (ier .ne. 0) then write(*,*) write(*,*) 'displaypix warning: ', $ 'could not clear frame ', iframe end if call cfgetfbconfig(ifb, nw, nh, nf) c write(*,*) 'nw, nh = ', nw, nh c call cfcompzscale(pix, nx, ny, ibitpix, c $ z1, z2) call cfsetztrans(izt) call cfzscaleimage(pix, nx, ny, ibitpix, $ z1, z2) a = 1. b = 0. c = 0. d = -1. tx = (nx/2.) - (nw/2.) + 1. ty = (nh/2.) + (ny/2.) call cfsetwcs(image, title, a, b, c, d, $ tx, ty, z1, z2, izt, ier) if (ier .ne. 0) then write(*,*) write(*,*) 'displaypix warning: ', $ 'could not set wcs for ', $ image end if lx = (nw / 2) - (nx / 2) ly = nh - ((nh / 2) + (ny / 2)) c write(*,*) 'lx, ly, nx, ny = ', lx, ly, nx, ny call cfwritesubras(lx, ly, nx, ny, pix, ier) if (ier .ne. 0) then write(*,*) write(*,*) 'displaypix warning: ', $ 'could not display image ', $ image end if return end c ************************************************************ subroutine dispap(maskname, directim, iap, cenx, ceny, $ dright, dleft, dety, $ z1, z2, izt, ibspec, ibdir, ibdir4, $ nxpad, nypad, nxbox, nybox, $ nspec, firstspec, mx, my, iyspec, $ ispecshift, ixim, iyim, ifb, $ xspecshift, yspecshift) real*4 cenx, ceny, dright, dleft, dety character*40 maskname, directim character*2 ch2(5) data ch2 /'U ', 'BB', 'g ', 'R ', 'I '/ integer*2 ibspec, ibdir, ibdir4(4), ibmed integer iap, izt, nxpad, nypad, nxbox, nybox integer nspec, mx, my, iyspec, ixim, iyim integer ispecshift, xspecshift, yspecshift real z1, z2 character*40 image integer iframe, ifb, izscale, ibitpix integer ixoff, idum integer m, nx, ny, ier, ix1, ix2, iy1, iy2 integer itx1, itx2, ity1, ity2 c parameter (m = 2048) parameter (m = 1024) integer mxsave, mysave integer*2 pix(m*m) integer*2 pix2(m*m), pix3(m*m) real*4 pixr(m*m) integer ndx, ndy, nxoff real*4 w1, w2, ws, disp, yoffset c spectrum location parameters for STIS parameter (w1 = 4400., w2 = 6300., ws = 5176.) parameter (disp = 5.2, yoffset = 68.) ccc usual setting c parameter (disp = 5.2, yoffset = 78.) cc parameter (disp = 5.2, yoffset = -55.) ccc for CFRS-filter "L" masks only real*4 p1, p2, y0 logical firstspec ndx = 25 - 1 ndy = 5 - 1 c nspec+1 is direct images do ii = 1, nspec+1 c for direct images if (ii .eq. nspec+1) then ix1 = nint(cenx) - nxbox ix2 = nint(cenx) + nxbox iy1 = nint(ceny) - nybox iy2 = nint(ceny) + nybox c for spectroscopic images else ix1 = nint( cenx - abs(dright) ) $ - nxpad ix2 = nint( cenx + abs(dleft) ) $ + nxpad y0 = ceny - dety + yoffset p1 = y0 + (w1-ws) / disp p2 = y0 + (w2-ws) / disp iy1 = nint(p1) - nypad + yspecshift iy2 = nint(p2) + nypad + yspecshift end if c for direct images if (ii .eq. nspec+1) then itx1 = 1 itx2 = imxsize ity1 = 1 ity2 = imysize image = directim c for individual spectroscopic images else if (ii .eq. 2 .or. ii .eq. 3) then itx1 = 25 itx2 = imxsize ity1 = 5 ity2 = imysize js = index(maskname, ' ') - 1 if (ii .eq. 2) then image = maskname(1:js)//'1' else if (ii .eq. 3) then image = maskname(1:js)//'2' end if c for coadded spectroscopic image else c ix1 = ix1 - ndx c ix2 = ix2 - ndx c iy1 = iy1 - ndy c iy2 = iy2 - ndy itx1 = 1 itx2 = imxsize ity1 = 1 ity2 = imysize image = maskname end if c write(*,*) c write(*,*) 'input x, y: ', ix1, ix2, iy1, iy2 c write(*,*) maskname c write(*,*) directim c write(*,*) image call getimsect(image, pix2, m, nx, ny, $ ix1, ix2, iy1, iy2, $ itx1, itx2, ity1, ity2, ier) if (ier .ne. 0) then write(*,*) $ 'error: cannot get image section for ', $ image write(*,*) 'x: ', ix1, ix2 write(*,*) 'y: ', iy1, iy2 stop end if c write(*,*) 'output x, y: ', ix1, ix2, iy1, iy2 nr = 0 do ix = ix1, ix2 do iy = iy1, iy2 nr = nr + 1 pixr(nr) = real(pix2(nr)) end do end do call sort(nr, pixr) ratio = real(ibdir4(1))/(pixr(nr*3/4)-pixr(nr/2)) ibmed = nint(pixr(nr/2)) $ - nint(real(ibdir)/ratio) c if (ii .eq. nspec+1) then c ibmed = nint(pixr(nr/2)) - ibdir c else c ibmed = nint(pixr(nr/2)) - ibspec c end if c write(*,'(a20,a21,i4)') c $ image(1:20), ' subtracting ibmed = ', ibmed if (ibmed .lt. 0) then ibmed = 0 c write(*,*) 'resetting ibmed = ', ibmed end if c write(*,*) 'ratio = ', ratio mx = ix2-ix1+1 my = iy2-iy1+1 c write(*,*) 'mx, my = ', mx, my if (ii .eq. 1) then mxsave = mx mysave = my nxoff = mx * nspec + 2*nxbox+1 ixoff = (ix2+ndx) - nint(cenx) + 1 end if if (ii .eq. nspec+1) then c jj = mysave/2 - (2*nybox+1) - 20 c jjj = mysave/2 + 20 jj = mysave/2 - (2*nybox+1) - 10 jjj = mysave/2 + 10 ixim = jjj + nybox iyim = nxbox + 1 do i = 1, 2*nxbox+1 do j = 1, mysave cc if (j .ge. jj+1 .and. cc $ j .le. jj+my .and. cc $ i .le. mx) then cc pix(i+(j-1)*nxoff+nspec*mxsave) cc $ = nint(real(pix2(i+(j-jj-1)*mx)-ibmed)*ratio) cc $ = pix2(i+(j-jj-1)*mx)-ibdir cc else if (j .ge. jjj+1 .and. if (j .ge. jjj+1 .and. $ j .le. jjj+my .and. $ i .le. mx) then pix(i+(j-1)*nxoff+nspec*mxsave) $ = nint(real(pix2(i+(j-jjj-1)*mx)-ibmed)*ratio) cc $ = pix2(i+(j-jjj-1)*mx)-ibdir else pix(i+(j-1)*nxoff+nspec*mxsave) = 0. end if end do end do jr = index(directim, 'R') do ifilt = 1, 5 image = directim(1:jr-1)//ch2(ifilt) call getimsect(image, pix3, m, nx, ny, $ ix1, ix2, iy1, iy2, $ itx1, itx2, ity1, ity2, ier) nr = 0 do ix = ix1, ix2 do iy = iy1, iy2 nr = nr + 1 pixr(nr) = real(pix3(nr)) end do end do call sort(nr, pixr) cc ibmed = nint(pixr(nr/2)) - ibdir ratio = real(ibdir4(1))/(pixr(nr*3/4)-pixr(nr/2)) ibmed = nint(pixr(nr/2)) $ - nint(real(ibdir)/ratio) c write(*,'(a20,a21,i4)') c $ image(1:20), ' subtracting ibmed = ', ibmed if (ibmed .lt. 0) then ibmed = 0 c write(*,*) 'resetting ibmed = ', ibmed end if c write(*,*) 'ratio = ', ratio if (ier .ne. 0) then write(*,*) $ 'error: cannot get image section for ', $ image write(*,*) 'x: ', ix1, ix2 write(*,*) 'y: ', iy1, iy2 stop end if c write(*,*) 'output x, y: ', ix1, ix2, iy1, iy2 if (ifilt .ge. 1 .and. ifilt .le. 3) then jj = mysave/2 - (2*nybox+1) - 10*(3-ifilt) $ - (13-4*ifilt)*nybox/2 else if (ifilt .eq. 4 .or. ifilt .eq. 5) then jj = mysave/2 + 10*(ifilt-3) $ + (4*ifilt-11)*nybox/2 end if do i = 1, 2*nxbox+1 do j = 1, mysave if (j .ge. jj+1 .and. $ j .le. jj+my .and. $ i .le. mx) then pix(i+(j-1)*nxoff+nspec*mxsave) $ = nint(real(pix3(i+(j-jj-1)*mx)-ibmed)*ratio) cc $ = pix3(i+(j-jj-1)*mx)-ibdir4(ifilt) end if end do end do end do c (ifilt = 1, 4) jj = mysave/2 - (2*nybox+1) - 10 jjj = mysave/2 + 10 ixim = jjj + nybox iyim = nxbox + 1 else if (ii .eq. 2 .or. ii .eq. 3) then do i = 1, mx do j = 1, my pix(i+(j-1)*nxoff+(ii-1)*mx) $ = (pix2(i+(j-1)*mx)-ibspec)*2 end do end do else do i = 1, mx do j = 1, my pix(i+(j-1)*nxoff+(ii-1)*mx) $ = pix2(i+(j-1)*mx) end do end do end if end do c (loop over nspec spectra to display) c c rotate image before displaying c c my = nspec*mxsave my = nxoff mx = mysave do i = 1, mx do j = 1, my pix2(i+(my-j+1)*mx) = pix(j+(i-1)*my) end do end do c write(*,*) 'mx, my = ', mx, my ibitpix = 16 izscale = 1 c write(*,*) c write(*,*) 'displaying aperture ', iap call displaypix(image, title, pix2, m, mx, my, $ ibitpix, iframe, ifb, $ izscale, z1, z2, izt, ier) if (firstspec) then c idum = system("imset colormap invert") c idum = system("imset zoom 3") firstspec = .false. end if iyspec = ixoff + 2*nxbox+1 + mxsave*(nspec-1) $ - ispecshift return end c ************************************************************ SUBROUTINE sort(n,arr) INTEGER n,M,NSTACK REAL arr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,ir,j,jstack,k,l,istack(NSTACK) REAL a,temp jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 12 j=l+1,ir a=arr(j) do 11 i=j-1,1,-1 if(arr(i).le.a)goto 2 arr(i+1)=arr(i) 11 continue i=0 2 arr(i+1)=a 12 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 temp=arr(k) arr(k)=arr(l+1) arr(l+1)=temp if(arr(l+1).gt.arr(ir))then temp=arr(l+1) arr(l+1)=arr(ir) arr(ir)=temp endif if(arr(l).gt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp endif if(arr(l+1).gt.arr(l))then temp=arr(l+1) arr(l+1)=arr(l) arr(l)=temp endif i=l+1 j=ir a=arr(l) 3 continue i=i+1 if(arr(i).lt.a)goto 3 4 continue j=j-1 if(arr(j).gt.a)goto 4 if(j.lt.i)goto 5 temp=arr(i) arr(i)=arr(j) arr(j)=temp goto 3 5 arr(l)=arr(j) arr(j)=a jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in sort' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END c ************************************************************