geometry/egs_box/egs_box.cpp

The following describes the implementation of a box geometry that can serve as a guide to users for developing their own geometry classes.

Our box geometry class will describe a parallelepiped of a user defined size that defines a single region filled with 1 medium.

As our geometry must be derived from EGS_BaseGeometry, we include its header file. We also want to be able to transform the box from being centered about the origin and having its faces parallel to the x-, y- and z-axis to any other location and/or orientation. For that we will employ affine transformations and therefore also include the transformations header file:

#include "egs_base_geometry.h"
#include "egs_transformations.h"

As with other geometries, we will be compiling our box class into a shared library. We have to therefore mark the symbols that will be exported from the library. For this purpose we define a macro EGS_BOX_EXPORT, which we define like this

#ifdef WIN32

    #ifdef BUILD_BOX_DLL
        #define EGS_BOX_EXPORT __declspec(dllexport)
    #else
        #define EGS_BOX_EXPORT __declspec(dllimport)
    #endif
    #define EGS_BOX_LOCAL

#else

    #ifdef HAVE_VISIBILITY
        #define EGS_BOX_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_BOX_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_BOX_EXPORT
        #define EGS_BOX_LOCAL
    #endif

#endif

We are now ready to develop our box geometry class by deriving from EGS_BaseGeometry:

class EGS_BOX_EXPORT EGS_Box : public EGS_BaseGeometry {

    EGS_Float           ax, ay, az;
    EGS_AffineTransform *T;
    static string       type;

public:

The data members ax,ay,az are for the box size, the type static data member will be set to "EGS_Box" and used and the T data member will point to the affine transformation associated with boxes not located about the origin and/or rotated boxes. We now need a constructor and a destructor for box objects. We provide two different constructors: one for boxes where all sides have the same size (i.e. cubes)

    EGS_Box(EGS_Float a, const EGS_AffineTransform *t = 0,
            const string &Name = "") : EGS_BaseGeometry(Name),
        ax(a), ay(a), az(a), T(0)  {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        nreg = 1;
    };

and one for boxes with 3 different sides

    EGS_Box(EGS_Float Ax, EGS_Float Ay, EGS_Float Az,
            const EGS_AffineTransform *t = 0,
            const string &Name = "") : EGS_BaseGeometry(Name),
        ax(Ax), ay(Ay), az(Az), T(0)  {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        nreg = 1;
    };

In the destructor we simply deallocate the memory used by the affine transformation, if there was a transformation defined:

    ~EGS_Box() {
        if (T) {
            delete T;
        }
    };

We can now start implementing the required geometry methods isInside() , isWhere() , howfar() and hownear() .

A point is inside a box about the origin if its x-position is between -ax/2 and ax/2, its y-position between -ay/2 and ay/2, etc. For a transformed box, the easiest way to determine if the position is inside the box is to apply the inverse transformation to the position and to then use the simple method for boxes about the origin. The isInside implementation therefore looks like this

    bool isInside(const EGS_Vector &x) {
        EGS_Vector xp = T ? x*(*T) : x;
        if (2*xp.x + ax < 0 || 2*xp.x - ax > 0) {
            return false;
        }
        if (2*xp.y + ay < 0 || 2*xp.y - ay > 0) {
            return false;
        }
        if (2*xp.z + az < 0 || 2*xp.z - az > 0) {
            return false;
        }
        return true;
    };

(multiplying the position vector with *T from the right is an operation equivalent to transforming x with the inverse of *T)

The implementation of the isWhere method is very easy: as we have a single region this function should return 0 if the position is inside and -1 if it is outside:

