# -*- coding: utf-8 -*-
"""
Created on Thu Jul 14 18:29:52 2022

@author: frede
"""

import numpy as np
from astropy.table import Table
import matplotlib.pyplot as plt

''' Constants '''
h = 6.62607015e-34 #Planck constant in J s
c = 299792458 #Light speed in m s-1

''' Conversions '''
ev = 1.602176634e-19 # 1 ev = 1.602176634e-19 J and value of electron charge

''' calculus of wavelength '''
wave_intensity = Table.read('./converted/new_habing1968.txt', format = 'ascii')
wavelength = wave_intensity['col1']

density_range = (wavelength >= h/ev * c *1e9/13.6) & (wavelength <= h/ev * c * 1e9/6)


habing_RF = np.trapz( ( (-25/6) * wavelength[density_range]**2/1e6 +\
                         (25/2) * wavelength[density_range]/1e4 -\
                       (13/300) ) * 1e-14, wavelength[density_range])

frac = wavelength/100   
wavelength_intensity = [] 
Draine_list = []
Mathis_list = []   
for i in frac:
    # test =  ( (c*1e2)/(2*np.pi) ) * ( (-25/6) * wavelength[i]**2/1e6 +\
    #                          (25/2) * wavelength[i]/1e4 -\
    #                        (13/300) ) * 1e-14
    # wavelength_intensity.append(test)
    if i >= 0.912:
        Draine = 4e-14 * 1.71 * (1/i**5) * (31.016*i**2 - 49.913*i + 19.897)
        Draine_list.append(Draine)
    
for i in wavelength:
    if i == 91.2 :
        Mathis_list.append(0)
    if i > 91.2 and i <= 110:
        Mathis = 1.287e-9 * (i*1e-3)**(4.4172)
        Mathis_list.append(Mathis)
    if i > 110 and i <= 134:
        Mathis = 6.825e-13 * (i*1e-3)
        Mathis_list.append(Mathis)
    if i > 134 and i <= 246:
        Mathis = 2.373e-14 * (i*1e-3)**(-0.6678)
        Mathis_list.append(Mathis)
    
Habing = [6.6e-14, 7e-14, 4e-14]

plt.figure(figsize=(10.90,8.93))
plt.axis()
plt.semilogy((h/ev *c*1e9)/wavelength, Draine_list,ls='dashdot', c='k',label='Draine (1978)')
plt.semilogy((h/ev *c*1e9)/wavelength[0:len(Mathis_list)], Mathis_list,c='k',label='MMP83')
plt.semilogy([(h/ev *c*1e9)/220, (h/ev *c*1e9)/140, (h/ev *c*1e9)/100], Habing,ls='--',c='k',marker="s",label='Habing (1968)')
plt.ylim([1e-14,2e-13])
plt.xlim([5,14])
plt.xlabel('Energy (eV)')
plt.xticks([6, 8, 10, 12, 13.6])
plt.ylabel('νu$_ν$ (erg cm$^{-3}$)')
plt.text(5.8, 1.6e-14, 'Local Interstellar UV', size=25)
plt.legend()


# pl = plt.plot([1, 4, 3, 7, 8], color = 'blue')
# plt.legend(pl, ['blue'])
# ax2 = plt.gca().twinx()
# pl2 = ax2.plot([80, 70, 30, 40, 10], color = 'red')
# plt.legend(pl2, ['red'])
# plt.show()
# document = Table([wavelength , wavelength_intensity] , names =('col1','col2'))
# document.meta['comments'] = ['wavelength (nm)' , 'intensity (erg s-1 cm-2 nm-1 sr-1)']
# document.write('habing1968.txt', format = 'ascii', overwrite=True)
    
print(habing_RF)