modules/post/d3plot/d3plot_automotive_assessment.mjs

import { D3PlotHelper } from "./d3plot.mjs";
import { JSPath } from "../../../modules/shared/path.mjs";
import { WorkflowUnits } from "../../../../modules/units.mjs";
import { Protocol, Protocols } from "../../shared/protocols.mjs";
import { ProtocolVehicle } from "../../shared/vehicle.mjs";
import { Structure } from "../../shared/structure.mjs";
import { AssessmentType } from "../../shared/assessment_types.mjs";

import { OasysProgram } from "../../../../modules/oasys.mjs";

export { ProtocolAssessment, DoStructureAssessmentStructureData };

/**
 * Class to do assessments for a protocol on D3PLOT<br><br>
 * It extends the Protocol class and adds the assessment functions<br><br>
 * It is separated like this because it uses D3PLOT JS-API functions
 * that if they were in the Protocol class would stop it from
 * working in PRIMER.
 * @extends Protocol
 * @param {string} regulation Regulation, e.g. Regulation.CNCAP
 * @param {string} crash_test Crash test, e.g. CrashTest.ODB
 * @param {string} version Version, e.g. "7.1.2"
 * @param {ProtocolVehicle} vehicle class containing info about valid occupant types
 * @example
 * let vehicle = Protocols.GetProtocolVehicle(regulation, crash_test, version);
 * let protocol = new ProtocolAssessment(Regulation.CNCAP, CrashTest.ODB, "7.1.2", vehicle, []);
 */
class ProtocolAssessment extends Protocol {
    constructor(regulation, crash_text, version, vehicle) {
        super(regulation, crash_text, version, vehicle, []);
    }

    /**
     * Create a ProtocolAssessment instance with the default Vehicle for
     * the given regulation, test and year
     * @param {string} regulation Protocol regulation
     * @param {string} crash_test Protocol test
     * @param {string} version Version, e.g. "7.1.2"
     * @returns {ProtocolAssessment}
     * @example
     * let p = ProtocolAssessment.CreateDefaultProtocol(Regulation.CNCAP, CrashTest.ODB, "7.1.2");
     */
    static CreateDefaultProtocol(regulation, crash_test, version) {
        let vehicle = Protocols.GetProtocolVehicle(regulation, crash_test, version);

        return new ProtocolAssessment(regulation, crash_test, version, vehicle);
    }

    /* Instance methods */

    /**
     * Do calculations on a Structure
     * @param {DoStructureAssessmentStructureData[]} structures_data Array of DoStructureAssessmentStructureData instances
     * @param {string[]} assessment_types Array of assessment types, e.g. AssessmentType.B_PILLAR_INTRUSION
     * @returns {DoAssessmentResults}
     * @example
     * let p = new ProtocolAssessment(Regulation.CNCAP, CrashTest.ODB, "7.1.2", []);
     * let m = Model.GetFromID(1);
     * let s = new Structure(Structure.B_PILLAR, []);
     * let extra = { "cut_section_method": "Constant X", "cut_section_nodes": [1000, 2000, 3000] }  // ...etc
     * let structures_data = [new DoStructureAssessmentStructureData(s, m, Workflow.UNIT_SYSTEM_U2, extra)];
     * let assessment_types = [AssessmentType.B_PILAR_INTRUSION];
     * let options = new DoAssessmentOptions("separate_pages");
     *
     * let assessment = p.DoStructureAssessment(structures_data,
     *                                         assessment_types);
     */
    DoStructureAssessment(structures_data, assessment_types) {
        if (!(structures_data instanceof Array)) {
            throw new Error(`<structures_data> must be an array in <ProtocolAssessment.DoStructureAssessment>`);
        }

        for (let structure_data of structures_data) {
            if (!(structure_data instanceof DoStructureAssessmentStructureData)) {
                throw new Error(
                    "Not all of the items in the <structures_data> array are DoStructureAssessmentStructureData instances in <ProtocolAssessment.DoStructureAssessment>"
                );
            }
        }

        /* Object to return values from assessments */

        let results = new DoAssessmentResults();

        for (let assessment_type of assessment_types) {
            /** @type {StructureAssessmentOutput} */
            let output = null;

            for (let structure_data of structures_data) {
                let structure = structure_data.structure;
                let model_id = structure_data.model_id;
                let unit_system = structure_data.unit_system;
                let extra = structure_data.extra;

                /* Only do the assessment if it's applicable to this structure */

                let structure_assessment_types = this.StructureAssessmentTypes(structure);

                if (!structure_assessment_types.includes(assessment_type)) continue;

                if (assessment_type == AssessmentType.B_PILLAR_INTRUSION) {
                    output = this.DoBPillarAssessment(model_id, extra, unit_system);
                } else if (assessment_type == AssessmentType.HEAD_EXCURSION) {
                    output = this.DoHeadExcursionAssessment(model_id, extra, unit_system);
                } else {
                    ErrorMessage(
                        `Unknown assessment type '${assessment_type}' in <ProtocolAssessment.DoOccupantAssessment>`
                    );
                }

                if (output) {
                    for (let p in output.values) {
                        results.output[`M${model_id} ${structure} ${p}`] = output.values[p];
                    }

                    results.image_filenames.push(output.image_filename);
                } else {
                    WarningMessage(
                        `No output from assessment type '${assessment_type}' for structure M${model_id} '${structure}'`
                    );
                }
            }
        }

        return results;
    }

    /**
     * Returns a list of the assessment types for the given structures
     * @param {Structure[]} structures Structure instances
     * @returns {string[]}
     */
    UniqueStructureAssessmentTypes(structures) {
        let assessment_types = new Set();

        for (let structure of structures) {
            let at = this.StructureAssessmentTypes(structure);

            for (let a of at) {
                assessment_types.add(a);
            }
        }

        return [...assessment_types];
    }

    /**
     * Object to return from the structure assesment functions
     * @typedef StructureAssessmentOutput
     * @property {Object} values Object with values for the output
     * @property {string} image_filename Filename of an image
     */

    /**
     * Object B-Pillar data
     * @typedef {Object} BPillarStructure
     * @property {string} cut_section_method Cut section method
     * @property {number[]} cut_section_nodes Cut section nodes
     * @property {number[]} pre_crash_parts Pre-crash parts
     * @property {number[]} post_crash_parts Post-crash parts
     * @property {number[]} shift_deform_nodes Shift deform nodes
     * @property {number} ground_z Ground Z coordinate
     * @property {number} seat_centre_y Seat centre Y coordinate
     * @property {number} h_point_z H-Point Z coordinate
     */