    int isWhere(const EGS_Vector &x) {
        return isInside(x) ? 0 : -1;
    };

The implementation of the howfar method is also easiest if we first transform the position and direction to a frame where the box is centered about the origin and its sides are parallel to the axis:

The howfar logic is different for positions inside (i.e. ireg=0 ) and outside. For inside points, we must calculate the distances to the x-planes, y-planes and z-planes from the transformed point xp along the transformed direction up and take the minimum t1 of the three. If t1 is less than the intended step t, we must set t to t1 and return -1 (i.e. the particle will exit the box at the end of the step). If normal is not null, our howfar function is being called from the visualization utility and we must also calculate the normal of the plane being intersected to exit the box with the normal pointing towards the position of the particle. The implementation therefore looks like this

When the position is outside, the logic is

1. if the particle intersects the top or bottom side, check if the intersection point is between the left-right and east-west sides.
2. If yes, this is the position where the trajectory is entering the box => return the required information
3. If no, repeat 1-2 for left-right sides and if necessary for the east-west sides.
4. If all of the above fails, the particle does not enter the box.

The code is lengthy and not reproduced here.

The next step is the implementation of the hownear function. As with the other methods, we first inverse transform the position and then calculate the minimum distance to a box about the origin

If the position is inside the boxe (i.e. ireg=0), the minimum distance to a boundary is the minimum distance to any of the 6 faces:

If the position is outside, the algorithm is slightly more difficult and is left to the reader to understand from the implementation:

We reimplement the getType function to return type so that the type of our geometry class is known:

We also reimplement the printInfo() function so that users of our geometry class can get a better description of the geometry object:

The actual implementation of printInfo is in the egs_box.cpp file, which includes egs_box.h and egs_input.h:

/*
###############################################################################
#
#  EGSnrc egs++ box 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
#
###############################################################################
*/


#include "egs_box.h"
#include "egs_input.h"

The latter is necessary to gain access to the definition of the EGS_Input class that is needed in the createGeometry function (see below). The printInfo implementation is simple:

void EGS_Box::printInfo() const {
    EGS_BaseGeometry::printInfo();
    egsInformation(" box size = %g %g %g\n",ax,ay,az);
    egsInformation("=======================================================\n");
}

We also need to declare the static data member type

string EGS_Box::type("EGS_Box");

According to verious sources, compilers can produce better code if text strings are static and local to the library. Hence, we define the various error messages and the key we will be using to define the size of a box object as such:

static char EGS_BOX_LOCAL ebox_message1[] = "createGeometry(box): %s\n";
static char EGS_BOX_LOCAL ebox_message2[] = "null input?";
static char EGS_BOX_LOCAL ebox_message3[] = "wrong/missing 'box size' input?";
static char EGS_BOX_LOCAL ebox_message4[] =
    "expecting 1 or 3 float inputs for 'box size'";
static char EGS_BOX_LOCAL ebox_key1[] = "box size";

The remaining task now is to provide a C-style function named createGeometry for the box shared library that will be used to dynamically create box objects based on the information stored in an EGS_Input object.

extern "C" {

    EGS_BOX_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {

Just to be sure, we check that input is not null and issue a warning and return null if it is:

        if (!input) {
            egsWarning(ebox_message1,ebox_message2);

We then check if the input object contains a box size property

            return 0;
        }
        vector<EGS_Float> s;
        int err = input->getInput(ebox_key1,s);

and if it doesn't (indicated by err not being zero), issue a warning and return null:

        if (err) {
            egsWarning(ebox_message1,ebox_message3);
            return 0;
        }

Next we check if the input contains a definition if a transformation for our box geometry:

        EGS_AffineTransform *t = EGS_AffineTransform::getTransformation(input);

Depending on whether there was one or three values to the box size key, we use the corresponding constructor to create the new box object, but refuse to create it if there was some other number of values (because, e.g., the user made a mistake in the input file):

        EGS_Box *result;
        if (s.size() == 1) {
            result = new EGS_Box(s[0],t);
        }

We must now delete the affine transformation pointed to by t as the box object has made its own copy in the constructor

        else if (s.size() == 3) {
            result = new EGS_Box(s[0],s[1],s[2],t);
        }
        else {
            egsWarning(ebox_message1,ebox_message4);
            if (t) {
                delete t;

We then use the setName and setMedia functions inherited from EGS_BaseGeometry to set the name of the newly created box and fill it with a medium

            }
            return 0;
        }
        if (t) {
            delete t;
        }
        result->setName(input);

to then return the box

        result->setBoundaryTolerance(input);
        result->setMedia(input);
        result->setLabels(input);
        return result;
    }

}

Based on this createGeometry implementation, users will be able to use our box geometry class by having the following input in the input file

:start geometry:
    name = some_name
    box size = 3.5 1 7
    :start transformation:
        input defining the transformation
    :stop transformation
    :start media input:
        media = H2O
    :stop media input
:stop geometry:

The transformation related input can be of course missing in which case the box will be centered about the origin and has its faces parallel to the x-, y- and z-axis.

If we would have been nice to the users of our box geometry class, we could have provided an easier way for defining the medium of the box by re-implementing the setMedia() function, e.g. something along these lines:

void EGS_Box::setMedia(EGS_Input *inp) {
    string medname; int err = inp->getInput("medium",medname);
    if( !err ) setMedium(medname);
}

That way, the user could define the medium of the box by simply using

medium = H2O

instead of :start media input: ... :stop media input: which is necessary for geometries with several regions and potentialy different media and therefore provided as default implementation.

The final step is to provide a Makefile for our geometry DSO. We include the config files containing various definitions for the make system:

include $(EGS_CONFIG)
include $(SPEC_DIR)egspp.spec
include $(SPEC_DIR)egspp_$(my_machine).conf

We then must add the BUILD_BOX_DLL macro to the list of defined macros,

DEFS = $(DEF1) -DBUILD_BOX_DLL

set the name of the DSO,

library = egs_box

define the source files needed to create the DSO,

lib_files = egs_box

and tell the make system that our class depends on the egs_transformations.h header file,

my_deps = egs_transformations.h
extra_dep = $(addprefix $(DSOLIBS), $(my_deps))

We can then use the generic rules for building the DSO:

include $(SPEC_DIR)egspp_libs.spec

and because our only object file needed to build the library depends on egs_box.cpp and egs_box.h, we can also use the generic dependence rule

$(make_depend)

That's all.

Here is the complete source code of the EGS_Box class.
The header file:

/*
###############################################################################
#
#  EGSnrc egs++ box 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_BOX_
#define EGS_BOX_

#include "egs_base_geometry.h"
#include "egs_transformations.h"


#ifdef WIN32

    #ifdef BUILD_BOX_DLL
        #define EGS_BOX_EXPORT __declspec(dllexport)
    #else
        #define EGS_BOX_EXPORT __declspec(dllimport)
    #endif
    #define EGS_BOX_LOCAL

#else

    #ifdef HAVE_VISIBILITY
        #define EGS_BOX_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_BOX_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_BOX_EXPORT
        #define EGS_BOX_LOCAL
    #endif

#endif

class EGS_BOX_EXPORT EGS_Box : public EGS_BaseGeometry {

    EGS_Float           ax, ay, az;
    EGS_AffineTransform *T;
    static string       type;

public:

    EGS_Box(EGS_Float a, const EGS_AffineTransform *t = 0,
            const string &Name = "") : EGS_BaseGeometry(Name),
        ax(a), ay(a), az(a), T(0)  {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        nreg = 1;
    };

    EGS_Box(EGS_Float Ax, EGS_Float Ay, EGS_Float Az,
            const EGS_AffineTransform *t = 0,
            const string &Name = "") : EGS_BaseGeometry(Name),
        ax(Ax), ay(Ay), az(Az), T(0)  {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        nreg = 1;
    };

    ~EGS_Box() {
        if (T) {
            delete T;
        }
    };

    bool isInside(const EGS_Vector &x) {
        EGS_Vector xp = T ? x*(*T) : x;
        if (2*xp.x + ax < 0 || 2*xp.x - ax > 0) {
            return false;
        }
        if (2*xp.y + ay < 0 || 2*xp.y - ay > 0) {
            return false;
        }
        if (2*xp.z + az < 0 || 2*xp.z - az > 0) {
            return false;
        }
        return true;
    };

    int isWhere(const EGS_Vector &x) {
        return isInside(x) ? 0 : -1;
    };

    int inside(const EGS_Vector &x) {
        return isWhere(x);
    };

    EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
                              const EGS_Vector &u) {
        EGS_Float t = veryFar;
        howfar(ireg,x,u,t);
        return t;
    };

    int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
               EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
        EGS_Vector xp(x), up(u);
        if (T) {
            xp = x*(*T);
            up = u*T->getRotation();
        }
        if (ireg == 0) {
            EGS_Float t1 = veryFar;
            int inew = 0;
            EGS_Vector n;
            if (up.x > 0) {
                t1 = (ax - 2*xp.x)/(2*up.x);
                n.x = -1;
            }
            else if (up.x < 0) {
                t1 = -(ax + 2*xp.x)/(2*up.x);
                n.x = 1;
            }
            if (t1 < t) {
                t = t1;
                inew = -1;
            }
            t1 = veryFar;
            if (up.y > 0) {
                t1 = (ay - 2*xp.y)/(2*up.y);
            }
            else if (up.y < 0) {
                t1 = -(ay + 2*xp.y)/(2*up.y);
            }
            if (t1 < t) {
                n.x = 0;
                n.y = up.y > 0 ? -1 : 1;
                t = t1;
                inew = -1;
            }
            t1 = veryFar;
            if (up.z > 0) {
                t1 = (az - 2*xp.z)/(2*up.z);
            }
            else if (up.z < 0) {
                t1 = -(az + 2*xp.z)/(2*up.z);
            }
            if (t1 < t) {
                t = t1;
                inew = -1;
                n.x = n.y = 0;
                n.z = up.z > 0 ? -1 : 1;
            }
            if (inew < 0) {
                if (newmed) {
                    *newmed = -1;
                }
                if (normal) {
                    if (T) {
                        *normal = T->getRotation()*n;
                    }
                    else {
                        *normal = n;
                    }
                }
            }
            return inew;
        }
        EGS_Float t1 = veryFar;
        if (2*xp.x + ax < 0 && up.x > 0) {
            t1 = -(2*xp.x + ax)/(2*up.x);
        }
        else if (2*xp.x - ax > 0 && up.x < 0) {
            t1 = -(2*xp.x - ax)/(2*up.x);
        }
        if (t1 < t) {
            EGS_Float y1 = xp.y + up.y*t1, z1 = xp.z + up.z*t1;
            if (2*y1 + ay > 0 && 2*y1 - ay < 0 &&
                    2*z1 + az > 0 && 2*z1 - az < 0) {
                t = t1;
                if (newmed) {
                    *newmed = med;
                }
                if (normal) {
                    if (T) {
                        const EGS_RotationMatrix &R = T->getRotation();
                        if (up.x > 0) {
                            *normal = R*EGS_Vector(-1,0,0);
                        }
                        else {
                            *normal = R*EGS_Vector(1,0,0);
                        }
                    }
                    else {
                        if (up.x > 0) {
                            *normal = EGS_Vector(-1,0,0);
                        }
                        else {
                            *normal = EGS_Vector(1,0,0);
                        }
                    }
                }
                return 0;
            }
        }
        t1 = veryFar;
        if (2*xp.y + ay < 0 && up.y > 0) {
            t1 = -(2*xp.y + ay)/(2*up.y);
        }
        else if (2*xp.y - ay > 0 && up.y < 0) {
            t1 = -(2*xp.y - ay)/(2*up.y);
        }
        if (t1 < t) {
            EGS_Float x1 = xp.x + up.x*t1, z1 = xp.z + up.z*t1;
            if (2*x1 + ax > 0 && 2*x1 - ax < 0 &&
                    2*z1 + az > 0 && 2*z1 - az < 0) {
                t = t1;
                if (newmed) {
                    *newmed = med;
                }
                if (normal) {
                    if (T) {
                        const EGS_RotationMatrix &R = T->getRotation();
                        if (up.y > 0) {
                            *normal = R*EGS_Vector(0,-1,0);
                        }
                        else {
                            *normal = R*EGS_Vector(0,1,0);
                        }
                    }
                    else {
                        if (up.y > 0) {
                            *normal = EGS_Vector(0,-1,0);
                        }
                        else {
                            *normal = EGS_Vector(0,1,0);
                        }
                    }
                }
                return 0;
            }
        }
        t1 = veryFar;
        if (2*xp.z + az < 0 && up.z > 0) {
            t1 = -(2*xp.z + az)/(2*up.z);
        }
        else if (2*xp.z - az > 0 && up.z < 0) {
            t1 = -(2*xp.z - az)/(2*up.z);
        }
        if (t1 < t) {
            EGS_Float x1 = xp.x + up.x*t1, y1 = xp.y + up.y*t1;
            if (2*x1 + ax > 0 && 2*x1 - ax < 0 &&
                    2*y1 + ay > 0 && 2*y1 - ay < 0) {
                t = t1;
                if (newmed) {
                    *newmed = med;
                }
                if (normal) {
                    if (T) {
                        const EGS_RotationMatrix &R = T->getRotation();
                        if (up.z > 0) {
                            *normal = R*EGS_Vector(0,0,-1);
                        }
                        else {
                            *normal = R*EGS_Vector(0,0,1);
                        }
                    }
                    else {
                        if (up.z > 0) {
                            *normal = EGS_Vector(0,0,-1);
                        }
                        else {
                            *normal = EGS_Vector(0,0,1);
                        }
                    }
                }
                return 0;
            }
        }
        return -1;
    };

    EGS_Float hownear(int ireg, const EGS_Vector &x) {
        EGS_Vector xp = T ? x*(*T) : x;
        if (ireg >= 0) {
            EGS_Float tmin;
            EGS_Float t1 = xp.x + 0.5*ax, t2 = 0.5*ax - xp.x;
            if (t1 < t2) {
                tmin = t1;
            }
            else {
                tmin = t2;
            }
            t1 = xp.y + 0.5*ay;
            if (t1 < tmin) {
                tmin = t1;
            }
            t1 = 0.5*ay - xp.y;
            if (t1 < tmin) {
                tmin = t1;
            }
            t1 = xp.z + 0.5*az;
            if (t1 < tmin) {
                tmin = t1;
            }
            t1 = 0.5*az - xp.z;
            if (t1 < tmin) {
                tmin = t1;
            }
            return tmin;
        }
        EGS_Float s1=0, s2=0;
        int nout = 0;
        if (2*xp.x + ax < 0) {
            EGS_Float t = -0.5*ax - xp.x;
            s1 += t;
            s2 += t*t;
            nout++;
        }
        else if (2*xp.x - ax > 0) {
            EGS_Float t = xp.x - 0.5*ax;
            s1 += t;
            s2 += t*t;
            nout++;
        }
        if (2*xp.y + ay < 0) {
            EGS_Float t = -0.5*ay - xp.y;
            s1 += t;
            s2 += t*t;
            nout++;
        }
        else if (2*xp.y - ay > 0) {
            EGS_Float t = xp.y - 0.5*ay;
            s1 += t;
            s2 += t*t;
            nout++;
        }
        if (2*xp.z + az < 0) {
            EGS_Float t = -0.5*az - xp.z;
            s1 += t;
            s2 += t*t;
            nout++;
        }
        else if (2*xp.z - az > 0) {
            EGS_Float t = xp.z - 0.5*az;
            s1 += t;
            s2 += t*t;
            nout++;
        }
        if (nout < 2) {
            return s1;
        }
        return sqrt(s2);
    };

    const string &getType() const {
        return type;
    };

    void printInfo() const;

};

#endif

The .cpp file:

/*
###############################################################################
#
#  EGSnrc egs++ box 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
#
###############################################################################
*/


#include "egs_box.h"
#include "egs_input.h"

void EGS_Box::printInfo() const {
    EGS_BaseGeometry::printInfo();
    egsInformation(" box size = %g %g %g\n",ax,ay,az);
    egsInformation("=======================================================\n");
}

string EGS_Box::type("EGS_Box");

static char EGS_BOX_LOCAL ebox_message1[] = "createGeometry(box): %s\n";
static char EGS_BOX_LOCAL ebox_message2[] = "null input?";
static char EGS_BOX_LOCAL ebox_message3[] = "wrong/missing 'box size' input?";
static char EGS_BOX_LOCAL ebox_message4[] =
    "expecting 1 or 3 float inputs for 'box size'";
static char EGS_BOX_LOCAL ebox_key1[] = "box size";

extern "C" {

    EGS_BOX_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
        if (!input) {
            egsWarning(ebox_message1,ebox_message2);
            return 0;
        }
        vector<EGS_Float> s;
        int err = input->getInput(ebox_key1,s);
        if (err) {
            egsWarning(ebox_message1,ebox_message3);
            return 0;
        }
        EGS_AffineTransform *t = EGS_AffineTransform::getTransformation(input);
        EGS_Box *result;
        if (s.size() == 1) {
            result = new EGS_Box(s[0],t);
        }
        else if (s.size() == 3) {
            result = new EGS_Box(s[0],s[1],s[2],t);
        }
        else {
            egsWarning(ebox_message1,ebox_message4);
            if (t) {
                delete t;
            }
            return 0;
        }
        if (t) {
            delete t;
        }
        result->setName(input);
        result->setBoundaryTolerance(input);
        result->setMedia(input);
        result->setLabels(input);
        return result;
    }

}

