====== Legendre decomposition of Heney Greenstein phase function ======

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
   

====== solution (Daniel, Benjamin, Florian) ======

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:

<file bash 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
</file>


<file python 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()

</file>


<file python 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()
</file>


\\
__solution:__\\
g=0    -  1 moment\\
g=0.5  -  9 moments\\
{{:teaching:radiative_transfer:phg_g0.5_moments9_.png|}} \\

g=0.7  -  19 moments\\
{{:teaching:radiative_transfer:phg_g0.7_moments19_.png|}}

g=0.8  -  27 moments\\
{{:teaching:radiative_transfer:phg_g0.8_moments27_.png|}}

g=0.85 -  43 moments\\
{{:teaching:radiative_transfer:phg_g0.85_moments43_.png|}}

g=0.9  -  75 moments\\
{{:teaching:radiative_transfer:phg_g0.9_moments75_.png|}}