function image_prof, im, cd, xc, yc, minr, maxr, nrad, verbose=verbose, rad=rad, dr=dr, proferr=proferr, area=area, brad=brad, expmap=expmap,plt=plt, tv=tv, fits=fits, bgd=bgd, zero=zero, noninter=noninter, bim=bim, perr=perr

	if n_params(0) lt 7 then begin
	  print,'prof = image_prof(im, cd, xc, yc, minr, maxr, nrad, rad=rad, proferr=proferr, area=area, brad = brad, expmap = expmap, plt=plt, tv=tv, fits=fits, bdg=bdd, zero=zero, noninter=noninter, bim=bim, perr=perr)'
	  print,'returns counts/arcmin^2 for nrad bins betweem minr and maxr
	  print,'from xc,yc'
	  print,'Radii returned in rad, area at each radius returned in area'
	  print,'h contains image header if a string array, otherwise = crdelt'
          print,'If brad is specified, [brad(0), brad(1)] gives range in arcmin for bgd.'
          print,'If expmap is given, do exposure map corrections:'
          print,'   - only use pixels where expmap > 0'
          print,'   - divide prof by expmap sum, return prof in counts/arcmin^2/s'
          print, 'expmap must be same size and scale as im'
          print, 'proferr return errors on prof
          print, 'If /plt, plot annuli'
          print, 'If /tv, use tvcircle rather than drawcirc for plotting'
          print, 'If /fits, subtract 0.5 from xc,yc before referencing im (and in tvcircle), but not'
          print, '  drawcirc (i.e., assume pixels were plotted with plot,plist.detx,plist.dety)'
          print, 'If brad is not given, background is taken from bgd, bgderr (def. = 0.0)'
          print, 'If bim is given (with pix in cts/s), determine count rate from source region
          print, 'If /zero, make first anullus be circular (i.e., r0=0)'
	  print, 'If /noninter, suppress all querying`
          print, 'If perr, use Gehrel''s approx. to Poissonian errors, i.e., 1+sqrt(N+0.75) '
	  retall
	endif

        if (n_elements(rad) gt 0 and n_elements(dr) gt 0) then begin
          nrad = n_elements(rad)
          minr = rad(0) - dr(0)
          maxr = rad(nrad-1) + dr(nrad-1)
        endif else begin
	  if keyword_set(verbose) then print, "Computing radii"
          stepsize = 1.*(maxr - minr)/nrad
          rad = (findgen(nrad) + 0.5)*stepsize
          dr = replicate(1., nrad)*stepsize/2
        endelse

        if keyword_set(verbose) then begin
          print,'minr, maxr, nrad = ',minr, maxr, nrad
        endif

        fitscor = keyword_set(fits)
        if keyword_set(verbose) then print, 'Adjusting (xc, yc) by ' + $
          strn(fitscor)
        x = xc-fitscor
        y = yc-fitscor

	sz = size(cd)
	if sz(sz(0)+1) eq 7 then cdelt = abs(sxpar(cd,'cdelt1')) else $
	  cdelt = cd
	cdelt = cdelt*60.		; get arcmin/pix
        if keyword_set(verbose) then print,'Plate scale: ' + strn(cdelt) + "'/pix"

; Exposure map?
        eflag = n_elements(expmap) gt 0

; Background image?
        bflag = n_elements(bim) gt 0; and n_elements(bgd) eq 0

; First check that the profile will fit into the supplied image
	sz = size(im)
	if n_elements(brad) gt 0 then npix = round(max([maxr,brad])/cdelt) $
          else npix = round(maxr/cdelt) ; no. pixels to outer radius
