;+
; NAME:
;        sb_disp2im
; PURPOSE:
;        Calculates the pixel index in an input image for each pixel
;        of the display.  Applies de-rotation by assuming that image
;        is on surface of a sphere.
; CATEGORY:
;        image processing
; CALLING SEQUENCE:
;        sb_disp2im,wcs,pix_im,pix_dis,disp=cimd[i].disp,obs=owcs
; INPUTS:
;        wcs : WCS structure of image (processed via sb_wcs)
; KEYWORDS (INPUT):
;        disp   = disp substructure of cimd (see sb_initcom)
;        obswcs = wcs of the observer (only using position and time tags)
;        Note: disp and obswcs keyword parameter required!
; OUTPUTS:
;        pix_im  : vector of pixel indices in image
;        pix_dis : vector of pixel indices in display
; KEYWORDS (OUTPUT):
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;        Many, including
;        - Only for 2D images
;        - Expects crval in either arcsec, deg, or solRad!
;          Other units give nonsensical results
;        - wcs.position.hgln_obs (and some others) must be present
; PROCEDURE:
; MODIFICATION HISTORY:
;        JPW, 6/1/2006
;        JPW, 2/27/2007 simple minded implementation of /sc and /helio
;                       proper way would be via expanded disp structure
;-

PRO sb_disp2im,wcs,ipix,dpix,disp=disp,obswcs=owcs ,sc=sc,helio=helio,swcs=swcs

; 0. Initialize output
ipix = [-1L]
dpix = [-1L]
solrad = 6.9599d8   ; m

; 1. Standardize input wcs (CD variant, cunit in arcsec, etc.)
;    and turn disp into a wcs as well
iwcs = sb_wcs(wcs)
;dwcs = sb_wcs(disp=disp,obswcs=owcs)
dwcs = sb_wcs(disp=disp,obswcs=owcs,helio=helio,sc=sc,swcs=swcs)

; 2. Calculate heliocentric-cartesian coords for display pixels
ndx = dwcs.naxis[0]
ndy = dwcs.naxis[1]
dww = lindgen(ndx*ndy)           ; pixel index vector
d0  = wcs_get_coord(dwcs)        ; d0(2,nx,ny), units in arcsec
a2s = dwcs.position.dsun_obs*!dpi/(solrad*1.8d2*3.6d3)
d0  = d0 * a2s                   ; arcsec2solrad
dx  = reform(d0[0,*,*])
dy  = reform(d0[1,*,*])

; 3. Find on-disk pixels, and project onto sphere
dr2  = dx*dx+dy*dy
ww  = where(dr2 lt 1.0,nww)
if nww gt 0L then begin
   dww = dww[ww]
   dx  = dx[ww]
   dy  = dy[ww]
   dz  = sqrt(1.0-dr2[ww])
endif else return

; 4. Apply B0 rotation to get onto sun's equator
b0 = dwcs.position.solar_b0 / !radeg
d0 = dy
dy = d0 * cos(b0) + dz * sin(b0)
dz = dz * cos(b0) - d0 * sin(b0)

; 5. Correct for observer longitude and solar rotation
a0 = ( iwcs.position.hgln_obs   - dwcs.position.hgln_obs   $
     + iwcs.position.carr_earth - dwcs.position.carr_earth) / !radeg
d0 = dz
dz = d0 * cos(a0) + dx * sin(a0)
dx = dx * cos(a0) - d0 * sin(a0)

; 6. Apply image B0 rotation to get to image observer
b0 = - iwcs.position.solar_b0 / !radeg
d0 = dy
dy = d0 * cos(b0) + dz * sin(b0)
dz = dz * cos(b0) - d0 * sin(b0)

; 7. Find visible pixels
ww = where(dz gt 0.0,nww)
if nww gt 0L then begin
   dww = dww[ww]
   dx  = dx[ww]
   dy  = dy[ww]
endif else return

; 8. Get image pixel coordinates
s2a = solrad*1.8d2*3.6d3/(iwcs.position.dsun_obs*!dpi)
d0  = transpose([[dx],[dy]]) * s2a                    ; solrad2arcsec
pxy = long(wcs_get_pixel(iwcs,d0) + 0.5)

; 9. Find valid pixels and calculate pixel indices
ww  = where(pxy[0,*] ge 0L and pxy[0,*] lt iwcs.naxis[0] and $
            pxy[1,*] ge 0L and pxy[1,*] lt iwcs.naxis[1],nww)
if nww gt 0L then begin
   dpix = dww[ww]
   ipix = reform(pxy[0,ww] + iwcs.naxis[0] * pxy[1,ww])
endif else return

END