#!/usr/bin/python

import fg_nla,math,numpy,sys

## __________ constants go here __________
m_e         =  0.000510998902               ## from fg.h (micheld)
m_mu        =  0.105658357                  ## "                 "
rem         =  m_e/m_mu
ame2amu2    =  rem*rem                      ##fortran param 5
alme        =  -math.log(ame2amu2)          ##"           " 6 
pi_md       =  3.14159265358979323846       ##as used in micheld
alpha_qed0  =  (1./137.03599976)/(2*pi_md)  ##fortran param 7
x0          =  2*rem/(1+ame2amu2)           ##fortran param 8
aldp        =  (1-0.)*math.log(rem)

## ____________________
def michel_spectrum(x,cos_theta,
                    mp={'rho':0.75,'eta':0.00,'xi':1.00,'delta':0.75},
                    RCs='all'):

#    if x[0]>1.: return 0.  #this line works for dealing with singularity, but NOT for derivatives

    # choose level of corrections
    if RCs == 'all':
        RC = {'base':1,'Oalpha':1,'Oalpha2L2':1,'Oalpha2L':1,'Oalpha3L3':1,
              'adhoc_exp':1, 'Oalpha2_vp':1, 'Oalpha2_sp':1}
    elif RCs == 'base':
        RC = {'base':1,'Oalpha':0,'Oalpha2L2':0,'Oalpha2L':0,'Oalpha3L3':0,
              'adhoc_exp':0, 'Oalpha2_vp':0, 'Oalpha2_sp':0}
    elif RCs == 'oa':
        RC = {'base':1,'Oalpha':1,'Oalpha2L2':0,'Oalpha2L':0,'Oalpha3L3':0,
              'adhoc_exp':0, 'Oalpha2_vp':0, 'Oalpha2_sp':0}
    else:
        return None

    # parameters must be in a numpy array to pass to a fortran shared-object
    params = numpy.array([mp['rho'],mp['eta'],mp['xi'],mp['delta'],
                          ame2amu2,alme,alpha_qed0,x0,0.0,0.0,0.0])

    # call the fortran functions, but MUST use lower case for func. names
    isotropic,asymmetric = 0.,0.
    if RC['base']:
        isotropic  += fg_nla.ffm0(x,params)
        ## for xi explanation, see https://twist.phys.ualberta.ca/forum/view.php?bn=twist_physics&key=1095896450
        asymmetric += (1/mp['xi'])*fg_nla.ggm0(x,params)
    if RC['Oalpha']:
        isotropic  += fg_nla.ffm1(x,params)
        asymmetric += fg_nla.ggm1(x,params)
    if RC['Oalpha2L2']:
        isotropic  += fg_nla.ffl2(x,params)
        asymmetric += fg_nla.ggl2(x,params)
    if RC['Oalpha2L']:
        isotropic  += fg_nla.fnl2(x,params)
        asymmetric += fg_nla.gnl2(x,params)
    if RC['Oalpha3L3']:
        isotropic  += fg_nla.ffl3(x,params)
        asymmetric += fg_nla.ggl3(x,params)
    if RC['adhoc_exp']:
        isotropic  += fg_nla.fahe(x,params)
        asymmetric += fg_nla.gahe(x,params)
    if RC['Oalpha2_vp']:
        isotropic  += fg_nla.ffvp(x,params)
        asymmetric += fg_nla.ggvp(x,params)
    if RC['Oalpha2_sp']:
        isotropic  += fg_nla.ffsp(x,aldp,params)
        asymmetric += fg_nla.ggsp(x,aldp,params)

    return isotropic + cos_theta*mp['xi']*asymmetric
##    return isotropic - cos_theta*mp['xi']*asymmetric