The Heney-Greenstein function p_HG=(1-g*g)/(1+g*g-2*g*mu)**(3/2) is often used to approximate scattering phase functions of aerosols or clouds. Show that the Legendre coefficients of p_HG are P_l=g^l. Investigate how many Legendre coefficients are required for asymmetry parameters of 0, 0.5, 0.7, 0.8, 0.85, 0.9.

Hint: Use the tool phase, which can be found in libRadtran/bin to calculate the phase function from the Legendre polynomials. Help is obtained by

 phase -h
 

We decided to use a graphical investigation of the problem. We need to find the number of required Legendre coefficients so that the two curves of P_HG and P_l converge. For this we write a shell script in combination with a python plot routine. Afterwards we look at the output to decide which number of Legendre coefficients is sufficient for individual asymmetry parameters.

Shell script:

rte_heneygreenstein.sh

#!/bin/sh -f
# shell skript
 
libradtran_path='/local/libRadtran-1.5-beta/bin/'
tmp_path='/home/users/merk/radtran_2010/rte/tmp/'
script_path='/home/users/merk/radtran_2010/rte/'
 
 
# create files for different g and l
echo 'create files for different g and l'
 
cd $script_path
 
 
# vary g and l to find a graphical solution
 
for g in 0.9 #0.7 0.8 0.85 0.9
   do echo 'g: ' $g
   for moments in 70 71 72 73 74 75 #10 20 50 70
	do python rte_heneygreenstein.py $moments $g 		# call python-program (param1=l, param2=g)
	echo '\n\n'
   done
done
 
###echo 'done. files created. \n \n'
 
echo 'run libradtran -> phase'
cd $libradtran_path
 
for g in 0.9 #0.7 0.8 0.85 0.9
   do echo 'g: ' $g
   for moments in 70 71 72 73 74 75 #10 20 50 70
	do ./phase -c $tmp_path"rte_g"$g"_moments"$moments"_.dat" > $tmp_path"phg_g"$g"_moments"$moments"_.dat"
	echo '\n\n'
   done
done
 
 
cd $script_path
echo 'make plots'
 
for g in 0.9 #0.7 0.8 0.85 0.9
   do echo 'g: ' $g
   for moments in 70 71 72 73 74 75 # 10 20 50 70
	do python rte_heneygreenstein_plot.py $moments $g 
	echo '\n\n'
   done
done

rte_heneygreenstein.py

from pylab import *
 
def moment(g,l):
	moment=pow(g,l)
	return moment
 
 
nrofmoments = int(sys.argv[1]) # uebergabeparameter 1
print('l=' + str(nrofmoments))
 
g = float(sys.argv[2]) # uebergabeparameter 2
print('g=' + str(g))
 
l=arange(0,nrofmoments,1) # momente 0 bis nrofmoments in 1er schritten
out = moment(g,l)
 
filename='tmp/rte_g'+str(g)+'_moments'+str(nrofmoments)+'_'+'.dat'
 
 
print 'open file\n write file ...'
 
f=open(filename, 'w')
for j in out:
	f.write(str(j)+'\n')
f.close
 
print 'close file ...'
 
 
 
exit()

rte_heneygreenstein_plot.py

from pylab import *
from matplotlib import pyplot
 
 
def pHG_analytical(mu,g):
	phg=(1.0-g*g)/(pow((1.0+g*g-2.0*g*mu),3.0/2.0))
	return phg
 
 
# parameter
nrofmoments = int(sys.argv[1]) # uebergabeparameter 1
print('l=' + str(nrofmoments))
 
g = float(sys.argv[2]) # uebergabeparameter 2
print('g=' + str(g))
 
# save or show? ### dont use show if run by shell-script!!!
save=1 # to show set 0, to save set 1
 
filename='tmp/phg_g'+str(g)+'_moments'+str(nrofmoments)+'_'+'.dat'
img_filename='tmp/phg_g'+str(g)+'_moments'+str(nrofmoments)+'_'+'.png'
 
print 'plot ...'
 
data=loadtxt(filename)
mu=data[:,0]
phg=data[:,1]
 
pyplot.plot(mu,phg,label='libradtran')
pyplot.plot(mu,pHG_analytical(mu,g),label='analytical')
 
pyplot.title('g='+str(g)+'    # moments='+str(nrofmoments))
pyplot.xlabel('mu=cos(theta)')
pyplot.ylabel('pHG')
 
if save==0:
   pyplot.show()
 
if save==1:
   print 'save plot ...'
   pyplot.savefig(img_filename)
 
 
 
exit()

solution:
g=0 - 1 moment
g=0.5 - 9 moments

g=0.7 - 19 moments

g=0.8 - 27 moments

g=0.85 - 43 moments

g=0.9 - 75 moments