## Take the cross setion data in an hdf5 file and output other ## versions of this cross section required in the database and some ## nice figures # from anh import * import lib_leiden import shutil from scipy import integrate import os import numpy as np from matplotlib import pyplot as plt output_directory_root = f'td/cross_sections/' for (species,input_h5_cross_section) in ( ## no problems 2020-09-09 # ("Al",None), ("AlH",None), ("C",None), ("C2",None), # ("C2H",None), ("C2H2",None), ("C2H4",None), # ("C2H5OH",None), # ("C2H6",None), ("C3",None), ("C3H7OH",None), ("Ca",None), # ("c-C3H",None), ("c-C3H2",None), ("CH",None), ("CH+",None), # ("CH2",None), ("CH2+",None), ("CH3",None), ("CH3CHO",None), # ("CH3CN",None), ("CH4+",None), ("CH3SH",None), ("CN",None), # ("Co",None), ("CO",None), ("CO+",None), ("CO2",None), # ("Cr",None), ("CS",None), ("CS2",None), ("H",None), # ("H-",None), ("H2",None), ("H2+",None), ("H2CO",None), # ("H2O2",None), ("H3+",None), ("HC3H",None), ("HC3N",None), # ("HCl+",None), ("HCN",None), ("HCO",None), ("HCO+",None), # ("HF",None), ("HO2",None), ("l-C3H",None), ("l-C3H2",None), # ("l-C4",None), ("l-C4H",None), ("LiH",None), ("MgH",None), # ("Mn",None), ("N2O",None), ("NaH",None), ("NH",None), # ("NH2",None), ("NH2CHO",None), ("Ni",None), ("NO",None), # ("NO2",None), ("O",None), ("O2",None), ("O2+",None), # ("O3",None), ("OCS",None), ("OH+",None), ("PH",None), # ("PH+",None), ("Rb",None), ("S",None), ("S2",None), # ("SH",None), ("SH+",None), ("SiH",None), ("SiH+",None), # ("SiO",None), ("SO",None), ("SO2",None), ("Ti",None), # ("Zn",None), # ## require hacks 2020-09-09 # ("H2S",None), ("H2O",None), ("CH3NH2",None), ("CH3OH",None), # ("CH4",None), ("Cl",None), ("Fe",None), ("HCl",None), # ("Ca+",None), ("K",None), ("l-C5H",None), ("Li",None), # ("Mg",None), ("N",None), ("N2",None), ("Na",None), # ("NaCl",None), ("NH+",None), ("NH3",None), ("OH",None), # ("P",None), ("Si",None), ): print( species) ## try to get hdf5 cross section file automatically if not given explicitly if input_h5_cross_section is None: for fn in ( f'~/data/species/{species}/cross_sections/all_cross_sections.h5', f'~/data/species/{species}/cross_sections/all_cross_sections.hdf5', ): if os.path.exists(os.path.expanduser(fn)): input_h5_cross_section = fn break else: raise Exception(f'Could not find input_h5_cross_section automatically.') ## load hdf5 into a dict input_h5_cross_section = lib_leiden.expand_path(input_h5_cross_section) data = lib_leiden.hdf5_to_dict(input_h5_cross_section) ## HACKS if 'header' in data and 'README' not in data: data['README'] = data.pop('header') if species in ('H2S','H2O','CH3NH2','CH3OH','Mg','N2'): for key in ('photoabsorption','photodissociation','photoionisation'): if key in data: data[key][data[key]<0] = 0 if species in ('H2O','Cl','Fe','HCl','Mg','N','NH3','Ca',): for key in list(data): if key not in ('README','wavelength','photoabsorption','photodissociation','photoionisation'): data.pop(key) if species in ('Na','P','Si'): data['photoabsorption'] = data.pop('photabsorption') if species in ('Ca','Cl','Fe','K','Mg','N','Na','P','Si',) and 'photoabsorption' not in data: data['photoabsorption'] = data['photoionisation'] ## some basic error checking assert 'README' in data,'README not in data' assert 'wavelength' in data,'wavelength not in data' assert 'photoabsorption' in data,'photoabsorption not in data' assert np.all(np.diff(data['wavelength'])>0),'Wavelength scale not monotonically increasing.' assert np.all(data['photoabsorption']>=0),'Photoabsorption cross section less than zero: '+str(np.min(data['photoabsorption'])) assert 'photodissociation' not in data or np.all(data['photodissociation']>=0),'Photodissociation cross section less than zero.: '+str(np.min(data['photodissociation'])) assert 'photoionisation' not in data or np.all(data['photoionisation']>=0),'Photoionisation cross section less than zero.: '+str(np.min(data['photoionisation'])) for key in data.keys(): assert key in ('README','wavelength','photoabsorption','photodissociation','photoionisation'),f'Unknown data key: {key}' data['README'] = str(data['README']) # avoids some kind of encoding problem products = [key for key in ('photoabsorption','photodissociation','photoionisation') if key in data] ## output h5 file output_directory = f'{output_directory_root}/{species}/' lib_leiden.mkdir_if_necessary(output_directory) lib_leiden.dict_to_hdf5(f'{output_directory}/{species}.h5',data) ## print integrated cross sections print('integrated cross section:') print(format('products','>10'),' '.join([format(key,'>18') for key in products])) print(format('h5','>10'),' '.join([format(integrate.trapz(data[key],data['wavelength']),'18.5e') for key in products])) ## output cross section into a multicolumn text file np.savetxt( f'{output_directory}/{species}.txt', np.column_stack([data[key] for key in ['wavelength'] + products]), fmt=['%12.7f']+['%17.5e' for t in products], header=data['README']+'\n\n'+' '.join([format('wavelength','>10s')] + [format(t,'>17s') for t in products]), ) ## output cross section into a multicolumn text file on a resampled wavelength grid dwl = 0.1 # the sampling spacing wl = np.arange( data['wavelength'][-1]-data['wavelength'][-1]%0.1, data['wavelength'][0]+1e-8, -dwl)[::-1] # do it backwards and reverse so threshold is not missed by arange d = {'wavelength':wl} for key in products: d[key] = lib_leiden.resample(data['wavelength'],data[key],wl) np.savetxt( f'{output_directory}/{species}_0.1nm.txt', np.column_stack([d[key] for key in ['wavelength'] + products]), fmt=['%12.7f']+['%17.5e' for t in products], header=data['README']+'\n\n'+' '.join([format('wavelength','>10s')] + [format(t,'>17s') for t in products]), ) print(format('0.1nm','>10'),' '.join([format(integrate.trapz(d[key],d['wavelength']),'18.5e') for key in products])) ## leiden database files, output pd and pi cross sections to text files ## ## RETHINK THIS -- DOES NOT SEPARATE LINES FROM CONTINUUM OR LIMIT TO 91.2nm print("warning: pd/pi files not very well implemented -- no lines") tspecies = lib_leiden.translate_species(species,'standard','leiden') if 'photodissociation'in data: lib_leiden.save_cross_section_leiden_photodissoc_database( filename=f'{output_directory}/{tspecies}.pd', header=f"Photodissociation cross section of {tspecies}. "+data['README'].replace('\n',' '), lines_wavelength=[], lines_integrated_cross_section=[], continuum_wavelength=data['wavelength'], continuum_cross_section=data['photodissociation'],) if 'photoionisation' in data: tspecies = lib_leiden.translate_species(species,'standard','leiden') lib_leiden.save_cross_section_leiden_photodissoc_database( filename=f'{output_directory}/{tspecies}.pi', header=f"Photoionisation cross section of {tspecies}. "+data['README'].replace('\n',' '), lines_wavelength=[], lines_integrated_cross_section=[], continuum_wavelength=data['wavelength'], continuum_cross_section=data['photoionisation'],) ## make a figure lib_leiden.presetRcParams('article_single_column') fig = plt.figure(1) fig.clf() ax = fig.gca() rasterized = len(d['wavelength'])>2000 # if too much data then rasterize for iproduct,product in enumerate(products): ax.plot(data['wavelength'],data[product],color=lib_leiden.newcolor(iproduct),label=product,rasterized=rasterized) lib_leiden.legend(loc='upper right') ax.set_xlabel(r'Wavelength (nm)') ax.set_ylabel(r'Cross section (cm$^2$)') ax.annotate(lib_leiden.translate_species(species,'standard','matplotlib'), (0.03,0.93),ha='left',va='top',xycoords='axes fraction') if rasterized: fig.savefig(f'{output_directory}/{species}.pdf',dpi=350) else: fig.savefig(f'{output_directory}/{species}.pdf') fig.savefig(f'{output_directory}/{species}.png',dpi=350) ## show last figure generated plt.show()