egs_prism.cpp

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ prism geometry
#  Copyright (C) 2015 National Research Council Canada
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Iwan Kawrakow, 2005
#
#  Contributors:    Frederic Tessier
#                   Marc Chamberland
#                   Reid Townson
#                   Hannah Gallop
#
###############################################################################
*/
 
 
#include "egs_prism.h"
#include "egs_input.h"
 
#include <vector>
using std::vector;
 
static EGS_PRISM_LOCAL string __prismX("EGS_PrismX");
static EGS_PRISM_LOCAL string __prismY("EGS_PrismY");
static EGS_PRISM_LOCAL string __prismZ("EGS_PrismZ");
static EGS_PRISM_LOCAL string __prism("EGS_Prism");
 
static bool EGS_PRISM_LOCAL inputSet = false;
 
extern "C" {
 
    static void setInputs() {
        inputSet = true;
 
        setBaseGeometryInputs();
 
        geomBlockInput->getSingleInput("library")->setValues({"egs_prism"});
 
        // Format: name, isRequired, description, vector string of allowed values
        auto typePtr = geomBlockInput->addSingleInput("type", true, "The type of prism", {"EGS_PrismX", "EGS_PrismY", "EGS_PrismZ", "EGS_Prism"});
 
        geomBlockInput->addSingleInput("closed", false, "Two inputs that define the distance from the top and bottom prism plane to the plane used to define the polygon");
        geomBlockInput->addSingleInput("points", true, "A list of 2D or 3D (for type=EGS_Prism) positions. For 3D points, they must reside on a plane. These create a polygon that define the base of the prism.");
    }
 
    EGS_PRISM_EXPORT string getExample() {
        string example;
        example = {
            R"(
    # Example of egs_prisms
 
    # EGS_PrismZ example
    #:start geometry:
        name = my_prism
        library = egs_prism
        type = EGS_PrismZ
        points = 1 1  -1 1  -1 -1  4 -1
        closed = 1 4
        :start media input:
            media = air
        :stop media input:
    :stop geometry:
)"};
        return example;
    }
 
    EGS_PRISM_EXPORT shared_ptr<EGS_BlockInput> getInputs() {
        if(!inputSet) {
            setInputs();
        }
        return geomBlockInput;
    }
 
    EGS_PRISM_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
 
        if (!input) {
            egsWarning("createGeometry(prism): null input?\n");
            return 0;
        }
        string type;
        int err = input->getInput("type",type);
        if (err) {
            egsWarning("createGeometry(prism): missing 'type' input\n");
            return 0;
        }
        EGS_BaseGeometry *g;
        vector<EGS_Float> p;
        err = input->getInput("points",p);
        if (err) {
            egsWarning("createGeometry(prism): missing 'points' input\n");
            return 0;
        }
        bool open = true;
        vector<EGS_Float> tmp;
        err = input->getInput("closed",tmp);
        if (!err && tmp.size() == 2) {
            open = false;
        }
        bool p_open = false;
        int itmp;
        err = input->getInput("open triangle",itmp);
        if (!err && itmp == 1) {
            p_open = true;
        }
        if (input->compare(type,__prismX) || input->compare(type,__prismY)
                || input->compare(type,__prismZ)) {
            int np = p.size()/2;
            if (np < 3) {
                egsWarning("createGeometry(prism): at least 3 points are required "
                           "to construct a prism\n");
                return 0;
            }
            vector<EGS_2DVector> points;
            for (int j=0; j<np; j++) {
                points.push_back(EGS_2DVector(p[2*j],p[2*j+1]));
            }
            if (open) {
                if (input->compare(type,__prismX))
                    g=new EGS_PrismX(new EGS_PolygonYZ(points,
                                                       EGS_XProjector(__prismX),p_open));
                else if (input->compare(type,__prismY))
                    g=new EGS_PrismY(new EGS_PolygonXZ(points,
                                                       EGS_YProjector(__prismY),p_open));
                else
                    g=new EGS_PrismZ(new EGS_PolygonXY(points,
                                                       EGS_ZProjector(__prismZ),p_open));
            }
            else {
                if (input->compare(type,__prismX))
                    g=new EGS_PrismX(new EGS_PolygonYZ(points,
                                                       EGS_XProjector(__prismX),p_open),tmp[0],tmp[1]);
                else if (input->compare(type,__prismY))
                    g=new EGS_PrismY(new EGS_PolygonXZ(points,
                                                       EGS_YProjector(__prismY),p_open),tmp[0],tmp[1]);
                else
                    g=new EGS_PrismZ(new EGS_PolygonXY(points,
                                                       EGS_ZProjector(__prismZ),p_open),tmp[0],tmp[1]);
            }
        }
        else {
            int np = p.size()/3;
            if (np < 3) {
                egsWarning("createGeometry(prism): at least 3 points are required"
                           " to construct a prism\n");
                return 0;
            }
            vector<EGS_Vector> points;
            for (int j=0; j<np; j++) {
                points.push_back(EGS_Vector(p[3*j],p[3*j+1],p[3*j+2]));
            }
            EGS_Vector aux(points[np-1]-points[0]);
            if (aux.length2() > epsilon) {
                points.push_back(points[0]);
            }
            np = points.size();
            EGS_Projector pro(points[0],points[1],points[np-2],__prism);
            vector<EGS_2DVector> p2;
            {
                for (int j=0; j<np; j++) {
                    p2.push_back(pro.getProjection(points[j]));
                    EGS_Float d = pro.distance(points[j]);
                    if (fabs(d) > epsilon) {
                        egsWarning("createGeometry(prism): "
                                   "points are not on a plane\n");
                        return 0;
                    }
                }
            }
            if (open) {
                g = new EGS_Prism(new EGS_Polygon(p2,pro,p_open));
            }
            else {
                g = new EGS_Prism(new EGS_Polygon(p2,pro,p_open),tmp[0],tmp[1]);
            }
        }
        g->setName(input);
        g->setBoundaryTolerance(input);
        g->setMedia(input);
        g->setLabels(input);
        return g;
    }
 
}