egs_shapes.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ shapes 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:    Frederic Tessier
#
###############################################################################
*/


#ifndef EGS_SHAPES_
#define EGS_SHAPES_

#include "egs_vector.h"
#include "egs_transformations.h"
#include "egs_rndm.h"
#include "egs_object_factory.h"

#include <string>
using std::string;

class EGS_Input;

class EGS_EXPORT EGS_BaseShape : public EGS_Object {

public:

    EGS_BaseShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_Object(Name,f), T(0) {
        otype = "base_shape";
    };
    virtual ~EGS_BaseShape() {
        if (T) {
            delete T;
        }
    };

    virtual EGS_Vector getRandomPoint(EGS_RandomGenerator *rndm) {
        if (T) {
            return (*T)*getPoint(rndm);
        }
        else {
            return getPoint(rndm);
        }
    };

    virtual EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        egsFatal("You need to implement the getPoint function in your "
                 "derived class\n");
        return EGS_Vector();
    };

    void setTransformation(EGS_Input *inp);

    void setTransformation(EGS_AffineTransform *t) {
        if (T) {
            delete T;
        }
        T = new EGS_AffineTransform(*t);
    };

    const EGS_AffineTransform *getTransform() const {
        return T;
    };

    static EGS_BaseShape *createShape(EGS_Input *inp);

    static EGS_BaseShape *getShape(const string &Name);

    virtual bool supportsDirectionMethod() const {
        return false;
    };

    virtual void getPointSourceDirection(const EGS_Vector &xo,
                                         EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        egsFatal("getPointSourceDirection: you have to implement this "
                 "method for the %s shape if you want to use it\n",otype.c_str());
    };

    virtual EGS_Float area() const {
        return 1;
    };

protected:

    EGS_AffineTransform *T; 

};

class  EGS_EXPORT EGS_SurfaceShape : public EGS_BaseShape {

public:

    EGS_SurfaceShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), A(1) {};
    ~EGS_SurfaceShape() {};
    bool supportsDirectionMethod() const {
        return true;
    };
    EGS_Float area() const {
        return A;
    };
    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Vector x = getPoint(rndm);
        u = x - xo;
        EGS_Float d2i = 1/u.length2(), di = sqrt(d2i);
        u *= di;
        wt = A*fabs(u.z)*d2i;
        if (T) {
            T->rotate(u);
        }
    };

protected:

    EGS_Float A;

};

class EGS_EXPORT EGS_PointShape : public EGS_BaseShape {

public:

    EGS_PointShape(const EGS_Vector &Xo = EGS_Vector(),
                   const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), xo(Xo) {
        otype = "point";
    };
    ~EGS_PointShape() { };
    EGS_Vector getPoint(EGS_RandomGenerator *) {
        return xo;
    };
    EGS_Object *createObject(EGS_Input *inp);

protected:

    EGS_Vector  xo; 

};

class EGS_EXPORT EGS_BoxShape : public EGS_BaseShape {

protected:

    EGS_Float  ax, ay, az; 

public:

    EGS_BoxShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), ax(1), ay(1), az(1) {
        otype="box";
    };
    EGS_BoxShape(EGS_Float A, const EGS_AffineTransform *t = 0,
                 const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), ax(A), ay(A), az(A) {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        otype="box";
    };
    EGS_BoxShape(EGS_Float Ax, EGS_Float Ay, EGS_Float Az,
                 const EGS_AffineTransform *t = 0,
                 const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), ax(Ax), ay(Ay), az(Az) {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        otype="box";
    };
    ~EGS_BoxShape() { };

    EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        EGS_Vector v(ax*(rndm->getUniform()-0.5),
                     ay*(rndm->getUniform()-0.5),
                     az*(rndm->getUniform()-0.5));
        return v;
    };

    EGS_Object *createObject(EGS_Input *);

    bool supportsDirectionMethod() const {
        return true;
    };

    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Float eta = rndm->getUniform()*area();
        if (eta < 2*ax*ay) {
            u.x = ax*(rndm->getUniform()-0.5);
            u.y = ay*(rndm->getUniform()-0.5);
            if (eta < ax*ay) {
                u.z = az/2;
                wt = u.z - xo.z;
            }
            else {
                u.z = -az/2;
                wt = xo.z - u.z;
            }
        }
        else if (eta < 2*(ax*ay + ax*az)) {
            u.x = ax*(rndm->getUniform()-0.5);
            u.z = az*(rndm->getUniform()-0.5);
            if (eta < 2*ax*ay + ax*az) {
                u.y = ay/2;
                wt = u.y - xo.y;
            }
            else {
                u.y = -ay/2;
                wt = xo.y - u.y;
            }
        }
        else {
            eta -= 2*(ax*ay + ax*az);
            u.y = ay*(rndm->getUniform()-0.5);
            u.z = az*(rndm->getUniform()-0.5);
            if (eta < ay*az) {
                u.x = ax/2;
                wt = u.x - xo.x;
            }
            else {
                u.x = -ax/2;
                wt = xo.x - u.x;
            }
        }
        u -= xo;
        EGS_Float d2 = u.length2(), d = sqrt(d2);
        u *= (1/d);
        wt *= (area()/(d2*d));
        if (T) {
            T->rotate(u);
        }
    };
    EGS_Float area() const {
        return 2*(ax*ay + ax*az + ay*az);
    };

};

class EGS_EXPORT EGS_SphereShape : public EGS_BaseShape {

protected:

    EGS_Float  R;     
    EGS_Vector xo;    

public:

