Commit 4937eb0e authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Merge branch 'script_devel' into 'master' 

Prepare v0.10.10

See merge request giacomo.mulas/np_tmcode!117
parents fca1a3df acce932a

src/scripts/guess_resources.py

0 → 100755
+795 −0

File added.

Preview size limit exceeded, changes collapsed.

src/scripts/inertia.py

+20 −10
Original line number Diff line number Diff line
@@ -25,11 +25,15 @@ 


import math

import numpy as np

#import pdb

import pdb



from sys import argv



## \brief Main execution code

## \cond

__version__ = "0.10.10"

## \endcond



## \brief Main execution code.

#

# `main()` is the function that handles the creation of the script configuration

# and the execution of the comparison. It returns an integer value corresponding
@@ -45,7 +49,9 @@ def main(): 
        print(ex)

        print("\nType \"./inertia.py --help\" to get more detailed help.")

        errors = 1

    if config['help_mode']:

    if config['version_mode']:

      print("inertia.py v%s."%__version__)

    elif config['help_mode']:

        config['help_mode'] = True

        print_help()

    else:
@@ -131,11 +137,11 @@ def get_centers(config): 
    for i in range(nsph):

        fline = geom_file.readline()

        elems = ingest_line(fline)

        value = float(elems[0])

        value = float(elems[0].replace('D', 'E').replace('d', 'E'))

        centers['x'].append(value)

        value = float(elems[1])

        value = float(elems[1].replace('D', 'E').replace('d', 'E'))

        centers['y'].append(value)

        value = float(elems[2])

        value = float(elems[2].replace('D', 'E').replace('d', 'E'))

        centers['z'].append(value)

    geom_file.close()

    return centers
@@ -397,8 +403,8 @@ def get_types(config): 
        while (read_layers < num_layers):

            fline = sc_file.readline()

            split_line = fline.split('(')[1].split(',')

            rval = float(split_line[0])

            ival = float(split_line[1][:-2])

            rval = float(split_line[0].replace('D', 'E').replace('d', 'E'))

            ival = float(split_line[1][:-2].replace('D', 'E').replace('d', 'E'))

            dc0 = (rval, ival)

            if not dc0 in found_dc0s:

                found_dc0s.append(dc0)
@@ -456,7 +462,8 @@ def parse_arguments(): 
        'scat_name': '',

        'geom_name': '',

        'specific': [1.0],

        'help_mode': False

        'help_mode': False,

        'version_mode': False

    }

    arg_index = 1

    skip_arg = False
@@ -482,6 +489,8 @@ def parse_arguments(): 
                config['specific'].append(float(vec_specifics[si]))

        elif (arg.startswith("--help")):

            config['help_mode'] = True

        elif (arg.startswith("--version")):

            config['version_mode'] = True

        else:

            raise ValueError("Unrecognized argument \'{0:s}\'".format(arg))

        arg_index += 1
@@ -501,6 +510,7 @@ def print_help(): 
    print("--geom GEOM              Geometry configuration file (mandatory).          ")

    print("--help                   Print this help and exit.                         ")

    print("--spew=SPEC_WEIGHT       Specific weight of materials (in g/cm3, optional).")

    print("--version                Print script version and exit.                    ")

    print("                                                                           ")



# ### PROGRAM EXECUTION ###

src/scripts/model_maker.py

+10 −10
Original line number Diff line number Diff line
@@ -46,7 +46,7 @@ except ModuleNotFoundError as ex: 
from pathlib import Path

from sys import argv



## \brief Main execution code

## \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.
@@ -73,7 +73,7 @@ def main(): 
            result = 1

    return result



## \brief Populate the dielectric constant data via interpolation

## \brief Populate the dielectric constant data via interpolation.

#

#  \param sconf: `dict` Scatterer configuration dictionary.

