"""
Some double beta decal rates
"""


from mass_table import MassTable
import numpy as np
from physconst import SEC, ME, MEV, CLIGHT

class GSDecayRate(CDatRec):
    def __init__(
            self, reaction, rate,
            Qnu = None,
            id = 0,
            loss = True,
            flags = 0,
            ):
        formula = 7
        reverse = False
        decay = True
        reaction =  Reaction(reaction)
        if Qnu is None:
            global mass
            try:
                mass
            except:
                mass = MassTable()
            Q = 0.
            nnu = 0
            for i,n in reaction.nuc_in.items():
                Q += mass[i]['mass_excess'] * n
            for i,n in reaction.nuc_out.items():
                Q -= mass[i]['mass_excess']
            for i,n in reaction.lep_in.items():
                if i.Z > 0:
                    assert i.flavor == 1
                    Q += ME * CLIGHT**2 / MEV * n * 2
                assert i.Z != 0
            for i,n in reaction.lep_out.items():
                if i.Z > 0:
                    assert i.flavor == 1
                    Q -= ME * CLIGHT**2 / MEV * n * 2
                else:
                    nnu += n
            if nnu == 0:
                Qnu = 0.
            elif nnu == 1:
                Qnu = Q / 3.
            elif nnu == 2:
                Qnu = Q * 2 / 3
            else:
                raise Exception('unknown nu mode')

        self.reaction = reaction
        self.reverse = reverse
        self.decay = decay
        self.formula = formula
        self.id = id
        self.loss = loss
        self.fcoef = list([rate, Qnu])
        self.icoef = list()
        self.flags = flags # decay and reverse could be flags
        self.label = '(Wikipedia)'

    def __str__(self):
        s = f'{self.formula}: {self.reaction}'
        if self.reverse:
            s += ', +rev'
        if self.loss == True:
            s += f', L={self.fcoef[-1]:7.3f}'
        if self.id != 0:
            s += f', id={self.id}'
        if self.flags != 0:
            s += f', flags={self.flags}'
        return f'{{{s}}}'


class DoubleBeta(CDat):
    """
    from https://en.wikipedia.org/wiki/Double_beta_decay (20210612)
    """
    def __init__(self):
        print(f'[{self.__class__.__name__}] should be mix-in instead')
        self.add_doublebeta()

    def add_doublebeta(self):
        if not hasattr(self, 'rates'):
            self.rates = list()
        rates = [
            GSDecayRate(' ca48(b-b-)', np.log(2) / (0.064e21 * SEC)),
            GSDecayRate(' ge76(b-b-)', np.log(2) / (1.926e21 * SEC)),
            GSDecayRate(' kr78(e-e-)', np.log(2) / (9.2e21 * SEC)),
            GSDecayRate(' se82(b-b-)', np.log(2) / (0.096e21 * SEC)),
            GSDecayRate(' zr96(b-b-)', np.log(2) / (0.0235e21 * SEC)),
            GSDecayRate('mo100(b-b-)', np.log(2) / (0.69e21 * SEC)),
            GSDecayRate('pd110(b-b-)', np.log(2) / (1.5e21 * SEC)),
            GSDecayRate('cd116(b-b-)', np.log(2) / (0.028e21 * SEC)),
            GSDecayRate('te128(b-b-)', np.log(2) / (2200e21 * SEC)),
            GSDecayRate('te130(b-b-)', np.log(2) / (0.82e21 * SEC)),
            GSDecayRate('xe124(e-e-)', np.log(2) / (18e21 * SEC)),
            GSDecayRate('xe136(b-b-)', np.log(2) / (2.165e21 * SEC)),
            GSDecayRate('ba130(e-e-)', np.log(2) / (1.6e21 * SEC)),
            GSDecayRate('nd150(b-b-)', np.log(2) / (0.00911e21 * SEC)),
            GSDecayRate(' u238(b-b-)', np.log(2) / (2.0e21 * SEC)),
        ]
        self.rates.extend(rates)