    EGS_SphereShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(1), xo() {
        otype="sphere";
    };
    EGS_SphereShape(EGS_Float r, const EGS_Vector &Xo = EGS_Vector(0,0,0),
                    const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(r), xo(Xo) {
        otype = "sphere";
    };
    ~EGS_SphereShape() {};

    EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        EGS_Float r = rndm->getUniform(), r1 = rndm->getUniform(),
                  r2 = rndm->getUniform();
        if (r1 > r) {
            r = r1;
        }
        if (r2 > r) {
            r = r2;
        }
        EGS_Float cost = 2*rndm->getUniform()-1;
        EGS_Float sint = sqrt(1-cost*cost);
        r1 = R*r*sint;
        EGS_Float cphi, sphi;
        rndm->getAzimuth(cphi,sphi);
        return xo + EGS_Vector(r1*cphi,r1*sphi,R*r*cost);
    };

    EGS_Object *createObject(EGS_Input *inp);

    bool supportsDirectionMethod() const {
        return true;
    };

    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Float cost = 2*rndm->getUniform()-1;
        EGS_Float sint = 1-cost*cost;
        EGS_Vector x;
        if (sint > epsilon) {
            EGS_Float cphi, sphi;
            rndm->getAzimuth(cphi,sphi);
            sint = R*sqrt(sint);
            x.x = sint*cphi;
            x.y = sint*sphi;
            x.z = R*cost;
        }
        else {
            x.z = R*cost;
        }
        u = (x + this->xo) - xo;
        EGS_Float di = 1/u.length();
        u *= di;
        wt = u*x*4*M_PI*R*di*di;
    };
    EGS_Float area() const {
        return 4*M_PI*R*R;
    };
};

class EGS_EXPORT EGS_CylinderShape : public EGS_BaseShape {

protected:

    EGS_Float  R;       
    EGS_Float  h;       
    EGS_Vector xo;      
    EGS_Vector a;       
    EGS_Float  phi_min;
    EGS_Float  phi_max;
    bool       has_phi; 

    inline void getPointInCircle(EGS_RandomGenerator *rndm, EGS_Float &x,
                                 EGS_Float &y) {
        if (!has_phi) {
            do {
                x = 2*rndm->getUniform()-1;
                y = 2*rndm->getUniform()-1;
            }
            while (x*x + y*y > 1);
            x *= R;
            y *= R;
        }
        else {
            EGS_Float r = R*sqrt(rndm->getUniform());
            EGS_Float eta = rndm->getUniform();
            EGS_Float phi = phi_min*(1-eta) + phi_max*eta;
            x = r*cos(phi);
            y = r*sin(phi);
        }
    };


public:

    EGS_CylinderShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(), R(1), h(1), xo(), a(0,0,1),
        phi_min(0), phi_max(2*M_PI), has_phi(false) {
        otype="cylinder";
    };
    EGS_CylinderShape(EGS_Float r, EGS_Float H,
                      const EGS_Vector &Xo = EGS_Vector(0,0,0),
                      const EGS_Vector &A = EGS_Vector(0,0,1),
                      const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(r), h(H), xo(Xo), a(A),
        phi_min(0), phi_max(2*M_PI), has_phi(false) {
        EGS_RotationMatrix rmat(a);
        if (xo.length2() > epsilon || !rmat.isI()) {
            T = new EGS_AffineTransform(rmat.inverse(),xo);
        }
        otype="cylinder";
    };
    EGS_CylinderShape(EGS_Float r, EGS_Float H, const EGS_AffineTransform *t,
                      const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(r), h(H),
        phi_min(0), phi_max(2*M_PI), has_phi(false) {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        otype="cylinder";
    };
    ~EGS_CylinderShape() { };

    void setPhiRange(EGS_Float Phi_min, EGS_Float Phi_max) {
        if (Phi_min < Phi_max) {
            phi_min = Phi_min;
            phi_max = Phi_max;
        }
        else                    {
            phi_min = Phi_max;
            phi_max = Phi_min;
        }
        if (phi_max - phi_min < 1.99999*M_PI) {
            has_phi = true;
        }
        else {
            has_phi = false;
        }
    };

    EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        EGS_Float x,y;
        getPointInCircle(rndm,x,y);
        EGS_Float z = h*(rndm->getUniform()-0.5);
        return EGS_Vector(x,y,z);
    };

    EGS_Object *createObject(EGS_Input *);

    EGS_Float getRadius() const {
        return R;
    };
    EGS_Float getHeight() const {
        return h;
    };

    bool supportsDirectionMethod() const {
        return true;
    };

    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Float eta = rndm->getUniform()*(R+h);
        EGS_Vector x;                               // point on cylinder with respect to midpoint
        EGS_Vector n;                               // normal to cylinder at point x
        if (eta < R) {
            getPointInCircle(rndm,x.x,x.y);
            if (2*eta < R) {
                x.z = h/2;    // top face: normal is up
                n.z = 1;
            }
            else {
                x.z = -h/2;    // bottom face: normal is down
                n.z = -1;
            }
        }
        else {
            EGS_Float cphi,sphi;
            rndm->getAzimuth(cphi,sphi);
            x.x = R*cphi;
            x.y = R*sphi;
            x.z = h*(rndm->getUniform()-0.5);
            n.x = x.x;
            n.y = x.y;                   // side face: normal is (x,y)
        }
        u = (x+this->xo) - xo;                      // direction vector from origin to cylinder point
        EGS_Float d2 = u.length2(), d = sqrt(d2);
        u *= (1/d);                                 // normalize direction vectors
        n.normalize();                              // normalize normal
        wt = u*n*area()/d2;
        if (T) {
            T->rotate(u);
        }
    };
    EGS_Float area() const {
        return 2*M_PI*R*(R+h);
    };
};

#endif