egs_dynamic_source.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ dynamic source 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:          Blake Walters, 2017
#
#  Contributors:    Reid Townson
#                   Alexandre Demelo
#
###############################################################################
*/
 
 
#ifndef EGS_DYNAMIC_SOURCE_
#define EGS_DYNAMIC_SOURCE_
 
#include "egs_vector.h"
#include "egs_base_source.h"
#include "egs_rndm.h"
#include "egs_shapes.h"
#include <string>
#include <sstream>
 
 
#ifdef WIN32
 
    #ifdef BUILD_DYNAMIC_SOURCE_DLL
        #define EGS_DYNAMIC_SOURCE_EXPORT __declspec(dllexport)
    #else
        #define EGS_DYNAMIC_SOURCE_EXPORT __declspec(dllimport)
    #endif
    #define EGS_DYNAMIC_SOURCE_LOCAL
 
#else
 
    #ifdef HAVE_VISIBILITY
        #define EGS_DYNAMIC_SOURCE_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_DYNAMIC_SOURCE_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_DYNAMIC_SOURCE_EXPORT
        #define EGS_DYNAMIC_SOURCE_LOCAL
    #endif
 
#endif
 
class EGS_DYNAMIC_SOURCE_EXPORT EGS_DynamicSource :
    public EGS_BaseSource {
 
public:
 
 
    struct EGS_ControlPoint {
 
        EGS_Vector iso; //isocentre position
        EGS_Float dsource; //source-isocentre distance
        EGS_Float theta; //angle of rotation about Y-axis
        EGS_Float phi;  //angle of rotation about Z-axis
        EGS_Float phicol; //angle of rotation in source plane
        EGS_Float time; //monitor unit index for control point
    };
 
    EGS_DynamicSource(EGS_BaseSource *Source, vector<EGS_ControlPoint> cpts,
                      const string &Name="", EGS_ObjectFactory *f=0) :
        EGS_BaseSource(Name,f), source(Source), valid(true) {
        //do some checks on cpts
        if (cpts.size()<2) {
            egsWarning("EGS_DynamicSource: not enough or missing control points.\n");
            valid = false;
        }
        else {
            if (cpts[0].time > 0.0) {
                egsWarning("EGS_DynamicSource: time index of control point 1 > 0.0.  This will generate many warning messages.\n");
            }
            int npts = cpts.size();
            for (int i=0; i<npts; i++) {
                if (i>0 && cpts[i].time < cpts[i-1].time) {
                    egsWarning("EGS_DynamicSource: time index of control point %i < time index of control point %i\n",i,i-1);
                    valid = false;
                }
                if (cpts[i].time<0.0) {
                    egsWarning("EGS_DynamicSource: time index of control point %i < 0.0\n",i);
                    valid = false;
                }
            }
            //normalize time values
            for (int i=0; i<npts-1; i++) {
                cpts[i].time /= cpts[npts-1].time;
            }
        }
        setUp();
    };
 
    EGS_DynamicSource(EGS_Input *, EGS_ObjectFactory *f=0);
 
    ~EGS_DynamicSource() {
        EGS_Object::deleteObject(source);
    };
 
    EGS_I64 getNextParticle(EGS_RandomGenerator *rndm,
                            int &q, int &latch, EGS_Float &E, EGS_Float &wt,
                            EGS_Vector &x, EGS_Vector &u) {
        int err = 1;
        EGS_ControlPoint ipt;
        EGS_I64 c;
        while (err) {
            c = source->getNextParticle(rndm,q,latch,E,wt,x,u);
            if (sync) {
                ptime = source->getTimeIndex();
                if (ptime<0) {
                    egsWarning("EGS_DynamicSource: You have selected synchronization of dynamic source with %s\n",source->getObjectName().c_str());
                    egsWarning("However, this source does not return time values for each particle.  Will turn off synchronization.\n");
                    sync = false;
                }
            }
            if (!sync) {
                ptime = rndm->getUniform();
            }
            setTimeIndex(ptime); //this is added for the storing of time index in a single location. Technically stored as a basesource attribute not dynamic source
            err = getCoord(ptime,ipt);
        }
 
        //translate source in Z
        x.z=x.z-ipt.dsource;
        //get the rotation matrices
        ipt.phicol *= M_PI/180;
        ipt.theta *= M_PI/180;
        ipt.phi *= M_PI/180;
        EGS_RotationMatrix Rcol(EGS_RotationMatrix::rotZ(ipt.phicol));
        EGS_RotationMatrix Rtheta(EGS_RotationMatrix::rotY(ipt.theta));
        EGS_RotationMatrix Rphi(EGS_RotationMatrix::rotZ(ipt.phi));
        //apply rotations in specific order and then translate relative
        //to the isocentre
        u=Rphi*Rtheta*Rcol*u;
        x=Rphi*Rtheta*Rcol*x + ipt.iso;
        return c;
    };
    EGS_Float getEmax() const {
        return source->getEmax();
    };
    EGS_Float getFluence() const {
        return source->getFluence();
    };
    bool storeState(ostream &data) const {
        return source->storeState(data);
    };
    bool setState(istream &data) {
        return source->setState(data);
    };
    bool addState(istream &data_in) {
        return source->addState(data_in);
    };
    void resetCounter() {
        source->resetCounter();
    };
 
    bool isValid() const {
        return (valid && source != 0);
    };
 
    void setSimulationChunk(EGS_I64 nstart, EGS_I64 nrun, int npar, int nchunk) {
        source->setSimulationChunk(nstart, nrun, npar, nchunk);
    };
 
    void containsDynamic(bool &hasdynamic);
 
protected:
 
    EGS_BaseSource *source; 
 
    vector<EGS_ControlPoint> cpts;  //control point
 
    int ncpts;  //no. of control points
 
    bool valid; //is this a valid source
 
    bool sync; //set to true if source motion synched with time read from
    //iaea phsp or beam simulation source
 
    int getCoord(const EGS_Float rand, EGS_ControlPoint &ipt);
 
    EGS_Float ptime; //time index corresponding to particle
    //could just be a random number.
 
    void setUp();
 
};
 
#endif