#  \return result: `int` An exit code (0 if successful).
@@ -341,7 +341,7 @@ def load_model(model_file): 
                    for j in range(expected_radii):

                        sconf['rcf'][i][j] = float(model['particle_settings']['rad_frac'][i][j])

        # Create the gconf dict

        use_refinement = True

        use_refinement = False

        debug_am = False

        dyn_orders = True

        inv_accuracy = 1.0e-07
@@ -352,11 +352,11 @@ def load_model(model_file): 
        try:

            use_refinement = False if int(model['runtime']['refinement']) == 0 else True

        except KeyError:

            use_refinement = True

            use_refinement = False

        try:

            debug_am = False if int(model['runtime']['debug_am']) == 0 else True

        except KeyError:

            debug_am = True

            debug_am = False

        try:

            dyn_orders = False if int(model['runtime']['dyn_orders']) == 0 else True

        except KeyError:
@@ -643,7 +643,7 @@ def parse_arguments(): 
def print_help():

    print("###############################################           ")

    print("#                                             #           ")

    print("#           NPTM_code MODEL_MAKER             #           ")

    print("#           NP_TMcode MODEL_MAKER             #           ")

    print("#                                             #           ")

    print("###############################################           ")

    print("                                                          ")
@@ -679,9 +679,9 @@ def print_model_summary(scatterer, geometry): 
    for i in range(scatterer['nsph']):

        sph_type_index = scatterer['vec_types'][i] - 1

        ros = scatterer['ros'][sph_type_index]

        avgX += (ros * geometry['vec_sph_x'][i])

        avgY += (ros * geometry['vec_sph_y'][i])

        avgZ += (ros * geometry['vec_sph_z'][i])

        avgX += geometry['vec_sph_x'][i]

        avgY += geometry['vec_sph_y'][i]

        avgZ += geometry['vec_sph_z'][i]

        R3tot += math.pow(ros, 3.0)

        if (ros > Rmax):

            Rmax = ros
@@ -1430,6 +1430,6 @@ def write_obj(scatterer, geometry, max_rad): 
    os.remove(str(Path(str(out_dir), "TMP_MODEL.obj")))

    os.remove(str(Path(str(out_dir), "TMP_MODEL.mtl")))



## \brief Exit code (0 for success)

## \brief Exit code (0 for success).

exit_code = main()

exit(exit_code)

src/scripts/parse_output.py

