## Take the cross setion data in an hdf5 file and output a text file on the same
## wavelength grid.  Output another two more text files in the .pd and
## .pi formats originally used in the Leiden database.
from anh import *
import spectra

# input_filename = '../all_cross_sections/hdf5/HF.hdf5'
input_filename = '../cross_sections/O2/O2.hdf5'
output_directory = 'td/'

species = my.rootname(input_filename)
species_matplotlib=r'O$_2$'
print( species)

## load hdf5 into a dict 
d = my.hdf5_to_dict(input_filename)

## output cross section into a multicolumn text file
README = d.pop('README')
keys,data = [],[]
for key in ('wavelength', 'photoabsorption', 'photodissociation', 'photoionisation'):
    if key not in d: continue
    keys.append(key)
    data.append(d[key])
my.array_to_file(
    output_directory+'/'+my.rootname(input_filename)+'.txt', # output filename
    *data,                               # data arrays in columns
    fmt=['%12.7f']+['%17.5e' for t in range(len(data)-1)],   # format of columns
    header=README+'\n\n'+' '.join(
        [format(keys[0],'>10s')] + 
        [format(t,'>17s') for t in keys[1:]])
)

## output cross section into a multicolumn text file on a resampled wavelength grid
dwl = 0.1                       # the sampling spacing
wl = np.concatenate((np.arange(d['wavelength'][0],d['wavelength'][-1],dwl), [d['wavelength'][-1]]))
keys,data = ['wavelength'],[wl]
for key in ('photoabsorption', 'photodissociation', 'photoionisation'):
    if key not in d: continue
    keys.append(key)
    xs = lib_molecules.resample_data(d['wavelength'],d[key],wl)
    data.append(xs)
    ## check integrated cross section not changed
    print(input_filename,'old',integrate.trapz(d[key],d['wavelength']))
    print(input_filename,'new',integrate.trapz(xs,wl))
    
my.array_to_file(
    output_directory+'/'+my.rootname(input_filename)+'_0.1nm.txt', # output filename
    *data,                               # data arrays in columns
    fmt=['%12.7f']+['%17.5e' for t in range(len(data)-1)],   # format of columns
    header=README+'\n\n'+' '.join(
        [format(keys[0],'>10s')] + 
        [format(t,'>17s') for t in keys[1:]])
)

## output pd cross section text file
if 'photodissociation' in d:
    lib_molecules.save_cross_section_leiden_photodissoc_database(
        filename=output_directory+'/'+my.rootname(input_filename)+'.pd',
        header="Photodissociation cross section of "+my.rootname(input_filename)+". See the Leiden photo website, http://home.strw.leidenuniv.nl/~ewine/photo/ "+my.date_string(),
        lines_wavelength=[],
        lines_integrated_cross_section=[],
        continuum_wavelength=d['wavelength'],
        continuum_cross_section=d['photodissociation'],
        )

## output pi cross section text file
if 'photoionisation' in d:
    lib_molecules.save_cross_section_leiden_photodissoc_database(
        filename=output_directory+'/'+my.rootname(input_filename)+'.pi',
        header="Photodissociation cross section of "+my.rootname(input_filename)+". See the Leiden photo website, http://home.strw.leidenuniv.nl/~ewine/photo/ "+my.date_string(),
        lines_wavelength=[],
        lines_integrated_cross_section=[],
        continuum_wavelength=d['wavelength'],
        continuum_cross_section=d['photoionisation'],
        )

## make a figure
my.presetRcParams('article_single_column')
fig = plt.figure(1)
ax = fig.gca()
## plot various cross sections
for iproduct,product in enumerate(('photoabsorption','photodissociation','photoionisation',)):
    if product not in d.keys(): continue
    ax.plot(d['wavelength'],d[product],color=my.newcolor(iproduct),label=product)
## finish and save figures
my.legend_colored_text()
ax.set_xlabel(r'Wavelength (nm)')
ax.set_ylabel(r'Cross section (cm$^2$)')
ax.annotate(species_matplotlib,(0.03,0.93),ha='left',va='top',xycoords='axes fraction')
# my.savefig(output_directory+'/cross_sections_'+species+'.pdf')
fig.savefig(output_directory+'/cross_sections_'+species+'.png',dpi=350,fig=fig)