import numpy as np
from becquerel.tools.isotope import Isotope
from becquerel.tools.isotope_qty import IsotopeQuantity, NeutronIrradiation
import datetime
from becquerel import Spectrum
import efficiencies as ef
from bs4 import BeautifulSoup
import urllib.request
import math
import NAA_Isotopes as na

"""
Explanation: Isotope_composition_calculator.py Tutorial

The purpose of this notebook is to demonstrate how Isotope_composition_calculator.py works including what the required inputs are needed and what it outputs
First, import the required documents:
End of explanation
"""


spec_S1 = Spectrum.from_file('/Users/jackiegasca/Documents/spectras/Sample4_24h_C.Spe')

"""
Explanation: The python file efficiencies is a script in charge of generating the efficiency spectra needed to calculate the efficiency values for the variety of energy peaks used for analysis. The efficiency spectra was created using a variety of sources whose imported spectra can be seen above
NAA_Isotopes is a script with a class class containing dictionaries of data for every isotope that is observed in neutron activation analysis (of marine life). This is used to identify the isotope that belongs to an observed energy peak
Load a spectra for observation:
(The spectra are uploaded onto the code itself, since an interface has yet to be created)
End of explanation
"""


back_spec = Spectrum.from_file('/Users/jackiegasca/Documents/2017.5.1_long_background.Spe')

"""
Explanation: Load a background spectra as well
End of explanation
"""


spec_S1_ener_spec = spec_S1.energies_kev[0:len(spec_S1)]
back_ener_spec = back_spec.energies_kev[0:len(back_spec)]

"""
Explanation: The code then converts bins into energy so that peaks can be identified:
End of explanation
"""


isotope_list = [na.Na_24, na.Br_82, na.Cs_134, na.K_40, na.Se_75, na.Sr_85,
                na.Ca_47, na.Rb_86, na.Zn_65, na.K_42, na.Fe_59, na.Au_198,
                na.Co_60, na.As_76, na.Hg_203, na.Cr_51, na.Sb_122, na.Sb_124]

"""
Explanation: For now, since the code is still under development, it needs to be given an isotope_list for it to use. This is basically telling the code what to look for in the given spectra. Since I previously identified some of the observed peaks, I provided the list of isotopes for which the peaks belong to:
End of explanation
"""


irr_start = '2017-04-27 14:02:00'
irr_stop = '2017-04-27 14:12:00'
flux = 3.1e11
N_0 = 6.02e23

"""
Explanation: For now, the code needs to be given the logistics of the irradiation, such flux and start time. Eventually, the code will be given a file it can parse to get this information. For now, this method must be used
End of explanation
"""


def IsotopeActivity():
    iso_name = []
    iso_energy = []
    iso_cps = []
    iso_br = []
    for iso in isotope_list:

        E = iso.energies['energy'][0]
        FWHM = ((2.355 * (0.09 * 0.00296 * E) ** 0.5) ** 2
                + (1.3) ** 2) ** 0.5 #keV
        start = E - 1 * FWHM
        end = E + 1 * FWHM
        bkgd_start = E - 2 * FWHM
        bkgd_end = E + 2 * FWHM

        en = (np.abs(spec_S1_ener_spec - E)).argmin()
        val1 = (np.abs(spec_S1_ener_spec - start)).argmin()
        val2 = (np.abs(spec_S1_ener_spec - end)).argmin()
        val3 = (np.abs(spec_S1_ener_spec - bkgd_start)).argmin()
        val4 = (np.abs(spec_S1_ener_spec - bkgd_end)).argmin()

        cps_values = spec_S1.cps_vals[val1:val2]
        max_cps_index = np.argmax(cps_values)
        ex_val = max_cps_index - (en - val1)
        val1 = val1 + ex_val
        val2 = val2 + ex_val
        val3 = val3 + ex_val
        val4 = val4 + ex_val

        peak_vals = spec_S1.cps_vals[val1:val2]
        back_vals = back_spec.cps_vals[val1:val2]
        peak_vals_sub = [a - b for a, b in zip(peak_vals, back_vals)]

        bkgd_vals1 = spec_S1.cps_vals[val3:val1 - 1]
        back_bkgd_vals1 = back_spec.cps_vals[val3:val1 - 1]

        bkgd_vals2 = spec_S1.cps_vals[val2 + 1:val4]
        back_bkgd_vals2 = back_spec.cps_vals[val2+1:val4]

        back_sub1 = [a - b for a, b in zip(bkgd_vals1, back_bkgd_vals1)]
        back_sub2 = [a - b for a, b in zip(bkgd_vals2, back_bkgd_vals2)]

        bkgd_cps = (sum(back_sub1) + sum(back_sub2)) / (len(back_sub1)
                    + len(back_sub2))
        peak_vals[:] = [x - bkgd_cps for x in peak_vals_sub]
        peak_cps = sum(peak_vals)

        name = '{0}_{1}'.format(iso.Symbol, iso.A + 1)
        iso_name.append(name)
        iso_energy.append(E)
        iso_cps.append(peak_cps)
        iso_br.append(iso.energies['branching_ratio'][0])
    #return(iso_name, iso_energy, iso_cps, iso_br)
    isotope_activities = []
    for j in range(len(iso_name)):
        ef_en = (np.abs(ef.x - iso_energy[j])).argmin()
        activ = iso_cps[j] / (ef.fit[ef_en] * iso_br[j])
        isotope_activities.append(activ)
    return(iso_name, iso_energy, iso_cps, iso_br, isotope_activities)
    for n in range(len(isotope_activities)):
        #statement = 'The activity of ' + '{0}'.format(iso_name[n]) + ' at ' + '{0}'.format(spec_S1.start_time) + ' is ' + '{0}'.format(iso_cps[n])
        statement = 'The activity of {0} at {1} is {2} bq'.format(iso_name[n], spec_S1.start_time, iso_cps[n])
        print(statement)


