#!/usr/bin/python
from math import sqrt

## __________ constants go here __________
m_e         =  0.000510998902               ## from fg.h (micheld)
m_mu        =  0.105658357                  ## "                 "
pi_md       =  3.14159265358979323846       ##as used in micheld
GF2         = 969.797746147
#GF2         =  940.5077 #er, brute force tuned by hand
#GF2         =   #er, brute force tuned by hand

## __________ preliminary calculations go here __________
Wem         = (m_mu**2+m_e**2)/(2*m_mu)
x0          = m_e/Wem
k_const     = (m_mu/4*(pi_md**3))*(Wem**4)*(GF2**2)

def checkx(x): return (x**2-x0**2)

def k(x): return k_const*sqrt(x**2-x0**2)
def FIS(x,mp): return x*(1-x)+(2./9.)*mp['rho']*(4*x*x-3*x-x0**2)+mp['eta']*x0*(1-x)
def FAS(x,mp): return (1./3.)*mp['xi']*sqrt(x**2-x0**2)*(1-x+(2./3.)*mp['delta']*(4*x-3+(sqrt(1-x0**2)-1)))

def approx_michel_spectrum(x,cos_theta,
                           mp={'rho':0.75,'eta':0.00,'xi':1.00,'delta':0.75}):
    if checkx(x)>0:
        return k(x)*(FIS(x,mp)-cos_theta*FAS(x,mp))
    else:
        return 0.
#    return k(x)*(FIS(x)+mp['xi']*cos_theta*FAS(x))

def Nu(x,a,b,mp={'rho':0.75,'eta':0.00,'xi':1.00,'delta':0.75}):
    if checkx(x)>0.:
        return k(x)*(FIS(x,mp)*(b-a)-0.5*(b**2-a**2)*FAS(x,mp))
    else:
        return 0
def Nd(x,a,b,mp={'rho':0.75,'eta':0.00,'xi':1.00,'delta':0.75}):
    if checkx(x)>0.:
        return k(x)*(FIS(x,mp)*(b-a)+0.5*(b**2-a**2)*FAS(x,mp))
    else:
        return 0

def asy(x,a,b,mp={'rho':0.75,'eta':0.00,'xi':1.00,'delta':0.75}):
    up = Nu(x,a,b,mp)
    do = Nd(x,a,b,mp)
    if up+do!=0.:
        return (do-up)/(up+do)
    else:
        return 0.