#!/usr/local/bin/python3 """ Command-line program to calculate alkali opacity from an expansion in density Inputs: temperature - gas temperature should match the file [kelvin] neutral_density - neutral particle density [cm-3] core_delta_nu - core wavenumber grid spacing [cm-1] n_core - number of points in the line core grid above and below center expansion file - name of expansion table file in the current working directory or with full path, with specific format for this species opacity file - base name of output opacity files Outputs: File with opacity versus wavelength [A] File with opacity versus delta wavenumber [cm^-1] File with Lorentzian core versus wavelength [A} File with Lorentzian core versus delta wavenumber [cm^-1] Contains: read_opacity_file(infile) - function to read and parse an alkali density expansion file Cautions: The tables are similar to but not formatted in the same way as those used in TLusty. A Lorentzian even with the core parameters from the table does not represent the core accurately due to an intrinsic core asymmetry that may differ from line to line. When possible, for the core use the frequencies in the table and interpolate. Do not use another Lorentzian calculated with ad hoc van der Waals coefficients. """ import os import sys import numpy as np # Set this True for information while the file is processed verbose_flag = True # Set this True for addtional diagnostics in case of failure diagnostic_flag = False # Read a file of opacity data return the contents def read_opacity_file(infile): # Read the expansion data from a formatted text table # Input: # Name (with path if needed) of the file # Return: # delta_f - numpy array of frequency differences [cm^-1] from the reference line center # expansion_data - numpy array of expansion coefficients for all delta_f # reference_wavelength - reference line center wavelength [A} from the table # reference_density - reference density at which the table was computed # volume_normalization - density normalization factor referred to as vn in Fortran code # opacity_normalization_pir0f - normalization factor from frequency-normalized profile # scales to physical opacity as cross-section [cm^2] # Handle the table file table_fp = open(infile, 'r') table_text = table_fp.readlines() n_lines = len(table_text) table_fp.close() # Parse the text header i_line = -1 # zero-based index to last processed from the table file reference_temperature = 273. # default value should be read from the text header later for line in table_text: i_line = i_line + 1 entry = line.strip().split() if diagnostic_flag: print("Reading line ", i_line, "with entry", entry[0]) if (entry[0] == "TK:"): # Test for the temperature so it can be passed on to applications reference_temperature = float(entry[1]) if verbose_flag: print("Temperature [K] from the text header is ", reference_temperature) continue if (entry[0] == "end"): # Test for an end to the text header and exit this text reader break if (i_line == n_lines): # Something went wrong print("The file ", infile, "is truncated without a text header. Check for an end to the text header.\n") exit() # The text header has been processed if i_line + 3 > n_lines: # Test that there are at least 3 more lines to process print("The file ", infile, "is truncated without a data header. Check for three data header lines.\n") exit() # Parse the three data header lines # Read the first data header line i_line = i_line + 1 line = table_text[i_line] entry = line.strip().split() reference_wavelength = float(entry[0]) reference_density = float(entry[1]) volume_normalization = float(entry[2]) opacity_normalization_pir0f = float(entry[3]) # Read the second data header line i_line = i_line + 1 line = table_text[i_line] entry = line.strip().split() n_blocks = int(entry[0]) n_coefficients = int(entry[1]) # Read the third data header line i_line = i_line + 1 line = table_text[i_line] entry = line.strip().split() core_width = float(entry[0]) core_shift = float(entry[1]) # The data header has been processed if verbose_flag: print("Reference wavelength [A]: ", reference_wavelength) print("Reference density [cm^-3]: ", reference_density) print("Volume normalization factor: ", volume_normalization) print("Opacity normalization factor: ", opacity_normalization_pir0f) # Create arrays to hold the table data delta_f = np.zeros(n_blocks) expansion_data = np.zeros((n_blocks,n_coefficients)) # Read the frequency and the coefficients one block at a time # Count blocks and coefficients as we go for i_block in range(n_blocks): # Within each block read delta_f first, then the coefficients # The length of lines within the data blocks may be inconsistent block to block # Each block contains n_coefficients values # Count them to know when to advance to the next block # For the K-H_2 tables there are exactly 5 lines per block # Other tables may differ i_coefficient = 0 # zero-based coefficient counter within this block while i_coefficient < n_coefficients: # Process a line from the block and split it into entries i_line = i_line + 1 line = table_text[i_line] entry = line.strip().split() n_entries = len(entry) # Parse the entries on this line into delta_f and coefficients # Each block is saved in delta_f[i] and expansion_data[i,j] # i -> block and j -> coefficient for i_entry in range(n_entries): if i_coefficient == 0 and i_entry == 0: # The first entry of a block is delta_f delta_f[i_block] = float(entry[i_entry]) if verbose_flag: print("Processing delta frequency ", delta_f[i_block]) else: # Other entries in this block are coefficients expansion_data[i_block, i_coefficient] = float(entry[i_entry]) i_coefficient = i_coefficient + 1 if verbose_flag: print("The expansion data in table file ", expansion_file, " have been read.\n") table_tuple = (delta_f, expansion_data, reference_wavelength, reference_density, reference_temperature, volume_normalization, opacity_normalization_pir0f, core_width, core_shift) return table_tuple # Main program begins here # Read and parse command line arguments if len(sys.argv) == 7: opacity_temperature = float(sys.argv[1]) opacity_density = float(sys.argv[2]) core_delta_nu = float(sys.argv[3]) n_core = int(sys.argv[4]) expansion_file = sys.argv[5] opacity_file_basename = sys.argv[6] else: print(" ") print("Usage: alkali_opacity_from_expansion.py temperature, neutral_density, core_delta_nu, n_core, expansion.txt, opacity_basename ") print(" ") sys.exit("Use the density expansion to find the opacity\n") # Read the expansion table (delta_f, expansion_data, reference_wavelength, reference_density, reference_temperature, volume_normalization, opacity_normalization_pir0f, core_width, core_shift) = read_opacity_file(expansion_file) (n_frequencies, n_coefficients) = np.shape(expansion_data) if verbose_flag: print("Reference wavelength [A}: ", reference_wavelength) print("Table contains data for ", n_frequencies, " explicit frequencies") print("Expansion table for temperature [K]: ", reference_temperature) print("Calculated opacity at temperature [K]: ", opacity_temperature) print("Exapansion table is for neutral atom density [cm^-3]: ", reference_density) print("Calculated opacity for neutral atom density [cm^-3]: ", opacity_density) print("\n") # Go through the table line by line and find the opacity for each entry # This routine uses the frequency offset assigned to that entry and does not interpolate between entries # Create arrays to hold the opacity and the intermediate density parameters # Return opacity for offset frequency [cm^-1] and vacuum wavelength [A] with explicit values calculated from the table # and other values inserted for line core on a mesh requested by the command line opacity_wn_file = opacity_file_basename + "_wn.dat" opacity_wl_file = opacity_file_basename + "_wl.dat" lorentzian_wn_file = opacity_file_basename + "_lorentzian_wn.dat" lorentzian_wl_file = opacity_file_basename + "_lorentzian_wl.dat" # Calculate the opacity from the table for each relative frequency and the wavelengths for those frequencies density_data = np.zeros(n_coefficients) n_data = n_frequencies opacity = np.zeros(n_data) wavenumber = np.zeros(n_data) wavelength = np.zeros(n_data) # The following Python code is adapted from the Fortran program lect_sig.f with these lines # vn1=dens_out/dens_in # vns=vn1*vn # xnorm=dexp(-vns) # cdens(1)=vn1 # # do i = 2,nexp # cdens(i)=cdens(i-1)*cdens(1) # enddo # do j=1,ntable # pcour=0 # do i=1,nexp # pp(i)= cdens(i)*plall(j,i) # pcour= pp(i)+pcour # enddo # pptot(j)=pcour # sig=pptot(j)*xnorm*pir0f # xl(j)=1.0/(1.0e-8*om(j)+1.0/xlam_c) # if(sig.gt.0.) then # write(11,'(2e15.7)') om(j),sig # write(10,'(2e15.7)') xl(j),sig # endif # enddo # Set up the density array once - it is the same for all table frequencies normalized_density = opacity_density/reference_density cross_section_normalization = np.exp(-volume_normalization*normalized_density) density_data = np.zeros(n_coefficients) for j in range(n_coefficients): # Successive powers of the normalized density for each coefficient in the table density_data[j] = normalized_density**j # Compute the opacity at each frequency # Uses numpy arrays opacity, wavenumber, and wavelength # Scalars reference_wavelength from the table, and reference_wavenumber computed for i in range(n_data): opacity[i] = np.sum(density_data[:]*expansion_data[i,:]) wavenumber[i] = delta_f[i] opacity = opacity_normalization_pir0f*cross_section_normalization*opacity # Assign a wavelength for each frequency reference_wavenumber = 1.e8 / reference_wavelength wavelength = 1.e8 / (reference_wavenumber + wavenumber) # Save the table spectra opacity_spectrum_wl = np.column_stack((wavelength,opacity)) opacity_spectrum_wn = np.column_stack((wavenumber,opacity)) np.savetxt(opacity_wl_file, opacity_spectrum_wl) np.savetxt(opacity_wn_file, opacity_spectrum_wn) # Calculate the analytical symmetrical Lorentzian core from the core parameters in the table # Half width at half maximum for this density and temperature [cm^-1]: hwhm # Shift of line center for this density and temperature [cm^-1]: shift width = core_width*normalized_density shift = core_shift*normalized_density if verbose_flag: print("Calculating a Lorentzian for the core.") print("Table Temperature [K]: ", reference_temperature) print("Requested Temperature [K]: ", opacity_temperature) print("Density [cm^-1]: ", opacity_density) print("Half width at half maximum [cm^-1]: ", width) print("Shift [cm^-1]: ", shift) print("\n") # Set up arrays to hold 2*n_core+1 values # lect_sig.f uses 2*n_core - 1 core_lorentzian = np.zeros(2*n_core + 1) core_wavelength = np.zeros(2*n_core + 1) core_wavenumber = np.zeros(2*n_core + 1) for n in range(2*n_core + 1): core_wavenumber[n] = (n - n_core)*core_delta_nu core_wavelength = 1.e8 / (reference_wavenumber + core_wavenumber) # Analytical symmetric Lorentzian with unit area integrated over wavenumber # lect_sig.f provides for an asymmetric core when those parameters are known core_lorentzian = (width/np.pi) / (width*width + (core_wavenumber - shift)* (core_wavenumber - shift)) # Normalize the Lorentzian to opacity for this system core_lorentzian = opacity_normalization_pir0f*core_lorentzian # Save the Lorentzian core lorentzian_wl = np.column_stack((core_wavelength,core_lorentzian)) lorentzian_wn = np.column_stack((core_wavenumber,core_lorentzian)) np.savetxt(lorentzian_wl_file, lorentzian_wl) np.savetxt(lorentzian_wn_file, lorentzian_wn) print("The %i opacities and a core Lorentizian for these conditions have been saved in: \n %s and %s" % (n_data, opacity_wn_file, opacity_wl_file)) print("\n") exit()