Commit d7f9bd25 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Add a python script to filter optical functions

parent b93429ed
Loading
Loading
Loading
Loading
+403 −0
Original line number Diff line number Diff line
#!/usr/bin/env python3

#   Copyright (C) 2026   INAF - Osservatorio Astronomico di Cagliari
#
#   This program is free software: you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#   
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#   
#   A copy of the GNU General Public License is distributed along with
#   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.

## @package filter_constants
#  \brief Script to create optical function grids for NPTM_code.
#
#  One of the factors that critically affect the duration of a calculation
#  with NPTM_code is the number of wavelengths for which the calculation
#  is executed. This number depends on the definition of the optical functions
#  of the involved materials. The filter_constants.py script can be used to
#  extract optimized grids of optical constants.
#
#  The script requires python3.

import math
import pdb
from sys import argv

## \cond
__version__ = "0.10.10"
allow_plots = True
## \endcond

try:
    import matplotlib.pyplot as plt
except ImportError:
    print("WARNING: MATPLOTLIB not found - disabling plots.")
    allow_plots = False

## \brief Main execution code.
#
# `main()` is the function that handles the creation of the code configuration.
# It returns an integer value as exit code, using 0 to signal successful execution.
#
# \returns result: `int` Exit code (0 = SUCCESS).
def main():
    result = 0
    config = parse_arguments()
    if config['version_mode']:
        print("filter_constants.py v%s."%__version__)
    elif (config['help_mode']):
        print_help()
    else:
        if (config['input_files'] == ''):
            print("ERROR: no input files provided (--in INPUT)!")
            result = 1
        else:
            split_input = config['input_files'].split(',')
            if (len(split_input) == 1):
                result = scan_single_file(config)
            else:
                result = scan_multiple_files(config)
    return result

def get_max_distance(wl_orig, reps_orig, ieps_orig, wl_filtered, reps_filtered, ieps_filtered):
    max_difference = 0.0
    max_value = 0.0
    max_wavelength = 0.0
    differing_set = ""
    filtered_index = 1
    for i in range(len(wl_orig)):
        wl = wl_orig[i]
        if (wl <= wl_filtered[0]):
            continue
        if (wl < wl_filtered[filtered_index]):
            dx = wl - wl_filtered[filtered_index - 1]
            dry = reps_filtered[filtered_index] - reps_filtered[filtered_index - 1]
            diy = ieps_filtered[filtered_index] - ieps_filtered[filtered_index - 1]
            dwl = wl_filtered[filtered_index] - wl_filtered[filtered_index - 1]
            rp = reps_filtered[filtered_index - 1] + dry * dx / dwl
            ip = ieps_filtered[filtered_index - 1] + diy * dx / dwl
            rdiff = rp - reps_orig[i]
            idiff = ip - ieps_orig[i]
            if (rdiff < 0.0):
                rdiff *= -1.0
            if (idiff < 0.0):
                idiff *= -1.0
            if (rdiff > max_difference):
                max_difference = rdiff
                max_value = reps_orig[i]
                max_wavelength = wl
                differing_set = "real"
            if (idiff > max_difference):
                max_difference = idiff
                max_value = ieps_orig[i]
                max_wavelength = wl
                differing_set = "imaginary"
        else:
            filtered_index += 1
            if (filtered_index == len(wl_filtered)):
                break # for i
        # end of wl step check
    result = {
        'max_wavelength': max_wavelength,
        'max_difference': max_difference,
        'max_value': max_value,
        'differing_set': differing_set
    }
    return result

## \brief Parse the command line arguments.
#
#  The script behaviour can be modified through a set of optional arguments.
#  The purpose of this function is to parse the command line in search for
#  such arguments and prepare the execution accordingly.
#
#  \returns config: `dict` A dictionary containing the script configuration.
def parse_arguments():
    config = {
        'input_files': '',
        'output_files': '',
        'help_mode': False,
        'version_mode': False,
        'force_peaks': True,
        'make_plots': allow_plots,
        'step': 5.0e-8,
        'threshold': 0.1,
        'wl_start': 0.0,
        'wl_end': 0.0,
        'wl_units': 'micrometers'
    }
    skip_arg = False
    dict_key = ''
    for arg in argv[1:]:
        if (skip_arg):
            if (dict_key != ''):
                config[dict_key] = arg
                dict_key = ''
                skip_arg = False
                continue
        if (arg.startswith("--help")):
            config['help_mode'] = True
        elif (arg.startswith("--version")):
            config['version_mode'] = True
        elif (arg.startswith("--no-peaks")):
            config['force_peaks'] = False
        elif (arg.startswith("--no-plots")):
            config['make_plots'] = False
        elif (arg.startswith("--in")):
            dict_key = 'input_files'
            skip_arg = True
        elif (arg.startswith("--out")):
            dict_key = 'output_files'
            skip_arg = True
        elif (arg.startswith("--step=")):
            split_arg = arg.split('=')
            config['step'] = float(split_arg[1])
        elif (arg.startswith("--threshold=")):
            split_arg = arg.split('=')
            config['threshold'] = float(split_arg[1])
        elif (arg.startswith("--wl-start=")):
            split_arg = arg.split('=')
            config['wl_start'] = float(split_arg[1])
        elif (arg.startswith("--wl-stop=")):
            split_arg = arg.split('=')
            config['wl_end'] = float(split_arg[1])
        elif (arg.startswith("--wl-units=")):
            known_units = [
                'nanometers', 'nm', 'micrometers', 'um', 'millimeters', 'mm',
                'centimeters', 'cm', 'decimeters', 'dm', 'meters', 'm'
            ]
            split_arg = arg.split('=')
            config['wl_units'] = split_arg[1]
            if (config['wl_units'] not in known_units):
                raise Exception("Unrecognized wavelength units %s!"%config['wl_units'])
        else:
            raise Exception("Unrecognized argument \"{0:s}\"!".format(arg))
    return config