    /**
     * Does the B-Pillar intrusion assessment
     * @param {number} model_id Model to do the assessment on
     * @param {BPillarStructure} b_pillar B-Pillar structure with data required for assessment
     * @param {number} unit_system Unit system
     * @returns {?StructureAssessmentOutput}
     */
    DoBPillarAssessment(model_id, b_pillar, unit_system) {
        /* Directories and files */

        let template_dir = `${JSPath.GetWorkflowsDirectory()}scripts/automotive_assessments/post/reporter/templates`;

        let first_state_coords_filename = `${template_dir}/first_state_coords.txt`;
        let last_state_coords_filename = `${template_dir}/last_state_coords.txt`;

        let template_file = `${template_dir}/b_pillar_assessment_image.ortx`;
        let log_file = `${template_dir}/b_pillar.log`;

        let image_filename = `${template_dir}/b_pillar_image_M${model_id}.png`;
        let variables_filename = `${template_dir}/b_pillar_variables`;

        Message(`Writing M${model_id} B-Pillar cut section coordinates to file...`);

        if (
            !this.WriteBPillarCutSectionCoordinates(
                model_id,
                b_pillar,
                unit_system,
                first_state_coords_filename,
                last_state_coords_filename
            )
        ) {
            WarningMessage(
                `Unable to write M${model_id} B-Pillar cut section coordinates to file in <ProtocolAssessment.Write.`
            );
            return null;
        }

        /* Run a REPORTER template to read the coordinates and create the image */

        Message(`Running REPORTER template to create M${model_id} B-Pillar image...`);

        if (!this.RunBPillarReporterTemplate(model_id, b_pillar, unit_system, template_file, log_file)) {
            WarningMessage(`Unable to run REPORTER template to create M${model_id} B-Pillar image.`);
            return null;
        }

        /* Get values from the REPORTER variables file */

        Message("Getting values from REPORTER variables file...");

        let variables = this.GetReporterVariableValues(variables_filename);

        if (!variables) {
            WarningMessage("Unable to get values from REPORTER variables file.");
            return null;
        }

        /* Delete unwanted files */

        let filenames = [first_state_coords_filename, last_state_coords_filename, variables_filename, log_file];

        for (let filename of filenames) {
            if (File.Exists(filename) && File.IsFile(filename)) {
                File.Delete(filename);
            }
        }

        /* Get here and everything seems to have worked so return the values and the image filename */

        let values = {};

        for (let variable of variables) {
            if (variable.name == "MAX_INTRUSION" || variable.name == "HEIGHT_ABOVE_GROUND") {
                values[variable.name] = `${parseFloat(variable.value).toPrecision(6)}mm`;
            } else {
                values[variable.name] = variable.value;
            }
        }

        return {
            values: values,
            image_filename: image_filename
        };
    }

    /**
     * Writes the cut section coordinates to two files; one at the first state, one at the last state.
     * The files are written to the same directory as the B-Pillar template.
     * @param {number} model_id Model to do the assessment on
     * @param {BPillarStructure} b_pillar B-Pillar structure with data required for assessment
     * @param {number} unit_system Unit system
     * @param {string} first_state_coords_filename Filename to write first state coordinates to
     * @param {string} last_state_coords_filename Filename to write last state coordinates to
     * @return {boolean}
     */
    WriteBPillarCutSectionCoordinates(
        model_id,
        b_pillar,
        unit_system,
        first_state_coords_filename,
        last_state_coords_filename
    ) {
        let window_id = D3PlotHelper.FirstWindowWithModel(model_id);

        if (!window_id) {
            WarningMessage(
                `Unable to find window with model ${model_id} in <ProtocolAssessment.WriteBPillarCutSectionCoordinates>`
            );
            return false;
        }

        let models = GetWindowModels(window_id);

        if (models.nm > 1) {
            WarningMessage(
                `Model M${model_id} shares a window with other models. It needs to be in a window by itself.`
            );
            return false;
        }

        /* Create cut section */

        if (
            !this.SetCutSectionByMethod(
                model_id,
                b_pillar.cut_section_method,
                0.0,
                b_pillar.cut_section_nodes[0],
                b_pillar.cut_section_nodes[1],
                b_pillar.cut_section_nodes[2]
            )
        )
            return;

        /* Shift deform nodes */

        this.SetBPillarShiftDeform(model_id, b_pillar);

        /* Get internal part IDS */

        let n_parts = GetNumberOf(PART);

        let internal_pre_crash_pids = [];
        let internal_post_crash_pids = [];

        for (let i = 1; i <= n_parts; i++) {
            let label = GetLabel(PART, i);

            for (let pre_crash_part of b_pillar.pre_crash_parts) {
                if (label == pre_crash_part) internal_pre_crash_pids.push(i);
            }

            for (let post_crash_part of b_pillar.post_crash_parts) {
                if (label == post_crash_part) internal_post_crash_pids.push(i);
            }
        }

        /* Scale factor to convert coordinates to mm (this is what the REPORTER template expects) */

        let scale_factor = 1.0 / WorkflowUnits.LengthToMillimetresFactor(unit_system);

        /* Loop over all shells in model.  If it is a shell in a selected parts
         * get the coordinates where the cut-section cuts it.
         *
         * Two passes - one at first state, one at last */

        let n_shells = GetNumberOf(SHELL);
        let n_states = GetNumberOf(STATE);
        /** @type {number[]} */
        let internal_part_ids;
        /** @type {string} */
        let filename;

        for (let pass = 0; pass < 2; pass++) {
            if (pass == 0) {
                Message(`Writing M${model_id} cut section coordinates at first state...`);
                SetCurrentState(1);
                internal_part_ids = internal_pre_crash_pids;
                filename = first_state_coords_filename;
            }
            if (pass == 1) {
                Message(`Writing M${model_id} cut section coordinates at last state...`);
                SetCurrentState(n_states);
                internal_part_ids = internal_post_crash_pids;
                filename = last_state_coords_filename;
            }

            let file = new File(filename, File.WRITE);

            for (let i = 1; i <= n_shells; i++) {
                let topology = GetTopology(SHELL, i);

                let valid = false;

                for (let internal_part_id of internal_part_ids) {
                    if (internal_part_id == topology.pid) {
                        valid = true;
                        break;
                    }
                }

                if (!valid) continue;

                /* Shell in one of the selected parts. Get the cut coords. */

                let o = GetCutCoords({ window_id: window_id, model_id: model_id, type_code: SHELL, item: i });

                let number_of_cuts = o.n;

                if (number_of_cuts <= 0) continue;

                file.Write(number_of_cuts);

                for (let j = 0; j < number_of_cuts; j++) {
                    file.Write(`,${o.x[j] * scale_factor},${o.y[j] * scale_factor},${o.z[j] * scale_factor}`);
                }

                file.Write("\n");
            }

            file.Close();
        }

        return true;
    }