;        if 2*npix+1 gt sz(1) then npix = (sz(1)-1)/2
;        if 2*npix+1 gt sz(2) then npix = (sz(2)-1)/2
        x0 = x - npix
        x1 = x + npix
        y0 = y - npix
        y1 = y + npix
        szx = x1 - x0
        szy = y1 - y0
	if (x0 lt 0 or y0 lt 0 or x1 gt sz(1)-1 or y1 gt sz(2)-1) then begin
	  print,'image_prof: warning... profile and/or bgd. exceeds image boundaries'
          if not eflag then begin
            print, 'WARNING: since no exposure map is supplied, an annulus will probably be a bad'
            print, 'approximation to the shape of the outer bins and/or bgd.'
          endif
          
          x0 = x0 > 0
          y0 = y0 > 0
          x1 = x1 < sz(1)-1
          y1 = y1 < sz(2)-1
        endif
; Now check to see if profile extends beyond exp. map, if nec.
        if eflag then begin
          esz = size(expmap)
          if x1 ge esz(1) then begin
            print, 'Warning: the maximum x-coordinate is being truncated from ' $
              + strn(x1) + ' to ' + strn(esz(1))
            print, 'due to a cut-off in the coverage of the exposure map'
            x1 = esz(1) - 1
          endif
          if y1 ge esz(2) then begin
            print, 'Warning: the maximum x-coordinate is being truncated from ' $
              + strn(y1) + ' to ' + strn(esz(2))
            print, 'due to a cut-off in the coverage of the exposure map'
            y1 = esz(2) - 1
          endif
          estamp = expmap(x0:x1, y0:y1)
          thresh = max(estamp)*0.1
          if keyword_set(verbose) then print, $
            'Exposure map threshold = ' + strn(thresh) + ' s'
        endif
        stamp = im(x0:x1, y0:y1)
	ssz = size(stamp)
        if bflag then bstamp = bim(x0:x1, y0:y1)

; Now stamp (and estamp, bstamp) contain properly limited "valid" regions for
; the entire profile

; Reset npix to max. dist.
        npix = max([x1 - x, x - x0, y1 - y, y0 - y])

; Create image to map distances to (x, y) in stamp
        dist_circle,circle, max([ssz(1), ssz(2)]), x - x0, y - y0
; Make sure circle has exact dim. of stamp (and estamp)
        circle = circle(0:ssz(1) - 1, 0:ssz(2) - 1)
        circle = circle*cdelt   ; now dist in arcmin

	psfout = fltarr(nrad) & area = psfout
        proferr = psfout
        bgderr = 0.
        if n_elements(bgd) eq 0 then bgd = 0.0
        if !d.name eq 'X' then color=100 else color=2
        if n_elements(brad) eq 2 then begin 
          w = where(circle ge brad(0) and circle le brad(1), barea)
          if keyword_set(plt) then begin
            if keyword_set(tv)  then begin
              tvcircle, brad(0)/cdelt, x, y, 100
              tvcircle, brad(1)/cdelt, x, y, 100
            endif else begin
              drawcirc, brad(0)/cdelt, xc, yc, line=2, color=color
              drawcirc, brad(1)/cdelt, xc, yc, line=2, color=color
            endelse
          endif
          if barea eq 0 then print, 'Warning, no pixels in bgd area' else begin
            bgdcnts = total(stamp(w))
            print, strn(bgdcnts) + ' bgd. counts'
            bgd = bgdcnts/barea/cdelt/cdelt
            if keyword_set(perr) then $
              bgderr = (1. + sqrt(bgdcnts + 0.75))/barea/cdelt/cdelt else $
              bgderr = sqrt(bgdcnts)/barea/cdelt/cdelt
            print, 'bgd = ' + strn(bgd) + ' counts/arcmin^2'
            print, 'bgd err = ' + strn(bgderr) + ' counts/arcmin^2'
            ans = ''
            if not eflag and not keyword_set(noninter) then $
              read,'Plot background area (y/n)? ',ans
            if ans eq 'y' then begin
              wy = w/ssz(1)
              wx = w-wy*ssz(1)
              oplot,wx + x0, wy + y0, psym=3
            endif
            if eflag then begin
              w = w(where(estamp(w) gt thresh, barea))
              if barea eq 0 then begin
                print,'No non-zero exp. map pixels in bgd area, setting bgd to zero'
                bgd = 0.0
              endif else begin
                bexp = total(estamp(w))/barea
                bgd = bgd/bexp
                bgderr = bgderr/bexp
                print,'bgd = '+strn(bgd)+' counts/s/arcmin^2'
                print,'bgd err = '+strn(bgderr)+' counts/s/arcmin^2'
                ans = ''
                if not keyword_set(noninter) then $
                  read,'Plot background area (y/n)? ',ans
                if ans eq 'y' then begin
                  wy = w/ssz(1)
                  wx = w - wy*ssz(1)
                  oplot,wx + x0, wy + y0, psym=3
                endif
              endelse
            endif
          endelse
        endif

        if bflag then begin
          w = where(circle le maxr, n)
          if eflag then w = w(where(estamp(w) gt thresh, n))
          bgd = total(bstamp(w))/n/cdelt/cdelt
          print, '(blank sky) bgd = ' + strn(bgd) + ' counts/s/arcmin^2'
