pro simrad,prfoz,prftemp,pres,sza,satza,phi,ref,bcw,fwhm,rad
;
; simrad calls mateer to simulate radiances for a finite
; bandwidth using solar flux weighting across a triangular
; slit function with fwhm in angstrom to the nearest half A,
; and band center wavelength in angstrom.
;
iter=[2,9]
rhoflag  = 1
;
;  read in coefficients data file
;
bass=fltarr(6,1601)
file_name='BSBR_V8.DAT'
openr,iunit, file_name, /get_lun
hdr=''
readf,iunit,hdr
readf,iunit,bass
close,iunit
free_lun, iunit
;
;  define set of bass wavelengths surrounding band center wavelength (bcw)
;
intgrl=where(abs(bass(0,*)-bcw) lt fwhm,nint)
wav=reform(bass(0,intgrl))
dlam=abs(wav-bcw)
slit=1.0-abs(dlam)/fwhm
c0=reform(bass(1,intgrl))
c1=reform(bass(2,intgrl))
c2=reform(bass(3,intgrl))
bet=reform(bass(4,intgrl))
rho=reform(bass(5,intgrl))
;
;  convert wav to nm
;
wav=wav/10.0
;
;  call tomrad code to calculate radiances at each wavelength
;
tomrad,pres,sza,satza,phi,ref,wav,c0,c1,c2,bet,rho,prfoz, $
 prftemp,iter,rhoflag,i0s,i1s,i2s,ts,sbs,itots,nval
rads=itots/!pi
;for i=0,22 do print,i,wav(i),i0s(i),itots(i),nval(i)
;
;  get solar flux for weighting
;
sols=fltarr(2,2047)
file_name='SOLSTICE'
openr,iunit, file_name, /get_lun
readf,iunit,sols
close,iunit
free_lun, iunit
flux=interpol(sols(1,*),sols(0,*),wav)
;
; calculate integral in straight forward fashion (this approach gives
; different answers at high sza and non-zero scan angles)
;
numsza      = n_elements(sza)
numsatza    = n_elements(satza)
rad=fltarr(numsza,numsatza)
for isza=0,numsza-1 do begin
for isatza=0,numsatza-1 do begin
rad(isza,isatza)=total(flux*slit*rads(*,isza,isatza))/total(flux*slit)
endfor
endfor
;
;  calculate weighted values for components of rad (as in tables)
;
;numsza      = n_elements(sza)
;numsatza    = n_elements(satza)
i0=dblarr(numsza,numsatza)
i1=dblarr(numsza,numsatza)
i2=dblarr(numsza,numsatza)
t=dblarr(numsza,numsatza)
for isza=0,numsza-1 do begin
for isatza=0,numsatza-1 do begin
i0(isza,isatza)=total(flux*slit*i0s(*,isza,isatza))/total(flux*slit)
i1(isza,isatza)=total(flux*slit*i1s(*,isza,isatza))/total(flux*slit)
i2(isza,isatza)=total(flux*slit*i2s(*,isza,isatza))/total(flux*slit)
t(isza,isatza)=total(flux*slit*ts(*,isza,isatza))/total(flux*slit)
endfor
endfor
sb=total(slit*sbs)/total(slit)
;
frac_ref = ref(0)/100.
fac=frac_ref / (1.-frac_ref*sb)
ground_term = t*fac
ez = i0 + i1*cos(phi(0)*!dtor) + i2*cos(2*phi(0)*!dtor)
itot = i0 + i1*cos(phi(0)*!dtor) + i2*cos(2*phi(0)*!dtor) + ground_term
rad=itot/!pi
;
;help,i0,i1,i2,t,sb,ground_term,itot,rad
;print,'i0/!pi,i1/!pi,i2/!pi,t/!pi,sb,ground_term/!pi,ez/!pi,rad,-100.0*alog10(rad)'
;print,i0/!pi,i1/!pi,i2/!pi,t/!pi,sb,ground_term/!pi,ez/!pi,rad,-100.0*alog10(rad)
;help,i0,t,sb,ez,rad
;print,'i0/!pi,t/!pi,sb,ez/!pi,rad,-100.0*alog10(rad)'
;print,i0/!pi,t/!pi,sb,ez/!pi,rad,-100.0*alog10(rad)
;help,i0,t,sb,ez,rad
;print,'i0,t,sb'
;print,i0,t,sb
;
return
;
end