    /**
     *
     * @param {number} model_id Model to set cut section on
     * @param {string} cut_section_method Cut section method ("Three Nodes", "Constant X", "Constant Y", "Constant Z")
     * @param {number} thickness Cut section thickness
     * @param {number} cut_section_node_1 Node ID of first node
     * @param {number} cut_section_node_2 Node ID of second node
     * @param {number} cut_section_node_3 Node ID of first node
     * @return {boolean}
     */
    SetCutSectionByMethod(
        model_id,
        cut_section_method,
        thickness,
        cut_section_node_1,
        cut_section_node_2,
        cut_section_node_3
    ) {
        let window_id = D3PlotHelper.FirstWindowWithModel(model_id);

        if (!window_id) {
            WarningMessage(
                `Unable to find window with model ${model_id} in <ProtocolAssessment.SetCutSectionByMethod>`
            );
            return false;
        }

        if (cut_section_method == "Three Nodes") {
            return SetCutSection({
                window_id: window_id,
                definition: N3,
                nodes: [-cut_section_node_1, -cut_section_node_2, -cut_section_node_3],
                space: BASIC,
                status: ON,
                thickness: thickness
            });
        } else {
            let definition = CONST_X;

            if (cut_section_method == "Constant X") {
                definition = CONST_X;
            } else if (cut_section_method == "Constant Y") {
                definition = CONST_Y;
            } else if (cut_section_method == "Constant Z") {
                definition = CONST_Z;
            } else {
                ErrorMessage(
                    `Invalid cut section method in <ProtocolAssessment.SetCutSectionByMethod>: ${cut_section_method}`
                );
                return false;
            }

            return SetCutSection({
                window_id: window_id,
                definition: definition,
                nodes: [-cut_section_node_1],
                space: BASIC,
                status: ON,
                thickness: thickness
            });
        }
    }

    /**
     * Set the shift deform nodes for the B-Pillar structure.
     * @param {number} model_id Model to do the assessment on
     * @param {BPillarStructure} b_pillar B-Pillar structure with data required for assessment
     */
    SetBPillarShiftDeform(model_id, b_pillar) {
        DialogueInput(
            `/CM ${model_id}`,
            "/DEFORM SHIFT_DEFORMED DEFINE",
            b_pillar.shift_deform_nodes[0].toString(),
            b_pillar.shift_deform_nodes[1].toString(),
            b_pillar.shift_deform_nodes[2].toString()
        );
    }
    /**
     * Runs the Reporter template to generate an image for the B-Pillar assessment using the
     * coordinates written out to files by D3PLOT.
     * Returns true if the template ran successfully, false otherwise.
     * @param {number} model_id Model to do the assessment on
     * @param {BPillarStructure} b_pillar B-Pillar structure with data required for assessment
     * @param {number} unit_system Unit system
     * @param {string} template_file Template file to run
     * @param {string} log_file Log file to write to
     * @return {boolean}
     */
    RunBPillarReporterTemplate(model_id, b_pillar, unit_system, template_file, log_file) {
        /* Get the reporter executable */

        let reporter_exe = OasysProgram.GetOaInstallExecutable(OasysProgram.REPORTER);

        if (!reporter_exe) {
            ErrorMessage(`Cannot find reporter executable in <RunBPillarReporterTemplate>.`);
            return false;
        }

        /* Scale factor to convert coordinates to mm (this is what the REPORTER template expects) */

        let scale_factor = 1.0 / WorkflowUnits.LengthToMillimetresFactor(unit_system);

        /* Variables to pass to template */

        let ground_z = b_pillar.ground_z * scale_factor;
        let seat_centre_y = b_pillar.seat_centre_y * scale_factor;
        let h_point_z = b_pillar.h_point_z * scale_factor;

        let var_ground_z = `-varGROUND_Z=${ground_z}`;
        let var_seat_centre_y = `-varSEAT_CENTRE_Y=${seat_centre_y}`;
        let var_h_point_z = `-varH_POINT_Z=${h_point_z}`;
        let var_model_id = `-varMODEL_ID=${model_id}`;

        /* Limits are in cm */

        let good_deformation = 18;
        let acceptable_deformation = 14;
        let marginal_deformation = 10;

        let var_good_deformation = `-varB_PILLAR_DEFORMATION_GOOD=${good_deformation}`;
        let var_acceptable_deformation = `-varB_PILLAR_DEFORMATION_ACCEPTABLE=${acceptable_deformation}`;
        let var_marginal_deformation = `-varB_PILLAR_DEFORMATION_MARGINAL=${marginal_deformation}`;

        /* TODO - option to select SEAT or VEHICLE from GUI */
        let var_centerline = `-varCENTERLINE=SEAT`;

        /* Run the template in batch to generate the image
         *
         * Use System() instead of Execute() because I couldn't get Execute() to work
         * when file paths contained spaces.
         */

        let cmd = `""${reporter_exe}" \
               -file="${template_file}" \
               ${var_ground_z} \
               ${var_seat_centre_y} \
               ${var_h_point_z} \
                ${var_model_id} \
               ${var_good_deformation} \
               ${var_acceptable_deformation} \
               ${var_marginal_deformation} \
               ${var_centerline} \
               -log="${log_file}" \
               -batch \
               -generate \
               -iconise \
               -exit"`;

        System(cmd);

        return true;
    }

    /**
     * Object to store Reporter variables
     * @typedef {Object} ReporterVariable
     * @property {string} name
     * @property {string} description
     * @property {string} value
     */

    /**
     * Reads the variables file written by REPORTER and returns an array
     * of objects containing the values of the variables. Returns null on error.
     * @param {string} filename Variables file to read
     * @returns {?ReporterVariable[]}
     */
    GetReporterVariableValues(filename) {
        if (!File.Exists(filename)) {
            ErrorMessage("Unable to find the variables file written by the REPORTER template.");
            return null;
        }

        let file = new File(filename, File.READ);

        /* Each line is of the form
         *
         * VAR <variable name> DESCRIPTION='<description>' VALUE='<value>'
         */

        /** @type {ReporterVariable[]} */
        let variables = [];

        let line = "";
        while ((line = file.ReadLongLine()) != undefined) {
            let regex = /VAR\s+(\w+)\s+DESCRIPTION='(.+)'\s*VALUE='(.+)'/;

            let match = line.match(regex);

            if (!match) {
                ErrorMessage(`Unable to parse line in variables file: ${line}`);
                return null;
            }

            variables.push({
                name: match[1],
                description: match[2],
                value: match[3]
            });
        }

        file.Close();

        return variables;
    }

    /**
     * Object to store data for Head Excursion
     * @typedef {Object} HeadExcursionStructure
     * @property {number} cut_section_thickness Cut section thickness
     * @property {number} cut_section_node Cut section node
     * @property {string} vehicle_direction Vehicle direction (either 'positive X' or 'negative X')
     * @property {number[]} head_parts Head parts
     * @property {number[]} barrier_parts Barrier parts
     * @property {number[]} shift_deform_nodes Shift deform nodes
     * @property {number} seat_centre_y Seat centre Y coordinate
     * @property {number} intrusion_from_seat_centre_y Intrusion From Seat centre Y coordinate
     * @property {string} countermeasure Countermeasure (either 'Yes' or 'No')
     * @property {string} keyword_file Keyword filename
     */