The Makefile:

###############################################################################
#
#  EGSnrc egs++ makefile to build box 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:
#
###############################################################################


include $(EGS_CONFIG)
include $(SPEC_DIR)egspp.spec
include $(SPEC_DIR)egspp_$(my_machine).conf

DEFS = $(DEF1) -DBUILD_BOX_DLL

library = egs_box
lib_files = egs_box
my_deps = egs_transformations.h
extra_dep = $(addprefix $(DSOLIBS), $(my_deps))

include $(SPEC_DIR)egspp_libs.spec

$(make_depend)

/*
###############################################################################
#
#  EGSnrc egs++ box 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
#
###############################################################################
*/


#include "egs_box.h"
#include "egs_input.h"

void EGS_Box::printInfo() const {
    EGS_BaseGeometry::printInfo();
    egsInformation(" box size = %g %g %g\n",ax,ay,az);
    egsInformation("=======================================================\n");
}

string EGS_Box::type("EGS_Box");

static char EGS_BOX_LOCAL ebox_message1[] = "createGeometry(box): %s\n";
static char EGS_BOX_LOCAL ebox_message2[] = "null input?";
static char EGS_BOX_LOCAL ebox_message3[] = "wrong/missing 'box size' input?";
static char EGS_BOX_LOCAL ebox_message4[] =
    "expecting 1 or 3 float inputs for 'box size'";
static char EGS_BOX_LOCAL ebox_key1[] = "box size";

extern "C" {

    EGS_BOX_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
        if (!input) {
            egsWarning(ebox_message1,ebox_message2);
            return 0;
        }
        vector<EGS_Float> s;
        int err = input->getInput(ebox_key1,s);
        if (err) {
            egsWarning(ebox_message1,ebox_message3);
            return 0;
        }
        EGS_AffineTransform *t = EGS_AffineTransform::getTransformation(input);
        EGS_Box *result;
        if (s.size() == 1) {
            result = new EGS_Box(s[0],t);
        }
        else if (s.size() == 3) {
            result = new EGS_Box(s[0],s[1],s[2],t);
        }
        else {
            egsWarning(ebox_message1,ebox_message4);
            if (t) {
                delete t;
            }
            return 0;
        }
        if (t) {
            delete t;
        }
        result->setName(input);
        result->setBoundaryTolerance(input);
        result->setMedia(input);
        result->setLabels(input);
        return result;
    }

}