"""
Explanation: IsotopeActivity() was created to search for the activity (in counts per second) in the energy peaks it was given from the list of isotopes. This function returns lists containing the isotope names and respective energies, cps, and branching ratios:
End of explanation
"""


lists = IsotopeActivity()
iso_cps = lists[2]
iso_name = lists[0]
iso_energy = lists[1]
iso_br = lists[3]
isotope_activities = lists[4]


"""
Explanation: The function then needs to be called to generate the data needed to proceed with the concentrtion calculations. The values are assigned to their corresponding list as follows:
End of explanation
"""


def Remover():
    for i in iso_cps:
        if i <= 0:
            n = iso_cps.index(i)
            iso_name.remove(iso_name[n])
            iso_energy.remove(iso_energy[n])
            iso_br.remove(iso_br[n])
            isotope_activities.remove(isotope_activities[n])
            iso_cps.remove(iso_cps[n])
        else:
            pass
    return(iso_name, iso_cps, iso_energy, iso_br, isotope_activities)

lists = Remover()
lists = Remover()
lists = Remover()


"""
Explanation: The function Remover() was created to remove any negative activities, since this implies the isotope isn't really there. It is then executed a couple of times to assure all the negative activities have been removed
End of explanation
"""


def Concentration():
    conc = []
    for i in range(len(iso_name)):
        c = iso_name[i].split('_')
        abb = c[0]
        A = int(c[1])
        A_0 = A - 1
        iso_2 = '{0}-{1}'.format(abb, A_0)
        iso_1 = '{0}-{1}'.format(abb, A)
        nuclide = Isotope(iso_1)
        def urlcreator(abb, A_0):
            A_num = str(A_0)
            if len(A_num) == 1:
                A_num = '00' + A_num
            elif len(A_num) == 2:
                A_num = '0' + A_num
            else:
                A_num = A_num
            url = 'http://wwwndc.jaea.go.jp/cgi-bin/Tab80WWW.cgi?/data' \
                    + '/JENDL/JENDL-4-prc/intern/' + abb + A_num + '.intern'
            html = urllib.request.urlopen(url)
            bslink = BeautifulSoup(html, 'lxml')

            return(bslink)

        bslink = urlcreator(abb, A_0)
        def tabledata(bslink):
            '''extracts data from the jaea website'''

            table = bslink.table
            table_rows = table.find_all('tr')
            for tr in table_rows:
                td = tr.find_all('td')
                row = [i.text for i in td]

                if len(row) == 7:
                    if row[0] == 'total       ':
                        x_sec = row[1]
                        x_sec_s = x_sec.split(' ')
                        x_val = float(x_sec_s[0])
                        barn = x_sec_s[1]
                        if barn[1] == 'k':
                            x_val = 10**(3) * x_val
                            return(x_val)
                        elif barn[1] == 'm':
                            x_val = 10**(-3) * x_val
                            return(x_val)

                        elif barn[1] == '&':
                            x_val = 10**(-6) * x_val
                            return(x_val)
                        else:
                            x_val = x_val
                            return(x_val)

                    else:
                        pass

                else:
                    pass
            return(x_val)

        x_val = tabledata(bslink)
        #print(x_val)

        quantity = IsotopeQuantity(nuclide, date=spec_S1.start_time, bq=isotope_activities[i])
        irrad_quan = quantity.bq_at(irr_stop)
        irrad_act = IsotopeQuantity(nuclide, date=irr_stop, bq=irrad_quan)
        ni = NeutronIrradiation(irr_start, irr_stop, n_cm2_s=flux)
        init_comp = ni.activate(x_val, initial=Isotope(iso_2), activated=irrad_act)
        init_comp.is_stable = True
        print(init_comp)
        

"""
Explanation: The function Concentration() calcualtes the composition of the isotope in the iso_name list using it's activity and branching ratio at the energy peak. The first part of the function converts the isotopes name into a format that is compatible with the becquerel code. The function urlcreator(abb, A_0) takes the name of the element and the preirradiated nuclide number to search through the jaea library to find the cross section of the preirradiated isotope. It returns a link with the cross section information. After the link is accessed by the webscraper, tabledata(bslink) is a function that accesses the link and converts the data from the web into a table that can be parsed. It looks for the thermal neutron cross section and returns that number. The if statement towards the bottom converts the 'k', 'm' and '&' values into their respective scientific notation. The function returns the approproate cross section. 
When the cross section is returned, Concentration() uses the tools from becquerel to calculate the concentration, which is returned at the end of the function.
End of explanation
"""


u = Concentration()
print(u)

"""
Explanation: Finally, the function is called to generate the results, which are then printed
End of explanation
"""