    /**
     * Does the Head Excursion assessment
     * @param {number} model_id Model to do the assessment on
     * @param {HeadExcursionStructure} head_excursion Head Excursion structure with data required for assessment
     * @param {number} unit_system Unit system
     * @returns {?StructureAssessmentOutput}
     */
    DoHeadExcursionAssessment(model_id, head_excursion, unit_system) {
        let EuroNCAPColourGREEN = "R38G155B41";
        let EuroNCAPColourYELLOW = "R255G204B0";
        let EuroNCAPColourORANGE = "R255G151B0";
        let EuroNCAPColourBROWN = "R117G63B42";
        let EuroNCAPColourRED = "R224G7B0";

        let EuroNCAPColourREDorBROWN = EuroNCAPColourRED; //default (changes based on intrusion level and countermeasure status)

        /* Set the orientation (-YZ or YZ) based on vehicle direction and then AC
         * (do this before reading in the zone parts, so it scales to vehicle */

        if (head_excursion.vehicle_direction == "Negative X") {
            DialogueInput("-YZ");
        } else {
            DialogueInput("+YZ");
        }

        DialogueInput("/AC");

        /* Create a PTF file to show the excursion zones */

        Message(`Creating Head Excursion zones for model M${model_id}...`);

        if (!this.CreateHeadExcursionZones(head_excursion, unit_system)) {
            WarningMessage(`Unable to create head excursion zones for model M${model_id}`);
            return null;
        }

        /* Get the PTF file and read it into D3PLOT */

        let output_dir = `${JSPath.GetWorkflowsDirectory()}scripts/automotive_assessments/post/primer/output`;
        let ptf_file = `${output_dir}/zones.ptf`;

        Message(`Reading in PTF file containing zones for model M${model_id}...`);

        let zone_model_id = D3PlotHelper.OpenModel(ptf_file);

        if (!zone_model_id) {
            WarningMessage(`Unable to read in PTF file containing zones for model M${model_id}`);
            return null;
        }

        /* Get the window ids containing the vehicle model and zone model and then
         * move the zones model into the same window as the vehicle model and close
         * the window that originally contained the zones model */

        let window_id = D3PlotHelper.FirstWindowWithModel(model_id);

        if (!window_id) {
            WarningMessage(
                `Unable to find window with model ${model_id} in <ProtocolAssessment.DoHeadExcursionAssessment>`
            );
            return null;
        }

        let zone_window_id = D3PlotHelper.FirstWindowWithModel(zone_model_id);

        if (!D3PlotHelper.PutModelsInWindow(window_id, [model_id, zone_model_id])) {
            WarningMessage(`Unable to move zones model into same window as model M${model_id}`);
            return null;
        }

        DeleteWindow([zone_window_id]);

        /* Create cut section */

        if (
            !this.SetCutSectionByMethod(
                model_id,
                "Constant X",
                head_excursion.cut_section_thickness,
                head_excursion.cut_section_node,
                0,
                0
            )
        )
            return null;

        /* Apply the shift deformed nodes */

        D3PlotHelper.SetShiftDeformed(
            model_id,
            head_excursion.shift_deform_nodes[0],
            head_excursion.shift_deform_nodes[1],
            head_excursion.shift_deform_nodes[2]
        );

        SetCurrentState(1);

        /* Define intrusion y based on intrusion_from_seat_center_y
         *
         * if intrusion_from_seat_center_y is > 0 then the intrusion does not reach the seat center
         * if intrusion_from_seat_center_y < 0 then the intrusion passes the seat center
         */

        let sign = -1;
        let intrusion_y = null;

        if (head_excursion.seat_centre_y > 0) {
            /* left hand-drive vehicle pointing in -X direction
             * OR
             * right hand-drive vehicle pointing in +X direction */

            /* in both cases passenger seat (seat nearest barrier) is at a +ve Y coordinate and
             * the dummy head moves in the +ve Y direction on impact */

            sign = 1;
            intrusion_y = head_excursion.seat_centre_y + head_excursion.intrusion_from_seat_centre_y;
        } else {
            /* right hand-drive vehicle pointing in -X direction
             * OR
             * left hand-drive vehicle pointing in +X direction */

            /* in both cases passenger seat (seat nearest barrier) is at a -ve Y coordinate and
             * the dummy head moves in the -ve Y direction on impact */
            sign = -1;
            intrusion_y = head_excursion.seat_centre_y - head_excursion.intrusion_from_seat_centre_y;
        }

        /* For each head part get extract all the nodes and store in internal_head_node_ids Array */

        let internal_head_node_ids = Array();

        for (let head_part of head_excursion.head_parts) {
            /** @type {GetElemsInPartReturn} */
            let a = null;
            if ((a = GetElemsInPart(-head_part))) {
                let nelems = a.nn;

                //get the nodes for each element and store in internal_head_node_ids Array
                for (let i = 0; i < nelems; i++) {
                    let top = GetTopology(a.type, a.list[i]);
                    internal_head_node_ids = internal_head_node_ids.concat(top.top);
                }
            }
        }

        /* Get rid of duplicate node definitions to avoid checking the same nodes multiple times */
        internal_head_node_ids = internal_head_node_ids.filter(ProtocolAssessment.only_unique);

        /* The origin is the first shift deform node and the relative Y distance is used to calculate
         * the head excursion distance the assumption is that the shift deform nodes do not deform
         * themselves in the impact but they may move due to global motion of the vehicle as a result
         * of the impact. */

        let origin_y = GetData(BY, NODE, -head_excursion.shift_deform_nodes[0]);

        /* We arbitrarily pick the first head node as the one to use to find the state of the maximum
         * excursion. Note that because the head will likely rotate the maximum head excursion may not
         * be fully captured by the node chosen so we need to do a second search of neighbouring states
         * once we have found the node with the maximum excursion and iterate until true max is found */

        let head_internal_node_id = internal_head_node_ids[0];

        SetCurrentState(1);

        /* get matrix to rotate n1->n2 vector to global y at initial state;
         * after first state car will start to move so we need to use a rotated 'local' y
         * which can be found by multiplying new n1->n2 vector by the roation matrix calculated below
         * local coordinate */

        let v_n1n2 = ProtocolAssessment.get_normalised_vector_from_coords(
            // @ts-ignore
            GetData(BV, NODE, -head_excursion.shift_deform_nodes[0]),
            GetData(BV, NODE, -head_excursion.shift_deform_nodes[1])
        );

        /* rotation matrix initially aligned with +ve y direction */
        let rotation_matrix = ProtocolAssessment.get_rotation_matrix(v_n1n2, [0, 1, 0]);

        let initial_local_datum = ProtocolAssessment.get_local_datum(head_excursion, rotation_matrix);

        Message("Local y = " + JSON.stringify(initial_local_datum.y_axis));
        let nstates = GetNumberOf(STATE);

        //initially set maximum excursion to large negative (we later calculate excursion working backwards from last state, but the assumption is max excursion is at or near last state)
        let max_excursion = -999999;
        let max_state = nstates;
        let threshold = 0.2;

        //for each state (work backwards from the end and break out when current_excursion is lower than max excursoin by > 20% of the distance from the middle of the car to the middle of the seat)
        for (let state = nstates; state >= 1; state--) {
            //Message("State " + state);
            SetCurrentState(state);
            let current_local_datum = ProtocolAssessment.get_local_datum(head_excursion, rotation_matrix);

            let current_excursion = ProtocolAssessment.get_current_excursion_from_seat_cl(
                // @ts-ignore
                origin_y,
                current_local_datum,
                head_excursion.seat_centre_y,
                head_internal_node_id,
                state
            );

            if (current_excursion > max_excursion) {
                max_excursion = current_excursion;
                max_state = state;
                Message(
                    "excursion for node " + head_internal_node_id + " at state " + state + " = " + current_excursion
                );
            } else if ((max_excursion - current_excursion) / Math.abs(head_excursion.seat_centre_y) > threshold) {
                Message(
                    "Max State = " +
                        max_state +
                        " Current State = " +
                        state +
                        ": Max excursion = " +
                        max_excursion +
                        ": Current excursion = " +
                        current_excursion
                );
                break; //break out of loop as max excursion already found and it is a waste of compute to continue checking earlier states
            }
        }

        //now we have state at which maximum excursion occurs we set the state and we calculate excursion for all nodes in head
        Message("Max excursion for N" + GetLabel(NODE, head_internal_node_id) + " found at state " + max_state);
        SetCurrentState(max_state);

        let local_datum = ProtocolAssessment.get_local_datum(head_excursion, rotation_matrix); //this only changes with state so no need to recalculate each iteration of loop
        let peak_node_id = head_internal_node_id;

        //can start with n=1 (i.e. second node in list as we have already calculated the peak excursion for the the first node above)
        for (let n = 1; n < internal_head_node_ids.length; n++) {
            let current_excursion = ProtocolAssessment.get_current_excursion_from_seat_cl(
                // @ts-ignore
                origin_y,
                local_datum,
                head_excursion.seat_centre_y,
                internal_head_node_ids[n],
                sign
            );

            if (current_excursion > max_excursion) {
                max_excursion = current_excursion;
                peak_node_id = internal_head_node_ids[n];
            }
        }

        //now we have the node which generate the peak excursion we need to check the surrounding states to be sure that
        //we have the true peak.

        //boolean used to determine if we have found the true pea
        let peak_found = false;

        // save the max state found above as it is overwritten below
        let peak_state = max_state;
        // we need to to this because this is used when we check in the other direction i.e. backwards in states

        while (peak_state < nstates && !peak_found) {
            //check the next above state for the 'peak node'
            let state = peak_state + 1;

            SetCurrentState(state);
            let current_excursion = ProtocolAssessment.get_current_excursion_from_seat_cl(
                // @ts-ignore
                origin_y,
                ProtocolAssessment.get_local_datum(head_excursion, rotation_matrix),
                head_excursion.seat_centre_y,
                peak_node_id,
                sign
            );

            if (current_excursion > max_excursion) {
                max_excursion = current_excursion;
                peak_state = state;
            } else {
                peak_found = true;
            }
        }

        //now check in the other direction

        let excursion = 0.0;

        if (peak_state == 1) {
            WarningMessage(
                "Maximum excursion found at first state\nCheck that near seat center y is defined correctly and that dummy moves towards far side impact."
            );
            max_state = peak_state; //1
            peak_found = true;
            excursion = max_excursion;
        } else {
            peak_found = false;

            while (max_state > 1 && !peak_found) {
                //check the next below state for the 'peak node'
                let state = max_state - 1;

                SetCurrentState(state);
                let current_excursion = ProtocolAssessment.get_current_excursion_from_seat_cl(
                    // @ts-ignore
                    origin_y,
                    ProtocolAssessment.get_local_datum(head_excursion, rotation_matrix),
                    head_excursion.seat_centre_y,
                    peak_node_id,
                    sign
                );

                if (current_excursion > max_excursion) {
                    max_excursion = current_excursion;
                    max_state = state;
                    peak_state = state;
                } else {
                    excursion = max_excursion;
                    peak_found = true;
                }
            }
        }

        let excursion_zone = ProtocolAssessment.get_excursion_zone(excursion, head_excursion, unit_system);

        //peak_state will now be the true peak state

        SetCurrentState(peak_state);

        //get excursion distance
        //present max_excursion value as distance from driver seat centreline rather than distance from occupant seat centreline
        //(we can assume driver_seat_centre_y = -head_excursion.seat_centre_y so distance is 2*Math.abs(oGlblData.seat_centre_y) then add the excursion
        //value excursion is +ve if it is further from the driver

        max_excursion = 2 * head_excursion.seat_centre_y + excursion;

        // TODO - first colour will be brown or red - need to calculate it

        this.ColourHeadExcursionZoneParts(
            zone_model_id,
            EuroNCAPColourREDorBROWN,
            EuroNCAPColourORANGE,
            EuroNCAPColourYELLOW,
            EuroNCAPColourGREEN
        );

        this.ColourHeadExcursionVehicleModel(model_id, head_excursion.head_parts, head_excursion.barrier_parts);

        /* Set up the display */

        this.SetupHeadExcursionDisplay(model_id, zone_model_id, peak_state);

        /* Take an image */

        let image_filename = `${output_dir}/head_excursion_image_M${model_id}.png`;

        DialogueInput("/IMAGE WHITE_BACKGROUND ON", "PNG_24BIT", image_filename);

        /* Get here and everything seems to have worked so return the values and the image filename */

        let values = {
            max_excursion: max_excursion,
            excursion: excursion,
            excursion_zone: excursion_zone
        };

        return {
            values: values,
            image_filename: image_filename
        };
    }

