Viewing contents of file '../idllib/contrib/buie/bidr2.pro'
;+
; NAME:
; bidr2
; PURPOSE: (one line)
; Compute the bi-directional reflectance (newer Hapke formula). This
; is coded from equation 12.55 on page 346 in Hapke's book, "Theory of
; Reflectance and Emittance Spectroscopy".
; DESCRIPTION:
; CATEGORY:
; Miscellaneous
; CALLING SEQUENCE:
; ans = bidr2(w,emu,imu,g,holes,pzero,b0,theta)
; INPUTS:
; w - Single scattering albedo.
; emu - Cosine of the emission angle.
; imu - Cosine of the incidence angle.
; g - Phase angle, in radians.
; holes - Compaction parameter value (1986 formalism).
; pzero - Value of the single particle phase function.
; theta - Surface roughness value. (radians)
; OPTIONAL INPUT PARAMETERS:
; None.
; KEYWORD PARAMETERS:
; None.
; OUTPUTS:
; Return value is the bi-directional reflectance.
; COMMON BLOCKS:
; None.
; SIDE EFFECTS:
; None.
; RESTRICTIONS:
; Any input may be a vector. If more than one is a vector then the
; lengths must match. The return will have the same dimensions as
; the input.
; PROCEDURE:
; MODIFICATION HISTORY:
; Written by Marc W. Buie, Lowell Observatory, 1997/08/21
; 97/09/18, MWB, added surface roughness parameter
;-
function bidr2,in_w,in_emu,in_imu,in_g,in_holes,in_pzero,in_b0,in_theta,H93=h93
check=[n_elements(in_w),n_elements(in_emu),n_elements(in_imu), $
n_elements(in_g),n_elements(in_holes),n_elements(in_pzero), $
n_elements(in_b0),n_elements(in_theta)]
h93 = keyword_set(h93)
tlen = max(check)
z=where(check ne 1 and check ne tlen,count)
IF count ne 0 THEN BEGIN
print,'BIDR2: Error, lengths of inputs must match or be 1.'
return,0.0
ENDIF
; Promote all inputs to same length
IF n_elements(in_w) eq 1 THEN w = replicate(in_w, tlen) ELSE w = in_w
IF n_elements(in_emu) eq 1 THEN emu = replicate(in_emu, tlen) ELSE emu = in_emu
IF n_elements(in_imu) eq 1 THEN imu = replicate(in_imu, tlen) ELSE imu = in_imu
IF n_elements(in_g) eq 1 THEN g = replicate(in_g, tlen) ELSE g = in_g
IF n_elements(in_holes) eq 1 THEN holes = replicate(in_holes,tlen) ELSE holes = in_holes
IF n_elements(in_pzero) eq 1 THEN pzero = replicate(in_pzero,tlen) ELSE pzero = in_pzero
IF n_elements(in_b0) eq 1 THEN b0 = replicate(in_b0, tlen) ELSE b0 = in_b0
IF n_elements(in_theta) eq 1 THEN theta = replicate(in_theta,tlen) ELSE theta = in_theta
gamma = sqrt(1-w)
bscat=replicate(0.,tlen)
m = where( g lt !pi/2 and b0 ne 0.0, count)
if count ne 0 then begin
bscat[m] = 1.0/(1.0 + 1.0/holes[m]*tan(g[m]/2.0))
end
bscat = b0*bscat
bidr = fltarr(tlen)
; This is the case where the surface roughness term is turned off.
z=where(theta eq 0.0,count)
IF count ne 0 THEN BEGIN
twoemu = emu[z] + emu[z]
twoimu = imu[z] + imu[z]
IF h93 THEN BEGIN
r0 = 2.0/(1.0+gamma[z])-1.0
lnt = replicate(0.0,count)
zn0 = where(emu[z] ne 0.0)
IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+emu[z[zn0]])/emu[z[zn0]])
hobs = 1.0/(1.0 - (1.0-gamma[z])*emu[z]*(r0+(1.0-0.5*r0-r0*emu[z])*lnt))
lnt = replicate(0.0,count)
zn0 = where(imu[z] ne 0.0)
IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+imu[z[zn0]])/imu[z[zn0]])
hsun = 1.0/(1.0 - (1.0-gamma[z])*imu[z]*(r0+(1.0-0.5*r0-r0*imu[z])*lnt))
ENDIF ELSE BEGIN
hobs = (1+twoemu) / (1 + gamma[z] * twoemu)
hsun = (1+twoimu) / (1 + gamma[z] * twoimu)
ENDELSE
bidr[z] = w[z]*imu[z]/(imu[z]+emu[z])*((1+bscat[z])*pzero[z]+hobs*hsun-1)/(4*!pi)
ENDIF
; With surface roughness
z=where(theta ne 0.0,count)
IF count ne 0 THEN BEGIN
emue = emu[z]
imue = imu[z]
emue0 = emu[z]
imue0 = imu[z]
sfun = emu[z]
; ancillary values
tanthe = tan(theta[z])
costhe = cos(theta[z])
cotthe = 1.0/tanthe
cotthe2 = cotthe^2
i = acos(imu[z])
sini = sin(i)
e = acos(emu[z])
sine = sin(e)
cosg = cos(g[z])
cosphi = replicate(1.0,count)
zz = where(i*e ne 0.0,countzz)
if countzz ne 0 then cosphi[zz]=(cosg - imu[z[zz]]*emu[z[zz]])/(sini[zz]*sine[zz])
phi = acos(cosphi)
sinphi2_2=sin(phi/2.0)^2
gold=1.0e-7
z0=where(abs(sini) lt gold,count0)
IF count0 ne 0 THEN sini[z0]=gold
z0=where(abs(sine) lt gold,count0)
IF count0 ne 0 THEN sine[z0]=gold
coti = imu[z]/sini
coti2 = coti^2
cote = emu[z]/sine
cote2 = cote^2
e1i = exp( -2.0/!pi*cotthe*coti ) ; eqn. 12.45b, p. 344
e2i = exp( -1.0/!pi*cotthe2*coti2 ) ; eqn. 12.45c, p. 344
e1e = exp( -2.0/!pi*cotthe*cote )
e2e = exp( -1.0/!pi*cotthe2*cote2 )
chi = 1.0/sqrt(1.0+!pi*tanthe^2) ; eqn. 12.45a, p. 344
fg = exp( -2.0 * tan(phi/2.0) ) ; eqn. 12.51, p. 345
emue0 = chi * ( emu[z] + sine * tanthe * e2e / ( 2.0 - e1e ) )
imue0 = chi * ( imu[z] + sini * tanthe * e2i / ( 2.0 - e1i ) )
; e >= i
zz = where(e ge i,countee)
IF countee ne 0 THEN BEGIN
denom = 2.0 - e1e[zz] - (phi[zz]/!pi)*e1i[zz]
imue[zz] = chi[zz] * ( imu[z[zz]] + sini[zz] * tanthe * $
( cosphi[zz]*e2e[zz] + sinphi2_2[zz]*e2i[zz] ) / denom )
emue[zz] = chi[zz] * ( emu[z[zz]] + sine[zz] * tanthe[zz] * $
( e2e[zz] - sinphi2_2[zz]*e2i[zz] ) / denom )
sfun[zz] = emue[zz]/emue0[zz] * imu[z[zz]]/imue0[zz] * chi[zz] / $
( 1.0 - fg[zz] + fg[zz]*chi[zz]*imu[z[zz]]/imue0[zz] )
ENDIF
; e < i
zz = where(e lt i,countee)
IF countee ne 0 THEN BEGIN
denom = 2.0 - e1i[zz] - (phi[zz]/!pi)*e1e[zz]
imue[zz] = chi[zz] * ( imu[z[zz]] + sini[zz] * tanthe[zz] * $
( e2i[zz] - sinphi2_2[zz]*e2e[zz] ) / denom )
emue[zz] = chi[zz] * ( emu[z[zz]] + sine[zz] * tanthe[zz] * $
( cosphi[zz]*e2i[zz] + sinphi2_2[zz]*e2e[zz] ) / denom )
sfun[zz] = emue[zz]/emue0[zz] * imu[z[zz]]/imue0[zz] * chi[zz] / $
( 1.0 - fg[zz] + fg[zz]*chi[zz]*emu[z[zz]]/emue0[zz] )
ENDIF
twoemu = emue + emue
twoimu = imue + imue
IF h93 THEN BEGIN
r0 = 2.0/(1.0+gamma[z])-1.0
lnt = replicate(0.0,count)
zn0 = where(emue ne 0.0)
IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+emue[zn0])/emue[zn0])
hobs = 1.0/(1.0 - (1.0-gamma[z])*emue*(r0+(1.0-0.5*r0-r0*emue)*lnt))
lnt = replicate(0.0,count)
zn0 = where(imue ne 0.0)
IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+imue[zn0])/imue[zn0])
hsun = 1.0/(1.0 - (1.0-gamma[z])*imue*(r0+(1.0-0.5*r0-r0*imue)*lnt))
ENDIF ELSE BEGIN
hobs = (1+twoemu) / (1 + gamma[z] * twoemu)
hsun = (1+twoimu) / (1 + gamma[z] * twoimu)
ENDELSE
bidr[z] = w[z]*imue/(imue+emue)*((1+bscat[z])*pzero[z]+hobs*hsun-1)/(4*!pi)*sfun
ENDIF
err=check_math()
IF (err AND '1B'X) ne 0 THEN $
print,'BIDR: Math error detected, code ',err
return,bidr
end