;	  stop
        endif

        if keyword_set(verbose) then begin
          print,format ='$(a,t10,a,t20,a,t30,a,t45,a,t60,a,t70,a)', 'rad low', $
            'rad hi', 'area', 'area', 'pi*(r2^2-r1^2)', 'exposure', 'counts'
          print,format ='$(a,t10,a,t20,a,t30,a,t45,a,t60,a)', '(arcmin)', $
            '(arcmin)','(pixels)','(arcmin^2)','(arcmin^2)','(secs)'
        endif

        for i=0,nrad-1 do begin
	  low = rad(i) - dr(i)
          if keyword_set(zero) and i eq 0 then low = 0.0
          hi = rad(i) + dr(i)
	  w = where(circle ge low and circle lt hi,c)
          if keyword_set(plt) then begin
            if keyword_set(tv)  then begin
              tvcircle, low/cdelt, x, y, 100
              tvcircle, hi/cdelt, x, y, 100
            endif else begin
              drawcirc, low/cdelt, xc, yc, color=color+1
              drawcirc, hi/cdelt, xc, yc, color=color+1
            endelse
          endif
          texp = 1.0
	  if c gt 0 and eflag then begin
            ew = where(estamp(w) gt thresh,c)
            if c gt 0 then w = w(ew) else stop
            if keyword_set(plt) then begin 
              wy = w/ssz(1)
              wx = w - wy*ssz(1)
              oplot, wx + x0, wy + y0, psym=3, color=(i mod 6) + 1
            endif
            texp = total(estamp(w))/c
;            print,'c, texp = ', c, texp
          endif
;         ww = where(stamp(w) gt 0)
;          yy = w/sz
;          xx = w - yy*sz
;          oplot, xx, yy, psym=3
          if c gt 0 then begin
            effarea = c*cdelt^2*texp
            bcnts = bgd*effarea
            cnts = total(stamp(w))
            psfout(i) = (cnts-bcnts)/effarea
;            proferr(i) = sqrt(total(stamp(w))+bcnts)/effarea
            if keyword_set(perr) then sigma = (1. + sqrt(cnts + $
              0.75))/effarea else sigma = sqrt(cnts)/effarea
            proferr(i) = sqrt(sigma^2 + bgderr^2)
            area(i) = effarea
	    if keyword_set(verbose) then $
              print,'$(f5.2,t10,f5.2,t20,i7,t30,g12.4,t45,g12.4,t60,f10.0,t70,i6)', $
                low,hi,c,c*cdelt^2,!pi*(hi^2-low^2),texp,cnts
;  else print, total(stamp(w)), bcnts, effarea/texp, texp
          endif else $
	    print,'Warning - no counts between r = ',strn(low),' and r = ',$
	      strn(hi)
	endfor
	return,psfout
	end