    /**
     * Create NULL shell parts in PRIMER to represent the zone bands in car space (relative to the driver)
     * and then save the parts to a PTF file for D3PLOT to read in
     * @param {HeadExcursionStructure} head_excursion Head Excursion structure with data required for assessment
     * @param {number} unit_system Unit system
     * @return {boolean}
     */
    CreateHeadExcursionZones(head_excursion, unit_system) {
        let primer_exe = OasysProgram.GetOaInstallExecutable(OasysProgram.PRIMER);

        if (!primer_exe) {
            ErrorMessage(`Cannot find primer executable in <CreateHeadExcursionZones>.`);
            return false;
        }

        /* Open primer, read in model and run the script to output shells to a PTF file, passing it
         * the values it requires using -js_arg= */

        let js = `${JSPath.GetWorkflowsDirectory()}scripts/automotive_assessments/post/primer/head_excursion_zone.js`;

        let keyword_file = head_excursion.keyword_file;

        if (!File.Exists(keyword_file)) {
            ErrorMessage(`Unable to find keyword file: ${keyword_file} in <CreateHeadExcursionZones>.`);
            return false;
        }

        /* The order that these are defined needs to match the order they are read in head_excursion_zone.js.
         * If they change here, make sure to update them in head_excursion_zone.js too. */

        let js_args = [
            head_excursion.cut_section_thickness,
            head_excursion.cut_section_node,
            head_excursion.shift_deform_nodes[0],
            head_excursion.shift_deform_nodes[1],
            head_excursion.shift_deform_nodes[2],
            head_excursion.seat_centre_y,
            head_excursion.intrusion_from_seat_centre_y,
            head_excursion.vehicle_direction,
            unit_system
        ];

        let cmd = `""${primer_exe}" \
-js="${js}" \
-js_arg="${js_args[0]}" \
-js_arg="${js_args[1]}" \
-js_arg="${js_args[2]}" \
-js_arg="${js_args[3]}" \
-js_arg="${js_args[4]}" \
-js_arg="${js_args[5]}" \
-js_arg="${js_args[6]}" \
-js_arg="${js_args[7]}" \
-js_arg="${js_args[8]}" \
"${keyword_file}" \
-exit"`;

        Message(`Running PRIMER with command: ${cmd}`);
        System(cmd);

        return true;
    }

