function scc_hi_diffuse,hdr,ipsum=ipsum
;+
;
; $Id: scc_hi_diffuse.pro,v 1.2 2013/10/24 15:34:42 crothers Exp $
;
; Project: STEREO-SECCHI
;
; Name: scc_hi_diffuse
;
; Purpose: Compute correction for diffuse sources arrising from changes in the solid angle in the optics
;
; Category: STEREO, SECCHI, HI
;
; Explanation: In the mapping of the optics the area of sky viewed is not equal off axis
;
; Syntax: image=image*scc_hi_diffuse(hdr[,ipsum=n])
;
; Inputs: hdr - the secchi header for the input image
;
; Keywords: ipsum - allows override of header ipsum value for use in L1 and beyond images
;
; Calls: none
;
; Common: None
;
; Restrictions: None 
;
; Side effects: None
;
; Prev. Hist.: None
;
; History: Written 16 October 2013 Jackie Davies
;
; $Log: scc_hi_diffuse.pro,v $
; Revision 1.2  2013/10/24 15:34:42  crothers
; removed scaling of cdelt from calibrated images since this had already been done
;
; Revision 1.1  2013/10/22 14:17:34  crothers
; Diffuse source correction code
;
;
;-

; extract distortion factor and plate scale from header if set, if not use defaults according to
; Brown, Bewsher & Eyles, Solar Phys (2009) 254: 185-335 table 2

; extract on board sum factor

if not keyword_set(ipsum) then ipsum=hdr.ipsum
summing=2^(ipsum-1)

; extract distortion coefficent mu, plate scale and focal length

if hdr.ravg ge 0 then begin
    mu=hdr.pv2_1
    cdelt=hdr.cdelt1*!dtor
endif else begin
    if hdr.detector eq 'HI1' then begin
        if hdr.OBSRVTRY eq 'STEREO_A' then begin
           mu=0.102422
           cdelt=35.96382/3600*!dtor*summing
        endif else if hdr.OBSRVTRY eq 'STEREO_B' then begin
           mu=0.095092
           cdelt=35.89977/3600*!dtor*summing
        endif
    endif else if hdr.detector eq 'HI2' then begin
        if hdr.OBSRVTRY eq 'STEREO_A' then begin
           mu=0.785486
           cdelt=130.03175/3600*!dtor*summing
        endif else if hdr.OBSRVTRY eq 'STEREO_B' then begin
           mu=0.68886
           cdelt=129.80319/3600*!dtor*summing
        endif
    endif
endelse

; from above compute pixel size in mm and paraxial focal length

pixelSize=0.0135*summing
fp=pixelSize/cdelt

; compute linear distance from centre of ccd

x=indgen(hdr.naxis1)-hdr.crpix1+hdr.dstart1
x=rebin(x,hdr.naxis1,hdr.naxis2)
y=indgen(hdr.naxis2)-hdr.crpix2+hdr.dstart2
y=rebin(reform(y,1,hdr.naxis2),hdr.naxis1,hdr.naxis2)
r=sqrt(x*x+y*y)*pixelSize

; get solid angle

gamma=fp*(mu+1.0)/r
cosalpha1=(-1.0*mu+gamma*sqrt(1.0-mu*mu+gamma*gamma))/(1.0+gamma*gamma)
; cosalpha2=(-1.0*mu-gamma*sqrt(1.0-mu*mu+gamma*gamma))/(1.0+gamma*gamma)

correct=( (mu+1.0)^2*(mu*cosalpha1+1.0) )/((mu+cosalpha1)^3)

return,correct

end