; Compute partial coding vs. solid angle
;
; This is a script, not a procedure, so input/output is simply done
; with IDL variables on the command line.
;
; Input
;   IMG - 2D IDL image contaning the partial coding map
;   HH  - IDL string array containing the header of the FITS file
;   (an example is provided below)
;
; Output
;   PP - array of partial coding values
;   OHM - array of solid angles (sr) corresponding to PP
; 
; C. Markwardt 2004 Aug 12
;

PCODERANGE = [0d, 1.05d]   ;; Range of partial coding fractions to test

; Read IMG if none is present 
if n_elements(img) EQ 0 then fxread, 'pcodemap.img', img, hh
; Convert FITS WCS header keywords to IDL structure
if n_elements(astr) EQ 0 then extast, hh, astr

;; Dimensions
sz = size(img)
nx = sz(1)
ny = sz(2)

;; Start in FITS pixel coords
i = dindgen(nx)+1d & j = dindgen(ny) + 1d
ii =  i # (j*0+1) 
jj =  (i*0+1) # j 

;; BAT image coordinates  (straight linear transformation)
imx = (ii-astr.crpix(0))*astr.cdelt(0) + astr.crval(0)
imy = (jj-astr.crpix(1))*astr.cdelt(1) + astr.crval(1)

;; Differential elements in IMX and IMY
dimx = astr.cdelt(0) & dimy = astr.cdelt(1)
dimx = abs(dimx) & dimy = abs(dimy)

;; Differential element of solid angle
dohm = dimx * dimy /(1d + imx*imx + imy*imy)^1.5d


;; Loop over partial coding fractions
pp = rvspan(pcoderange, 105)
ohm = pp*0
for i = 0, n_elements(pp)-1 do begin
    wh = where(img GE pp(i), ct)
    if ct EQ 0 then continue

    ;; Perform the integral over solid angle
    ohm(i) = total(1.0 * dohm(wh))
end


end