    /**
     * Colour the zone parts in a model
     * @param {number} model_id Model ID
     */
    ColourHeadExcursionZoneParts(model_id, part_2_colour, part_3_colour, part_4_colour, part_5_colour) {
        SetCurrentModel(model_id);
        let nparts = GetNumberOf(PART);

        for (let i = 1; i <= nparts; i++) {
            /* note part 1 is blanked so not important to set colour but had issue setting colour of first part called
             * "PART", "1", "BLACK", so left this in so that colour of parts that matter are set properly */
            switch (GetLabel(PART, i).toString()) {
                case "1":
                    DialogueInput(`CM ${model_id}`, "BLANK", "PART", "1");
                    break;
                case "2":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "2", part_2_colour);
                    break;
                case "3":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "3", part_3_colour);
                    break;
                case "4":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "4", part_4_colour);
                    break;
                case "5":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "5", part_5_colour);
                    break;
                case "6":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "6", "BLUE");
                    break;
                case "7":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "7", "GREY");
                    break;
                case "8":
                    DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "8", "WHITE");
                    break;
            }
        }
    }

    /**
     * Colour the vehicle, occupant head and barrier parts
     * @param {number} model_id Model ID
     * @param {number[]} head_parts Head part IDs
     * @param {number[]} barrier_parts Barrier part IDs
     */
    ColourHeadExcursionVehicleModel(model_id, head_parts, barrier_parts) {
        let D3PLOT_CAR_COLOUR = "R200G200B200";
        let D3PLOT_HEAD_COLOUR = "BLUE";
        let D3PLOT_BARRIER_COLOUR = "DARK_GREY";

        SetCurrentModel(model_id);

        /* Call shaded again to force redraw as sometimes it was blank */
        DialogueInput(`CM ${model_id}`, "SHADED");

        /* Set the colour of model (i.e. car + occupant + barrier) */
        DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", "ALL", D3PLOT_CAR_COLOUR);

        /* Colour barrier to D3PLOT_BARRIER_COLOUR = "DARK_GREY" */
        if (barrier_parts && barrier_parts.length > 0) {
            DialogueInput(
                `CM ${model_id}`,
                "PROPERTIES",
                "COLOUR",
                "PART",
                barrier_parts.join(" "),
                D3PLOT_BARRIER_COLOUR
            );
        } else {
            Message("No barrier parts provided. This may because none are prenent in the model.");
        }

        /* Colour the head to D3PLOT_HEAD_COLOUR = "BLUE" */
        if (head_parts && head_parts.length > 0) {
            DialogueInput(`CM ${model_id}`, "PROPERTIES", "COLOUR", "PART", head_parts.join(" "), D3PLOT_HEAD_COLOUR);
        } else {
            Message(`No head parts provided.`);
        }
    }

    /**
     * Sets up the display for the head excursion assessment
     * @param {number} vehicle_model_id Vehicle model ID
     * @param {number} zone_model_id Zone model ID
     * @param {number} peak_state Peak state
     */
    SetupHeadExcursionDisplay(vehicle_model_id, zone_model_id, peak_state) {
        DialogueInput("/SHADED");

        DialogueInput(`CM ${vehicle_model_id}`, "DISPLAY", "HEADER_SWITCH", "OFF");
        DialogueInput(`CM ${vehicle_model_id}`, "DISPLAY", "CLOCK_SWITCH", "OFF");

        /* Set overlay lines to off (like pressing 'y' shortcut) */

        DialogueInput(`CM ${vehicle_model_id}`, "DISPLAY", "OD", "OFF");

        DialogueInput(
            `CM ${vehicle_model_id}`,
            "DISPLAY",
            "ENTITY",
            "BEAM",
            "OFF" /* beam elements off */,
            "SPRING",
            "OFF" /* spring elements off */,
            "X",
            "OFF" /* cut section off */
        );

        /* Turn off clock, header and display lines in zones model too */

        DialogueInput(`CM ${zone_model_id}`, "DISPLAY", "HEADER_SWITCH", "OFF");
        DialogueInput(`CM ${zone_model_id}`, "DISPLAY", "CLOCK_SWITCH", "OFF");
        DialogueInput(`CM ${zone_model_id}`, "DISPLAY", "OD", "OFF");

        /* Set to the peak state */

        DialogueInput(`STATE ${peak_state}`);

        /* Autocentre model and redraw by doing a shaded plot */

        DialogueInput(`CM ${vehicle_model_id}`, "AC");
        DialogueInput("/SH");
    }

    // TODO - these functions have been ported from the old template.
    //        they should go into a utilities module. Should also add checking of args.
    //        Add error checking to normalised - to check for zero magnitude vector

    /**
     * Filter function to return only unique values in an array
     * @param {any} value Value of current element
     * @param {number} index Index of current element
     * @param {any[]} arr Array to filter
     * @returns {boolean}
     */
    static only_unique(value, index, arr) {
        return arr.indexOf(value) === index;
    }

    /**
     * Object to return from the get_local_datum function
     * @typedef GetLocalDatumReturn
     * @property {number[]} origin Origin coordinate
     * @property {number[]} y_axis Y-axis vector
     */
    /** For each time step calculate the dot product of the head node displacement (just pick one for efficiency)
     * and the vector that defines the cross car direction (note this can change because the car moves during the analysis too)
     * assumption is that n1 and n2' define car x axis are the datum and Y offsett of n1 at t=0 (first state) is calculated
     * as the Y_n1 - Y_seat_centre = Y_offset. After state 1, Y_offset is used as the datum and excursion is measured as
     * Y distance perpendicular to n1 and n2' (note that n2' is coordinate of n2 rotated back to global X axis at t=0 )
     * @param {HeadExcursionStructure} head_excursion Head Excursion structure with data required for assessment
     * @param {number[][]} rotation_matrix Rotation matrix
     * @return {GetLocalDatumReturn}
     */
    static get_local_datum(head_excursion, rotation_matrix) {
        /* note n1 and n2 are external node ids so must be -ve */
        let n1_coords = GetData(CV, NODE, -head_excursion.shift_deform_nodes[0]);
        let n2_coords = GetData(CV, NODE, -head_excursion.shift_deform_nodes[1]);

        /** @type {number[]} */
        let origin_local = [n1_coords[0], n1_coords[1], n1_coords[2]];

        /* local coordinate */

        let v_n1n2 = ProtocolAssessment.get_normalised_vector_from_coords(
            [n1_coords[0], n1_coords[1], n1_coords[2]],
            [n2_coords[0], n2_coords[1], n2_coords[2]]
        );

        /* y_axis_local is used to take dot product with to find excursion distance (relative to origin_local) */

        let y_axis_local = ProtocolAssessment.vector_x_matrix(v_n1n2, rotation_matrix);

        return { origin: origin_local, y_axis: y_axis_local };
    }

    /**
     * Calculates the normalised vector from one coordinate to another.
     * Expects coordinates to be in the form [x, y, z]
     * @param {number[]} coords_1 First coordinates
     * @param {number[]} coords_2 Second coordinates
     * @returns {number[]}
     */
    static get_normalised_vector_from_coords(coords_1, coords_2) {
        let v_n1n2 = new Array(3);

        /* calculate n1->n2 vector */
        for (let i = 0; i < 3; i++) v_n1n2[i] = coords_2[i] - coords_1[i];

        let mag = ProtocolAssessment.magnitude(v_n1n2);

        /* normalise n1->n2 vector */
        for (let i = 0; i < 3; i++) v_n1n2[i] /= mag;

        return v_n1n2;
    }

    /**
     * Returns the current excursion from the seat centre line
     * @param {number} origin_y Origin Y coordinate
     * @param {GetLocalDatumReturn} local Local datum
     * @param {number} seat_centre_y Set centre Y coordinate
     * @param {number} head_internal_node_id Hed internal node id
     * @param {number} sign Sign -1 or +1
     * @returns {number}
     */
    static get_current_excursion_from_seat_cl(origin_y, local, seat_centre_y, head_internal_node_id, sign) {
        var head_node_coords = GetData(CV, NODE, head_internal_node_id);

        //vector from current n1 to current head node
        var origin_vector = [
            head_node_coords[0] - local.origin[0],
            head_node_coords[1] - local.origin[1],
            head_node_coords[2] - local.origin[2]
        ];

        var distance_of_head_from_origin_in_local_y = ProtocolAssessment.dot(local.y_axis, origin_vector);

        /*   Example of vector sum

                <<--------------------LOCAL Y-----------------

                Seat       Head   C/L    First shift deformed node
                350        100     0             -700
                x           H      |              N1
                            <_______________________   distance_of_head_from_origin_in_local_y
                            +
                                    _______________>   oGlblData.origin_y
                            -
                <__________________                    oGlblData.seat_centre_y       
                            *
                            1                          sign              
                    
                            ==
                ____________>        
                
                ((100-(-700)) + (-700) - (350)) * 1 = -250 => Yellow/Green Zone boundary


                <<--------------------LOCAL Y-----------------

                        C/L  Head   Seat     First shift deformed node (intentionally picked on barrier side to confirm it still works)
                            0    -100  -350     -600
                            |      H     x       N1
                                <______________    distance_of_head_from_origin_in_local_y
                            +
                            ______________________>   oGlblData.origin_y
                            -
                            _____________>            oGlblData.seat_centre_y       
                            *
                        -1                         sign              
                    
                            ==
                                ______>        
                    
                ((-100-(-600)) + (-600) - (-350)) * -1 = -250 => Yellow/Green Zone boundary


        */

        //note that return value is +ve if head excursion passes seat centreline and -ve if it does not
        //i.e. the smaller (more negative) the return value, the less the head moves towards the barrier side

        return (distance_of_head_from_origin_in_local_y + origin_y - seat_centre_y) * sign;
    }

    /**
     * Returns the excursion zone
     * @param {number} excursion Excursion value
     * @param {HeadExcursionStructure} head_excursion Head excursion
     * @param {number} unit_system
     * @returns {string}
     */
    static get_excursion_zone(excursion, head_excursion, unit_system) {
        //get excursion line reference
        //excursion is relative to seat centre y
        //a positive excursion value is worse than negative because it means head goes past seat centreline

        //convert to mm so that comparison is consistent

        excursion = excursion / WorkflowUnits.LengthToMillimetresFactor(unit_system);

        Message("get_excursion_zone => excursion in mm = " + excursion);

        let red_excursion_line_from_passenger_seat_centre_y = head_excursion.intrusion_from_seat_centre_y;
        let orange_excursion_line_from_passenger_seat_centre_y = 0;
        let yellow_excursion_line_from_passenger_seat_centre_y = -125;
        let green_excursion_line_from_passenger_seat_centre_y = -250;

        // protocol specifies that red band should be brown if the side intrusion is >125mm outboard from seat cl and the vehicle has suitable countermeasure

        //convert countermeasure string to boolean
        let countermeasure = false;
        if (/yes/i.test(head_excursion.countermeasure)) countermeasure = true;

        if (countermeasure && red_excursion_line_from_passenger_seat_centre_y > 125) {
            Message(
                "Red band is higher as it is >125mm from seat centreline and vehicle has countermeasures so colour zone brown"
            );
            // TODO - need to set this somehow
            //EuroNCAPColourREDorBROWN = EuroNCAPColourBROWN;
        }

        if (excursion > red_excursion_line_from_passenger_seat_centre_y) {
            return "CAPPING";
        } else if (excursion > orange_excursion_line_from_passenger_seat_centre_y) {
            // if > 0
            if (countermeasure && red_excursion_line_from_passenger_seat_centre_y > 125) {
                Message(
                    "Red band is higher as it is >125mm from seat centreline and vehicle has countermeasures so colour zone brown"
                );
                return "BROWN"; //"RED_GT125";
            }

            return "RED"; //countermeasure = false or RED_LT125
        } else if (excursion > yellow_excursion_line_from_passenger_seat_centre_y) {
            return "ORANGE";
        } else if (excursion > green_excursion_line_from_passenger_seat_centre_y) {
            return "YELLOW";
        } else {
            return "GREEN";
        }
    }

    //matrix operations:

    /**
     * get_rotation_matrix function is based on logic from the following stackexchange question:
     * https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/897677#897677
     *
     * the rotation matrix is used to rotate the n1->n2 vector used to define the shift deform
     * to the global y axis at the first state. The same rotation matrix is then used to calculate
     * the local y axis which is the assumed direction to measure excursion in.
     *
     * e.g. if the car rotates on impact it will not be correct/meaningful to calculate the y
     * displacement of the head nodes as the excursion will be in the 'cross-car' direction which
     * is what we call local y because the initial cross car direction is assumed to be the same
     * as global y (i.e. the car is pointing in the x direction)
     * @param {number[]} a
     * @param {number[]} b
     * @returns
     */
    static get_rotation_matrix(a, b) {
        // let a = [1, 1, 1];
        // let b = [Math.sqrt(0.5), Math.sqrt(0.5), 0];//[0,0,1];

        let v = ProtocolAssessment.cross(a, b);
        let s = ProtocolAssessment.magnitude(v);
        let c = ProtocolAssessment.dot(a, b);
        let vx = [
            [0, -v[2], v[1]],
            [v[2], 0, -v[0]],
            [-v[1], v[0], 0]
        ];
        let I = [
            [1, 0, 0],
            [0, 1, 0],
            [0, 0, 1]
        ];
        let mult = (1 - c) / (s * s);
        let w = [
            [mult, 0, 0],
            [0, mult, 0],
            [0, 0, mult]
        ];
        let r = ProtocolAssessment.matrix_3x3_Add([
            I,
            vx,
            ProtocolAssessment.multiply_matrices(w, ProtocolAssessment.multiply_matrices(vx, vx))
        ]);
        let rotated_vector = ProtocolAssessment.vector_x_matrix(a, r);

        Message("rotation_matrix  = " + JSON.stringify(r));
        Message("b_rotated  = " + JSON.stringify(rotated_vector));
        Message("b_original = " + JSON.stringify(b));

        return r;
    }

    /**
     * Returns the multiplication of a vector and matrix
     * @param {number[]} v
     * @param {number[][]} m
     * @returns {number[]}
     */
    static vector_x_matrix(v, m) {
        if (m[0][0] == null) {
            //if the matrix has null elements then original n1->n2 vector was parallel to global Y
            //so no rotation required and just return v
            Message(
                "Warning: First element in matrix is null which will occur if the choice of n1->n2 is parallel to global y. "
            );
            Message("The normalised n1->n2 vector will therefore be used so we return it (v).");
            return v;
        }
        let result = [0, 0, 0];
        for (let row = 0; row < 3; row++) {
            for (let i = 0; i < 3; i++) {
                result[row] += m[row][i] * v[i];
            }
        }

        return result;
    }

    /**
     * Return the matrix addition of an array of 3x3 matrices
     * @param {number[][][]} matrices
     * @returns {number[][]}
     */
    static matrix_3x3_Add(matrices) {
        let total_matrix = Array(3);

        for (let matrix = 0; matrix < matrices.length; matrix++) {
            for (let row = 0; row < 3; row++) {
                if (matrix == 0) total_matrix[row] = [0, 0, 0];
                for (let column = 0; column < 3; column++) {
                    total_matrix[row][column] += matrices[matrix][row][column];
                }
            }
        }

        return total_matrix;
    }

    /**
     * Return the matrix multiplication of two matrices
     * @param {number[][]} m1
     * @param {number[][]} m2
     * @returns {number[][]}
     */
    static multiply_matrices(m1, m2) {
        let result = [];
        for (let i = 0; i < m1.length; i++) {
            result[i] = [];
            for (let j = 0; j < m2[0].length; j++) {
                let sum = 0;
                for (let k = 0; k < m1[0].length; k++) {
                    sum += m1[i][k] * m2[k][j];
                }
                result[i][j] = sum;
            }
        }
        return result;
    }

    /**
     * Return the cross product of two vectors
     * @param {number[]} A
     * @param {number[]} B
     * @returns {number[]}
     */
    static cross(A, B) {
        let a1, a2, a3, b1, b2, b3;
        [a1, a2, a3] = A;
        [b1, b2, b3] = B;
        return [a2 * b3 - a3 * b2, a3 * b1 - a1 * b3, a1 * b2 - a2 * b1];
    }

    /**
     * Return the dot product of two vectors
     * @param {number[]} A
     * @param {number[]} B
     * @returns {number}
     */
    static dot(A, B) {
        let a1, a2, a3, b1, b2, b3;
        [a1, a2, a3] = A;
        [b1, b2, b3] = B;
        return a1 * b1 + a2 * b2 + a3 * b3;
    }

    /**
     * Return the magnitude of a vector
     * @param {number[]} A Vector
     * @returns {number}
     */
    static magnitude(A) {
        return Math.sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
    }
}