## \brief Print a command-line help summary.
def print_help():
    print("###############################################           ")
    print("#                                             #           ")
    print("#          NPTM_code FILTER_CONSTANTS         #           ")
    print("#                                             #           ")
    print("###############################################           ")
    print("                                                          ")
    print("Filter the optical functions for a model material.        ")
    print("                                                          ")
    print("Usage: \"./filter_constants.py --in INPUT [options]\"     ")
    print("                                                          ")
    print("Valid options are:                                        ")
    print("--in                  Comma separated list of input files.")
    print("                      (mandatory).                        ")
    print("--help                Print this help and exit.           ")
    print("--version             Print script version and exit.      ")
    print("                                                          ")

## \brief Filter a single file based on the configuration options.
#
#  \param[in] config: `dict` A dictionary containing the script configuration.
#  \return result: `int` An integer exit code (0 if successful).
def scan_single_file(config):
    result = 0
    try:
        file_name = config['input_files']
        input_file = open(file_name, 'r')
        file_line = input_file.readline()
        num_read_lines = 1
        num_orig_data = 0
        num_filtered_data = 0
        output_name = config['output_files']
        if (output_name == ''):
            output_name = file_name.split('.')[0] + "_filtered.csv"
        split_output = output_name.split('/')
        if (len(split_output) > 1):
            output_name = split_output[-1]
        output_file = open(output_name, 'w')
        wl_factor = 1.0e6
        wl_units = config['wl_units']
        if (wl_units in ["nanometers", "nm"]):
            wl_factor = 1.0e9
        elif (wl_units in ["millimeters", "mm"]):
            wl_factor = 1.0e3
        elif (wl_units in ["centimeters", "cm"]):
            wl_factor = 1.0e2
        elif (wl_units in ["decimeters", "dm"]):
            wl_factor = 1.0e1
        elif (wl_units in ["meters", "m"]):
            wl_factor = 1.0
        wl_orig = []
        reps_orig = []
        ieps_orig = []
        wl_filtered = []
        reps_filtered = []
        ieps_filtered = []
        step = config['step']
        threshold = config['threshold']
        wl0 = 0.0
        wl1 = 0.0
        reps0 = 0.0
        reps1 = 0.0
        ieps0 = 0.0
        ieps1 = 0.0
        last_dreps = 0.0
        last_dieps = 0.0
        while (file_line != ""):
            if (file_line.startswith('#')):
                output_file.write(file_line)
                file_line = input_file.readline()
                num_read_lines += 1
                continue
            split_line = file_line.split(',')
            if (len(split_line) == 3):
                # parse the line
                if (wl0 == 0.0):
                    # always parse the first data line
                    wl0 = float(split_line[0])
                    reps0 = float(split_line[1])
                    ieps0 = float(split_line[2])
                    wl_orig.append(wl0 * wl_factor)
                    reps_orig.append(reps0)
                    ieps_orig.append(ieps0)
                    num_orig_data += 1
                    if (wl0 >= config['wl_start']):
                        breakpoint()
                        output_file.write(file_line)
                        wl_filtered.append(wl0 * wl_factor)
                        reps_filtered.append(reps0)
                        ieps_filtered.append(ieps0)
                        num_filtered_data += 1
                else:
                    can_write = True
                    wl1 = float(split_line[0])
                    reps1 = float(split_line[1])
                    ieps1 = float(split_line[2])
                    dreps = reps1 - reps0
                    dieps = ieps1 - ieps0
                    wl_orig.append(wl1 * wl_factor)
                    reps_orig.append(reps1)
                    ieps_orig.append(ieps1)
                    num_orig_data += 1
                    if (wl1 < config['wl_start']):
                        num_read_lines += 1
                        file_line = input_file.readline()
                        continue # while loop
                    if (config['wl_end'] > config['wl_start'] and wl1 > config['wl_end']):
                        break # while loop
                    if (step > 0.0):
                        if (wl1 - wl0 >= step):
                            #breakpoint()
                            # write a line if step is enabled and satisfied
                            output_file.write(file_line)
                            wl_filtered.append(wl1 * wl_factor)
                            reps_filtered.append(reps1)
                            ieps_filtered.append(ieps1)
                            num_filtered_data += 1
                            can_write = False
                            wl0 = wl1
                            reps0 = reps1
                            ieps0 = ieps1
                    # end of step > 0.0 check
                    if (config['force_peaks']):
                        rpeak = (dreps * last_dreps < 0.0)
                        ipeak = (dieps * last_dieps < 0.0)
                        if ((rpeak or ipeak) and can_write):
                            # write a line if peaks are enabled and satisfied
                            output_file.write(file_line)
                            can_write = False
                            wl_filtered.append(wl1 * wl_factor)
                            reps_filtered.append(reps1)
                            ieps_filtered.append(ieps1)
                            num_filtered_data += 1
                            wl0 = wl1
                            reps0 = reps1
                            ieps0 = ieps1
                    # end of force_peaks check
                    last_dreps = dreps
                    last_dieps = dieps
                    if (reps0 != 0.0):
                        rel_dreps = (reps0 + dreps) / reps0 if dreps > 0.0 else (reps0 - dreps) / reps0
                        if (rel_dreps < 0.0):
                            rel_dreps *= -1.0
                        if ((rel_dreps > 1.0 + threshold or rel_dreps < 1.0 - threshold) and can_write):
                            # write a line if tolerance is violated
                            output_file.write(file_line)
                            can_write = False
                            wl_filtered.append(wl1 * wl_factor)
                            reps_filtered.append(reps1)
                            ieps_filtered.append(ieps1)
                            num_filtered_data += 1
                            wl0 = wl1
                            reps0 = reps1
                            ieps0 = ieps1
                    # end of reps0 != 0.0 check
                    if (ieps0 != 0.0):
                        rel_dieps = (ieps0 + dieps) / ieps0 if dieps > 0.0 else (ieps0 - dieps) / ieps0
                        if (rel_dieps < 0.0):
                            rel_dieps *= -1.0
                        if ((rel_dieps > 1.0 + threshold or rel_dieps < 1.0 - threshold) and can_write):
                            # write a line if tolerance is violated
                            output_file.write(file_line)
                            can_write = False
                            wl_filtered.append(wl1 * wl_factor)
                            reps_filtered.append(reps1)
                            ieps_filtered.append(ieps1)
                            num_filtered_data += 1
                            wl0 = wl1
                            reps0 = reps1
                            ieps0 = ieps1
                    # end of reps0 != 0.0 check
                # end of wl0 == 0.0 check
            else:
                print("ERROR: invalid input file %s at line %d!"%(file_name, num_read_lines))
                result = 1
                break # while loop
            # end of len(split_line) check
            file_line = input_file.readline()
            num_read_lines += 1
        # end of while loop
        input_file.close()
        output_file.close()
        # INFO section
        msg_filtered = "{0:d} filtered lines".format(num_filtered_data) if num_filtered_data != 1 else "1 filtered line"
        msg_orig = "{0:d} input lines".format(num_orig_data) if num_orig_data != 1 else "1 input line"
        max_distance = get_max_distance(wl_orig, reps_orig, ieps_orig, wl_filtered, reps_filtered, ieps_filtered)
        msg_distance1 = "INFO: maximum absolute distance was {0:.5g} at {1:.5e} {2:s} in the {3:s} part".format(
            max_distance['max_difference'],
            max_distance['max_wavelength'],
            config['wl_units'],
            max_distance['differing_set']
        )
        msg_distance2 = "      (actual data value is {0:.5g}).".format(max_distance['max_value'])
        print("INFO: extracted %s out of %s."%(msg_filtered, msg_orig))
        print(msg_distance1)
        print(msg_distance2)
        # Plot making section
        if (config['make_plots']):
            plt.plot(wl_orig, reps_orig, color='red', marker='', ls='-', label=r"Original $\mathfrak{Re}(\varepsilon)$")
            plt.plot(wl_filtered, reps_filtered, color='red', marker='o', ls='', label=r"Filtered $\mathfrak{Re}(\varepsilon)$")
            plt.plot(wl_orig, ieps_orig, color='blue', marker='', ls='-', label=r"Original $\mathfrak{Im}(\varepsilon)$")
            plt.plot(wl_filtered, ieps_filtered, color='blue', marker='o', ls='', label=r"Filtered $\mathfrak{Im}(\varepsilon)$")
            plt.xlabel("Wavelength ({0:s})".format(config['wl_units']))
            plt.ylabel(r"$\mathfrak{Re}(\varepsilon)$|$\mathfrak{Im}(\varepsilon)$")
            plt.legend(loc="best")
            plt.show()
        # end plot making section
    except FileNotFoundError as ex:
        print("ERROR: file not found %s!"%config['input_files'])
        result = 1
    return result

def scan_multiple_files(config):
    result = 0
    return result

## \brief Exit code (0 for success).
exit_code = main()
exit(exit_code)