egs_prism.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ prism geometry headers
#  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:    Reid Townson
#
###############################################################################
*/


#ifndef EGS_PRISM_
#define EGS_PRISM_

#include "egs_base_geometry.h"
#include "egs_polygon.h"
#include "egs_functions.h"

#ifdef WIN32

    #ifdef BUILD_PRISM_DLL
        #define EGS_PRISM_EXPORT __declspec(dllexport)
    #else
        #define EGS_PRISM_EXPORT __declspec(dllimport)
    #endif
    #define EGS_PRISM_LOCAL

#else

    #ifdef HAVE_VISIBILITY
        #define EGS_PRISM_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_PRISM_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_PRISM_EXPORT
        #define EGS_PRISM_LOCAL
    #endif

#endif

template <class T>
class EGS_PRISM_EXPORT EGS_PrismT : public EGS_BaseGeometry {

protected:

    T      *p;             
    EGS_Vector a;          
    EGS_Float
    d1, 
    d2; 

    bool   open; 

public:

    EGS_PrismT(T *P, const string &Name="") :
        EGS_BaseGeometry(Name), p(P), a(p->getNormal()), open(true) {
        is_convex = p->isConvex();
        nreg = 1;
    };

    EGS_PrismT(T *P, EGS_Float D1, EGS_Float D2, const string &Name="") :
        EGS_BaseGeometry(Name), p(P), a(p->getNormal()),
        d1(D1), d2(D2), open(false) {
        if (d1 > d2) {
            d1 = D2;
            d2 = D1;
        }
        is_convex = p->isConvex();
        nreg = 1;
    };

    ~EGS_PrismT() {
        delete p;
    };

    bool isInside(const EGS_Vector &x) {
        if (!open) {
            EGS_Float d = p->distance(x);
            if (d < d1 || d > d2) {
                return false;
            }
        }
        return p->isInside2D(x);
    };

    int isWhere(const EGS_Vector &x) {
        if (isInside(x)) {
            return 0;
        }
        else {
            return -1;
        }
    };
    int inside(const EGS_Vector &x) {
        return isWhere(x);
    };

    int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
               EGS_Float &t, int *newmed = 0, EGS_Vector *normal = 0) {
        EGS_2DVector v;
        EGS_2DVector *pv = normal ? &v : 0;
        if (open) {
            bool hit = p->howfar2D(ireg == 0 ? true : false,x,u,t,pv);
            if (!hit) {
                return ireg;
            }
            if (newmed) {
                *newmed = !ireg ? -1 : med;
            }
            if (normal) {
                *normal = p->getNormal(v);
            }
            return !ireg ? -1 : 0;
        }
        EGS_Float up = a*u, d = p->distance(x);
        if (!ireg) {  // inside
            EGS_Float tt = veryFar;
            int inew = ireg;
            if (up > boundaryTolerance) {
                tt = (d2 - d)/up;
            }
            else if (up < -boundaryTolerance) {
                tt = (d1 - d)/up;
            }
            else {
                tt = 0;
            }
            if (tt <= t) {
                if (tt > 0) {
                    t = tt+boundaryTolerance;
                }
                else {
                    t = tt;
                }
                inew = -1;
                if (normal) {
                    *normal = up>0 ? a*(-1) : a;
                }
                if (newmed) {
                    *newmed = -1;
                }
            }
            bool hit = p->howfar2D(true,x,u,t,pv);
            if (!hit) {
                return inew;
            }
            if (newmed) {
                *newmed = -1;
            }
            if (normal) {
                *normal = p->getNormal(v);
            }
            return -1;
        }
        if (d < d1 || d > d2) {
            EGS_Float tt = veryFar;
            if (d < d1 && up > boundaryTolerance) {
                tt = (d1 - d)/up;
            }
            else if (d > d2 && up < -boundaryTolerance) {
                tt = (d2 - d)/up;
            }
            if (tt < t) {
                EGS_Vector xp(x + u*tt);
                if (p->isInside2D(xp)) {
                    t = tt;
                    if (newmed) {
                        *newmed = med;
                    }
                    if (normal) {
                        *normal = up>0 ? a*(-1) : a;
                    }
                    return 0;
                }
            }
        }
        EGS_Float tt = t;
        bool hit = p->howfar2D(false,x,u,tt,pv);
        if (!hit) {
            return ireg;
        }
        d = p->distance(x+u*tt);
        if (d >= d1 && d <= d2) {
            t = tt;
            if (newmed) {
                *newmed = med;
            }
            if (normal) {
                *normal = p->getNormal(v);
            }
            return 0;
        }
        return ireg;
    };

    EGS_Float hownear(int ireg, const EGS_Vector &x) {
        EGS_Float tperp = p->hownear2D(ireg == 0 ? true : false,x);
        if (open) {
            return tperp;
        }
        EGS_Float d = p->distance(x);
        if (!ireg) {  // inside
            EGS_Float t = d2 - d;
            if (t < tperp) {
                tperp = t;
            }
            t = d - d1;
            if (t < tperp) {
                tperp = t;
            }
        }
        else {
            EGS_Float t;
            if (d < d1) {
                t = d1 - d;
            }
            else if (d > d2) {
                t = d - d2;
            }
            else {
                return tperp;
            }
            if (p->isInside2D(x)) {
                tperp = t;
            }
            else {
                tperp = sqrt(tperp*tperp + t*t);
            }
        }
        return tperp;
    };

    //const string &getType() const { return type; };
    const string &getType() const {
        return p->getType();
    };

    void printInfo() const {
        EGS_BaseGeometry::printInfo();
        if (open) {
            egsInformation("  open\n");
        }
        else {
            egsInformation("  closed with planes at %g and %g\n",d1,d2);
        }
        egsInformation("===================================================\n");
    };

};

typedef EGS_PrismT<EGS_PolygonYZ> EGS_PrismX;
typedef EGS_PrismT<EGS_PolygonXZ> EGS_PrismY;
typedef EGS_PrismT<EGS_PolygonXY> EGS_PrismZ;
typedef EGS_PrismT<EGS_Polygon> EGS_Prism;

#endif