/**
 * Class to hold data for a structure its model ID, unit system and extra data
 * Used as an argument in the ProtocolAssessment.DoStructureAssessment function
 */
class DoStructureAssessmentStructureData {
    /**
     * @param {Structure} structure Structure instance
     * @param {number} model_id Model ID
     * @param {number} unit_system Unit system, e.g. Workflow.UNIT_SYSTEM_U2
     * @param {object} [extra=null] Extra data
     * @example
     * let sd = new DoStructureAssessmentStructureData(s, 1, Workflow.UNIT_SYSTEM_U2);
     */
    constructor(structure, model_id, unit_system, extra = null) {
        this.structure = structure;
        this.model_id = model_id;
        this.unit_system = unit_system;
        this.extra = extra;
    }

    /* Instance property getter and setters */

    /**
     * Structure instance
     * @type {Structure} */
    get structure() {
        return this._structure;
    }
    set structure(new_structure) {
        if (!(new_structure instanceof Structure)) {
            throw new Error(
                `<structure> is not a Structure instance in <DoStructureAssessmentStructureData> constructor`
            );
        }

        this._structure = new_structure;
    }
    /**
     * Model ID
     * @type {number} */
    get model_id() {
        return this._model_id;
    }
    set model_id(new_model_id) {
        if (typeof new_model_id != "number") {
            throw new Error(`<model_id> is not a number in <DoStructureAssessmentStructureData> constructor`);
        }

        this._model_id = new_model_id;
    }
    /**
     * Unit system, e.g. Workflow.UNIT_SYSTEM_U2
     * @type {number} */
    get unit_system() {
        return this._unit_system;
    }
    set unit_system(new_unit_system) {
        if (typeof new_unit_system != "number") {
            throw new Error(`<unit_system> is not a number in <DoStructureAssessmentStructureData> constructor`);
        }

        this._unit_system = new_unit_system;
    }
    /**
     * Object for extra data
     * @type {Object} */
    get extra() {
        return this._extra;
    }
    set extra(new_extra) {
        if (!(new_extra instanceof Object)) {
            throw new Error(`<extra> is not a Object instance in <DoStructureAssessmentStructureData> constructor`);
        }

        this._extra = new_extra;
    }
}

/**
 * Class used to return results from the ProtocolAssessment.DoStructureAssessment function
 */
class DoAssessmentResults {
    constructor() {
        this.output = {};
        this.image_filenames = [];
    }

    /* Instance property getter and setters */

    /**
     * Output object. This can store any results you want, but the properties are generally
     * set to something that describes the result, e.g "M1 NECK Driver-front-left max_extension"
     * @type {object} */
    get output() {
        return this._output;
    }
    set output(new_output) {
        if (typeof new_output != "object") {
            throw new Error(`<output> is not an object in <DoAssessmentResults>`);
        }

        this._output = new_output;
    }

    /**
     * Image filename
     * @type {string[]} */
    get image_filenames() {
        return this._image_filenames;
    }
    set image_filenames(new_image_filenames) {
        if (!(new_image_filenames instanceof Array)) {
            throw new Error(`<image_filenames> must be an array in <DoAssessmentResults>`);
        }

        for (let image_filename of new_image_filenames) {
            if (typeof image_filename != "string") {
                throw new Error(
                    "Not all of the items in the <image_filenames> array are strings in <DoAssessmentResults>"
                );
            }
        }

        this._image_filenames = new_image_filenames;
    }
}