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

Implement np_sphere output parser

parent 5e3b1843
Loading
Loading
Loading
Loading
+150 −0
Original line number Diff line number Diff line
@@ -65,6 +65,12 @@ def main():
                    except Exception as ex:
                        print(ex)
                        errors += 1
                if (config['application'] == 'SPH'):
                    try:
                        errors += parse_legacy_osph(config)
                    except Exception as ex:
                        print(ex)
                        errors += 1
    return errors

## \brief Parse the command line arguments.
@@ -293,6 +299,150 @@ def parse_legacy_oclu(config):
        errors += 1
    return errors

## \brief Parse a legacy output file of np_sphere.
#
#  \param config: `dict` A dictionary containing the script configuration.
#  \return errors: `int` The number of encountered errors.
def parse_legacy_osph(config):
    errors = 0
    osph_name = config['result_file']
    root_name = config['output_name']
    out20 = None # Sphere average cross-sections
    out30 = None # Sphere integrated asymmetry parameter and radiation pressure
    out40 = None # Sphere differential radiation forces
    try:
        osph_file = open(osph_name, "r") # open the OSPH file for reading
        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")
        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")
        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")

        # Define the quantities that you need to extract
        alam = 0.0
        scaleOnXi = False
        vk = 0.0

        # Read the output file preamble
        for i in range(2):
            file_line = osph_file.readline()
        nsph = int(file_line[0:6])
        if (nsph != 1):
            errors += 1
            print("ERROR: number of spheres is not 1!")
            return errors
        for i in range(2):
            file_line = osph_file.readline()
        thifirst = float(file_line[0:11].replace("D", "E"))
        thistep = float(file_line[12:21].replace("D", "E"))
        thilast = float(file_line[22:31].replace("D", "E"))
        thsfirst = float(file_line[32:41].replace("D", "E"))
        thsstep = float(file_line[42:51].replace("D", "E"))
        thslast = float(file_line[52:61].replace("D", "E"))
        nthi = 1 if thistep == 0.0 else 1 + int((thilast - thifirst) / thistep)
        nths = 1 if thsstep == 0.0 else 1 + int((thslast - thsfirst) / thsstep)
        for i in range(2):
            file_line = osph_file.readline()
        phifirst = float(file_line[0:11].replace("D", "E"))
        phistep = float(file_line[12:21].replace("D", "E"))
        philast = float(file_line[22:31].replace("D", "E"))
        phsfirst = float(file_line[32:41].replace("D", "E"))
        phsstep = float(file_line[42:51].replace("D", "E"))
        phslast = float(file_line[52:61].replace("D", "E"))
        nphi = 1 if phistep == 0.0 else 1 + int((philast - phifirst) / phistep)
        nphs = 1 if phsstep == 0.0 else 1 + int((phslast - phsfirst) / phsstep)
        ndirs = nthi * nths * nphi * nphs

        while ("JXI =" not in file_line):
            file_line = osph_file.readline() # read the next OSPH file line
            if ("XI IS SCALE FACTOR FOR LENGTHS" in file_line):
                scaleOnXi = True
                vk = float(file_line[5:20].replace('D', 'E'))
        
        # Parsing loop until the end of the OSPH file
        while (file_line != ""):
            file_line = osph_file.readline() # read the next OSPH file line
            if (scaleOnXi):
                if (file_line.startswith("  XI=")):
                    xi = float(file_line[5:20].replace('D', 'E'))
                    alam = 2.0 * math.pi * xi / vk
            else:
                if (file_line.startswith("  VK=")):
                    # we found VK, so we calculate lambda
                    # extract VK as a number from a string section, after
                    # replacing FORTRAN's D with E
                    vk = float(file_line[5:20].replace("D", "E"))
                    alam = 2.0 * math.pi / vk
            if ("SPHERE  1" 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 = osph_file.readline()
                    if ("----- SCS ----- ABS ----- EXS ----- ALBEDS --" in file_line):
                        file_line = osph_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 (out20 is not None): out20.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 = osph_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 (out30 is not None): out30.write(output_line)
                        found_averages = True # terminate the inner loop
                # the averages were written. We look for the differential
                # section
                found_differentials = False
                while (not found_differentials):
                    for di in range(ndirs):
                        while ("JTH =" not in file_line):
                            file_line = osph_file.readline()
                        file_line = osph_file.readline()
                        tidg = float(file_line[7:17].replace("D", "E"))
                        pidg = float(file_line[24:34].replace("D", "E"))
                        tsdg = float(file_line[41:51].replace("D", "E"))
                        psdg = float(file_line[58:68].replace("D", "E"))
                        while ("SPHERE  1" not in file_line):
                            file_line = osph_file.readline()
                        if ("SPHERE  1" in file_line):
                            # we found SPHERE section. We know that the forces
                            # will be after 3 more lines
                            for i in range(3):
                                file_line = osph_file.readline()
                                # the following check is needed to parse C++ output
                                if ("INSERTION" in file_line):
                                    file_line = osph_file.readline()
                            fx = float(file_line[5:20].replace("D", "E"))
                            fy = float(file_line[25:40].replace("D", "E"))
                            fz = float(file_line[45:60].replace("D", "E"))
                            # we can write the differential values, similarly to fort.31
                            output_line = "{0:.7E},{1:.3E},{2:.3E},{3:.3E},{4:.3E},{5:.7E},{6:.7E},{7:.7E}\n".format(alam, tidg, tsdg, pidg, psdg, fx, fy, fz)
                            if (out40 is not None): out40.write(output_line)
                    found_differentials = True # terminate the inner loop
            # The parsing loop ends here

        if (out20 is not None): out20.close()
        if (out30 is not None): out30.close()
        if (out40 is not None): out40.close()
        osph_file.close() # close the OSPH file
    except Exception as ex:
        print(ex)
        errors += 1
    return errors

## \brief Print a command-line help summary.
def print_help():
    print("                                                  ")