Hi,
The following is the procedure for computation of legendre
moments. I am trying to compile the examples indicated in the
procedure, but they were showing that the arguments of LEGENDRE shall
be equal to or greater than 0. Normally, the arguments will be between
-1 and 1 in my case.
Can anyone try to run these example in the procedure and let me know
whether this procedure work or not? If this is showing the same errors
as mine, please suggest me the rectification for the same.
------------------------------------------------------------ ------------------------------------------------------------ -----
pro pmom,mu,phase,nm,pm,view=view
;+
; ROUTINE: pmom
;
; PURPOSE: compute legendre moments of scattering phase function
;
; USEAGE: pmom,mu,phase,nm,pm
;
; INPUT:
; mu cosine of polar angle, in range -1, 1
; phase phase function at each polar angle, normalization is
arbitrary
;
; NOTE: the accuracy of the integration depends on the
; number values provided for mu and phase. For
; example, 10000 steps were required to obtain better
; than 1% accuracy for all 21 moments of the
; Henyey-Greenstein function with g=0.7. More steps
; may be required when higher moments or more sharply
; peaked functions are considered.
;
; nm number of of legendre moments. moments are numbered from 0
to nm
; for a total of nm+1 values.
;
; OUTPUT:
; pm legendre moments (nm+1) values
;
; DISCUSSION:
; This procedure generates the legendre moments of an arbitrary
; function defined in the range -1 to 1. These coefficients are
; the input quantities required by DISORT for the specification of
; of scattering phase function or BRDF. The output quantity, pm,
; is defined as,
;
; / / /
; pm(i) = | f(mu) P(i,mu) d mu / | f(mu) d mu
; / / /
;
; where P(i,mu) is the legendre polynomial, f is the phase function,
and
; the range of the integrals are from -1 to 1.
;
;
; LIMITATIONS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; EXAMPLE: find moments of Henyey-Greenstein phase function
; The analytic result has pmom=1/g^indgen(nm+1)
;
; ns=10000 & g=.7 & nm=20
;
; mu=findrng(-1,1,ns)
; phase=(1.-g^2)/(1.+g^2-2*g*mu)^1.5
; pmom,mu,phase,nm,pm,/view
; print,f='(a/(10f8.4))','pmom',pm(1:*)
; print,f='(a/(10f8.4))','pmom/g^i',pm(1:*)/g^(1+indgen(nm))
;
;
; EXAMPLE: find moments of Rayleigh scattering phase function
; the analytic result has pmom(0)=1 and pmom(2)=.1
;
; ns=1000 & nm=10
;
; mu=findrng(-1,1,ns)
; phase=(1.+mu^2)
; pmom,mu,phase,nm,pm
; print,f='(a/10f8.4)','pmom',pm
;
; phs=phase<0.
; norm=integral(mu,phase)
; for i=0,nm do phs=phs+.5*norm*(2*i+1)*pm(i)*legendre(i,mu)
; plot,mu,phase
; oplot,mu,phs
;
; AUTHOR: Paul Ricchiazzi 29 Jul 98
; Institute for Computational Earth System Science
; University of California, Santa Barbara
; paul@icess.ucsb.edu
;
; REVISIONS:
;
;-
;
if not keyword_set(view) then view=0
pm=fltarr(nm+1)
norm=integral(mu,phase)
phs=phase<0.
for i=0,nm do begin
p=legendre(i,mu)
pm(i)=integral(mu,phase*p)/norm
if view then phs=phs+.5*norm*(2*i+1)*pm(i)*p
endfor
if view then begin
colors=[!p.color,.5*!d.n_colors]
plot_io,mu,phase,title='phase function'
legend,'Phase Function\\Legendre Fit',color=colors,spos='tl'
oplot,mu,phs,color=.5*!d.n_colors
endif
return
end
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