+87 −31
Original line number Diff line number Diff line
@@ -133,25 +133,31 @@ def parse_legacy_oclu(config): 
    errors = 0

    oclu_name = config['result_file']

    root_name = config['output_name']

    out20 = None # Cluster average cross-sections

    out21 = None # Cluster average cross-sections in state -1

    out22 = None # Cluster average cross-sections in state +1

    out31 = None # Cluster differential cross-sections in state -1

    out32 = None # Cluster differential cross-sections in state +1

    out40 = None # Cluster integrated asymmetry parameter and radiation pressure

    out41 = None # Cluster integrated asymmetry parameter and radiation pressure in state -1

    out42 = None # Cluster integrated asymmetry parameter and radiation pressure in state +1

    out51 = None # Cluster differential asymmetry parameter and radiation pressure forces in state -1

    out52 = None # Cluster differential asymmetry parameter and radiation pressure forces in state +1

    oclu_file = open(oclu_name, "r") # open the OCLU file for reading

    file_line = "first line" # a string to parse the OCLU file lines

    if ('ALL' in config['selection'] or 'ICS' in config['selection']):

        out20 = open(root_name + "_ics.csv", "w") # open a file for integrated cross sections

        out20.write("Wavelength_m,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

        out21 = open(root_name + "_ics1.csv", "w") # open a file for integrated cross sections in state -1

        out22 = open(root_name + "_ics2.csv", "w") # open a file for integrated cross sections in state +1

        out21.write("Wavelength_m,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

        out22.write("Wavelength_m,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

    if ('ALL' in config['selection'] or 'DCS' in config['selection']):

        out31 = open(root_name + "_dcs1.csv", "w") # open a file for differential cross-sections in state -1

        out31.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

        out32 = open(root_name + "_dcs2.csv", "w") # open a file for differential cross-sections in state +1

        out32.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

    if ('ALL' in config['selection'] or 'IRP' in config['selection']):

        out40 = open(root_name + "_irp.csv", "w") # open a file for integrated radiation pressure forces

        out40.write("Wavelength_m,CosAv,RaPr_m2\n")

        out41 = open(root_name + "_irp1.csv", "w") # open a file for integrated radiation pressure forces in state -1

        out42 = open(root_name + "_irp2.csv", "w") # open a file for integrated radiation pressure forces in state +1

        out41.write("Wavelength_m,CosAv,RaPr_m2\n")

        out42.write("Wavelength_m,CosAv,RaPr_m2\n")

    if ('ALL' in config['selection'] or 'DRP' in config['selection']):

        out51 = open(root_name + "_drp1.csv", "w") # open a file for differential radiation pressure forces in state -1

        out51.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,CosAv,RaPr_m2,Fl_m2,Fr_m2,Fk_m2,Fx_m2,Fy_m2,Fz_m2,TQEl_m2,TQEr_m2,TQEk_m2,TQEx_m2,TQEy_m2,TQEz_m2,TQSl_m2,TQSr_m2,TQSk_m2,TQSx_m2,TQSy_m2,TQSz_m2\n")
@@ -213,13 +219,13 @@ def parse_legacy_oclu(config): 
                alam = 2.0 * math.pi / vk

        if ("CLUSTER (ENSEMBLE AVERAGE, " in file_line):

            # we are in average section. We know the average cross-sections

            # are after 3 more lines.

            # are after 3 more lines for state -1.

            for i in range(3):

                file_line = oclu_file.readline()

                # the following check is needed to parse C++ output

                if ("INSERTION" in file_line):

                    file_line = oclu_file.readline()

            # we know we are in LIN -1 because it is the first one

            # we know we are in state -1 because it is the first one

            # we also know that the next line contains the values

            # we are looking for, so we read it

            scasm = float(file_line[1:15].replace("D", "E"))
@@ -228,7 +234,7 @@ def parse_legacy_oclu(config): 
            # we can write the average values, similarly to fort.22

            # note that \n puts a new line at the end of the string

            output_line = "{0:.7E},{1:.7E},{2:.7E},{3:.7E}\n".format(alam, scasm, abssm, extsm)

            if (out20 is not None): out20.write(output_line)

            if (out21 is not None): out21.write(output_line)

            # we know that the asymmetry parameter and the radiation

            # pressure forces will be after 8 more lines

            for i in range(8):
@@ -236,7 +242,28 @@ def parse_legacy_oclu(config): 
            cosav = float(file_line[8:23].replace("D", "E"))

            rapr = float(file_line[31:46].replace("D", "E"))

            output_line = "{0:.7E},{1:.7E},{2:.7E}\n".format(alam, cosav, rapr)

            if (out40 is not None): out40.write(output_line)

            if (out41 is not None): out41.write(output_line)

            # now we look for state +1 cross-sections, after 4 lines

            for i in range(4):

                file_line = oclu_file.readline()

                # the following check is needed to parse C++ output

                if ("INSERTION" in file_line):

                    file_line = oclu_file.readline()

            scasm = float(file_line[1:15].replace("D", "E"))

            abssm = float(file_line[17:30].replace("D", "E"))

            extsm = float(file_line[32:45].replace("D", "E"))

            # we can write the average values, similarly to fort.22

            # note that \n puts a new line at the end of the string

            output_line = "{0:.7E},{1:.7E},{2:.7E},{3:.7E}\n".format(alam, scasm, abssm, extsm)

            if (out22 is not None): out22.write(output_line)

            # we know that the asymmetry parameter and the radiation

            # pressure forces will be after 8 more lines

            for i in range(8):

                file_line = oclu_file.readline()

            cosav = float(file_line[8:23].replace("D", "E"))

            rapr = float(file_line[31:46].replace("D", "E"))

            output_line = "{0:.7E},{1:.7E},{2:.7E}\n".format(alam, cosav, rapr)

            if (out42 is not None): out42.write(output_line)

            # the averages were written. We look for differential section

            for di in range(ndirs):

                while ("JTH =" not in file_line):
@@ -369,10 +396,12 @@ def parse_legacy_oclu(config): 
            # closes di for loop

        # The parsing loop ends here



    if (out20 is not None): out20.close()

    if (out21 is not None): out21.close()

    if (out22 is not None): out22.close()

    if (out31 is not None): out31.close()

    if (out32 is not None): out32.close()

    if (out40 is not None): out40.close()

    if (out41 is not None): out41.close()

    if (out42 is not None): out42.close()

    if (out51 is not None): out51.close()

    if (out52 is not None): out52.close()

    oclu_file.close() # close the OCLU file
@@ -386,31 +415,37 @@ def parse_legacy_oinclu(config): 
    errors = 0

    oclu_name = config['result_file']

    root_name = config['output_name']

    out20 = None # Inclusion average cross-sections

    out21 = None # Inclusion average cross-sections in state -1

    out22 = None # Inclusion average cross-sections in state +1

    out31 = None # Inclusion differential cross-sections in state -1

    out32 = None # Inclusion differential cross-sections in state +1

    out40 = None # Inclusion integrated asymmetry parameter and radiation pressure

    out41 = None # Inclusion integrated asymmetry parameter and radiation pressure in state -1

    out42 = None # Inclusion integrated asymmetry parameter and radiation pressure in state +1

    out51 = None # Inclusion differential asymmetry parameter and radiation pressure forces in state -1

    out52 = None # Inclusion differential asymmetry parameter and radiation pressure forces in state +1

    try:

        oclu_file = open(oclu_name, "r") # open the OINCLU file for reading

        file_line = "first line" # a string to parse the OINCLU file lines

        if ('ALL' in config['selection'] or 'ICS' in config['selection']):

            out20 = open(root_name + "_ics.csv", "w") # open a file for integrated cross sections

            out20.write("Wavelength,ScaSec,AbsSec,ExtSec\n")

            out21 = open(root_name + "_ics1.csv", "w") # open a file for integrated cross sections in state -1

            out21.write("Wavelength_m,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

            out22 = open(root_name + "_ics2.csv", "w") # open a file for integrated cross sections in state +1

            out22.write("Wavelength_m,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

        if ('ALL' in config['selection'] or 'DCS' in config['selection']):

            out31 = open(root_name + "_dcs1.csv", "w") # open a file for differential cross-sections in state -1

            out31.write("Wavelength,THi,THs,PHi,PHs,ScaSec,AbsSec,ExtSec\n")

            out31.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

            out32 = open(root_name + "_dcs2.csv", "w") # open a file for differential cross-sections in state +1

            out32.write("Wavelength,THi,THs,PHi,PHs,ScaSec,AbsSec,ExtSec\n")

            out32.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

        if ('ALL' in config['selection'] or 'IRP' in config['selection']):

            out40 = open(root_name + "_irp.csv", "w") # open a file for integrated radiation pressure forces

            out40.write("Wavelength,CosAv,RaPr\n")

            out41 = open(root_name + "_irp1.csv", "w") # open a file for integrated radiation pressure forces in state -1

            out41.write("Wavelength_m,CosAv,RaPr_m2\n")

            out42 = open(root_name + "_irp1.csv", "w") # open a file for integrated radiation pressure forces in state +1

            out42.write("Wavelength_m,CosAv,RaPr_m2\n")

        if ('ALL' in config['selection'] or 'DRP' in config['selection']):

            out51 = open(root_name + "_drp1.csv", "w") # open a file for differential radiation pressure forces in state -1

            out51.write("Wavelength,THi,THs,PHi,PHs,CosAv,RaPr,Fl,Fr,Fk,Fx,Fy,Fz,TQEl,TQEr,TQEk,TQEx,TQEy,TQEz,TQSl,TQSr,TQSk,TQSx,TQSy,TQSz\n")

            out51.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,CosAv,RaPr_m2,Fl_m2,Fr_m2,Fk_m2,Fx_m2,Fy_m2,Fz_m2,TQEl_m2,TQEr_m2,TQEk_m2,TQEx_m2,TQEy_m2,TQEz_m2,TQSl_m2,TQSr_m2,TQSk_m2,TQSx_m2,TQSy_m2,TQSz_m2\n")

            out52 = open(root_name + "_drp2.csv", "w") # open a file for differential radiation pressure forces in state +1

            out52.write("Wavelength,THi,THs,PHi,PHs,CosAv,RaPr,Fl,Fr,Fk,Fx,Fy,Fz,TQEl,TQEr,TQEk,TQEx,TQEy,TQEz,TQSl,TQSr,TQSk,TQSx,TQSy,TQSz\n")

            out52.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,CosAv,RaPr_m2,Fl_m2,Fr_m2,Fk_m2,Fx_m2,Fy_m2,Fz_m2,TQEl_m2,TQEr_m2,TQEk_m2,TQEx_m2,TQEy_m2,TQEz_m2,TQSl_m2,TQSr_m2,TQSk_m2,TQSx_m2,TQSy_m2,TQSz_m2\n")



        # Define the quantities that you need to extract

        alam = 0.0
@@ -463,14 +498,14 @@ def parse_legacy_oinclu(config): 
                    # replacing FORTRAN's D with E

                    vk = float(file_line[5:20].replace("D", "E"))

                    alam = 2.0 * math.pi / vk

            if ("ENSEMBLE AVERAGE, MODE 0" in file_line):

            if ("ENSEMBLE AVERAGE, MODE" in file_line):

                # we are in average section. We start a nested loop to

                # extract the average values

                found_averages = False

                while (not found_averages):

                    file_line = oclu_file.readline()

                    if ("----- SCS ----- ABS ----- EXS ----- ALBEDS --" in file_line):

                        # we know we are in LIN -1 because it is the first one

                        # we know we are in state -1 because it is the first one

                        # we also know that the next line contains the values

                        # we are looking for, so we read it

                        file_line = oclu_file.readline()
@@ -481,7 +516,7 @@ def parse_legacy_oinclu(config): 
                        # we can write the average values, similarly to fort.22

                        # note that \n puts a new line at the end of the string

                        output_line = "{0:.7E},{1:.7E},{2:.7E},{3:.7E}\n".format(alam, scasm, abssm, extsm)

                        if (out20 is not None): out20.write(output_line)

                        if (out21 is not None): out21.write(output_line)

                        # we know that the asymmetry parameter and the radiation

                        # pressure forces will be after 5 more lines

                        for i in range(5):
@@ -489,7 +524,26 @@ def parse_legacy_oinclu(config): 
                        cosav = float(file_line[8:23].replace("D", "E"))

                        rapr = float(file_line[31:46].replace("D", "E"))

                        output_line = "{0:.7E},{1:.7E},{2:.7E}\n".format(alam, cosav, rapr)

                        if (out40 is not None): out40.write(output_line)

                        if (out41 is not None): out41.write(output_line)

                        # we know that state +1 starts after 4 more lines

                        for i in range(4):

                            file_line = oclu_file.readline()

                        # we now extract the values from string sections

                        scasm = float(file_line[1:15].replace("D", "E"))

                        abssm = float(file_line[17:30].replace("D", "E"))

                        extsm = float(file_line[32:45].replace("D", "E"))

                        # we can write the average values, similarly to fort.22

                        # note that \n puts a new line at the end of the string

                        output_line = "{0:.7E},{1:.7E},{2:.7E},{3:.7E}\n".format(alam, scasm, abssm, extsm)

                        if (out22 is not None): out22.write(output_line)

                        # we know that the asymmetry parameter and the radiation

                        # pressure forces will be after 5 more lines

                        for i in range(5):

                            file_line = oclu_file.readline()

                        cosav = float(file_line[8:23].replace("D", "E"))

                        rapr = float(file_line[31:46].replace("D", "E"))

                        output_line = "{0:.7E},{1:.7E},{2:.7E}\n".format(alam, cosav, rapr)

                        if (out42 is not None): out42.write(output_line)

                        found_averages = True # terminate the inner loop

                # the averages were written. We look for SINGLE SCATTERER section

                # using another inner loop
@@ -609,10 +663,12 @@ def parse_legacy_oinclu(config): 
                    found_differentials = True # terminate the inner loop

            # The parsing loop ends here



        if (out20 is not None): out20.close()

        if (out21 is not None): out21.close()

        if (out22 is not None): out22.close()

        if (out31 is not None): out31.close()

        if (out32 is not None): out32.close()

        if (out40 is not None): out40.close()

        if (out41 is not None): out41.close()

        if (out42 is not None): out42.close()

        if (out51 is not None): out51.close()

        if (out52 is not None): out52.close()

        oclu_file.close() # close the OINCLU file
@@ -637,13 +693,13 @@ def parse_legacy_osph(config): 
        file_line = "first line" # a string to parse the OSPH file lines

        if ('ALL' in config['selection'] or 'ICS' in config['selection']):

            out20 = open(root_name + "_ics.csv", "w") # open a file for integrated cross sections

            out20.write("Wavelength,ScaSec,AbsSec,ExtSec\n")

            out20.write("Wavelength_m,ScaSec_m2,AbsSec_m2,ExtSec_m2\n")

        if ('ALL' in config['selection'] or 'IRP' in config['selection']):

            out30 = open(root_name + "_irp.csv", "w") # open a file for integrated radiation pressure forces

            out30.write("Wavelength,CosAv,RaPr\n")

            out30.write("Wavelength_m,CosAv,RaPr_m2\n")

        if ('ALL' in config['selection'] or 'DRP' in config['selection']):

            out40 = open(root_name + "_drp.csv", "w") # open a file for differential radiation pressure forces in state -1

            out40.write("Wavelength,THi,THs,PHi,PHs,Fx,Fy,Fz\n")

            out40.write("Wavelength_m,THi_deg,THs_deg,PHi_deg,PHs_deg,Fx_m2,Fy_m2,Fz_m2\n")



        # Define the quantities that you need to extract

        alam = 0.0

src/scripts/pywiscombe.py

+13 −2
Original line number Diff line number Diff line
@@ -29,6 +29,11 @@ 
import math

from sys import argv



## \cond

__version__ = "0.10.10"

## \endcond





## \brief Main execution code

#

# `main()` is the function that handles the creation of the script configuration
@@ -41,7 +46,9 @@ def main(): 
    errors = 0

    try:

        config = parse_arguments()

        if config['help_mode']:

        if config['version_mode']:

            print("pywiscombe.py v%s."%__version__)

        elif config['help_mode']:

            print_help()

        else:

            if (config['mode'] == ''):
@@ -83,7 +90,8 @@ def parse_arguments(): 
        'wavelength': 0.0,

        'radius': 0.0,

        'refraction' : 1.0,

        'help_mode': False

        'help_mode': False,

        'version_mode': False

    }

    arg_index = 1

    skip_arg = False
@@ -104,6 +112,8 @@ def parse_arguments(): 
            config['refraction'] = float(split_arg[1])

        elif (arg.startswith("--help")):

            config['help_mode'] = True

        elif (arg.startswith("--version")):

            config['version_mode'] = True

        else:

            raise ValueError("Unrecognized argument \'{0:s}\'".format(arg))

        arg_index += 1
@@ -124,6 +134,7 @@ def print_help(): 
    print("--wave=WAVELENGTH        Radiation wavelength in meters (mandatory).               ")

    print("--rad=RADIUS             Particle radius in meters (mandatory).                    ")

    print("--help                   Print this help and exit.                                 ")

    print("--version                Print script version and exit.                            ")

    print("                                                                                   ")



# ### PROGRAM EXECUTION ###