#Program to test ideas about the asymmetry factor

import math
from ClimateUtilities import *

def P(muRel):
    return math.exp(2.*(1+muRel)*(1+muRel))

def Pcos(muRel):
    return .5*muRel*P(muRel)

m = romberg(P)
norm = m([-1.,1.])/2.
m = romberg(Pcos)
ghat = m([-1.,1.])/norm
print 'Asymmetry factor',ghat

def H(mu):
    if mu>=0:
        return 1.
    else:
        return -1.

def f(phip,args):
    costhetap = args[0]
    costheta = args[1]
    sintheta = math.sqrt(1.-costheta**2)
    sinthetap = math.sqrt(1.-costhetap**2)
    mu = costheta*costhetap +\
      sintheta*sinthetap*math.cos(phip)
    return (1.*P(mu)+-1.*P(-mu))/(norm*4.*math.pi)

def intphi(costhetap,costheta):
    m = romberg(f)
    return m([0.,2.*math.pi],[costhetap,costheta])

def g(costheta):
    m = romberg(intphi)
    return m([0.,1.],costheta)

def intg1(costheta):
    return g(costheta)

m = romberg(intg1)
print 'Hemi asymm',m([0.,1.])

costhetaL = [-1. + .1*i for i in range(21)]
gL = [g(costheta) for costheta in costhetaL]
gCos = [gL[-1]*costheta for costheta in costhetaL]
c = Curve()
c.addCurve(costhetaL)
c.addCurve(gL)
c.addCurve(gCos)
plot(c)