EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_dynamic_source.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ dynamic source
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Blake Walters, 2017
25 #
26 # Contributors:
27 #
28 ###############################################################################
29 */
30 
31 
37 #include "egs_dynamic_source.h"
38 #include "egs_input.h"
39 
41  EGS_ObjectFactory *f) : EGS_BaseSource(input,f), source(0), valid(true) {
42  EGS_Input *isource = input->takeInputItem("source",false);
43  if (isource) {
45  delete isource;
46  }
47  if (!source) {
48  string sname;
49  int err = input->getInput("source name",sname);
50  if (err)
51  egsWarning("EGS_DynamicSource: missing/wrong inline source "
52  "definition and missing wrong 'source name' input\n");
53  else {
55  if (!source) egsWarning("EGS_DynamicSource: a source named %s"
56  " does not exist\n",sname.c_str());
57  }
58  }
59  //now read inputs relevant to dynamic source
60  //see if user wants to synchronize source with mu read from
61  //iaea phsp or beam simulation source
62  vector<string> sync_options;
63  sync_options.push_back("no");
64  sync_options.push_back("yes");
65  sync = input->getInput("synchronize motion",sync_options,0);
66  if (sync && source->getObjectType()!="IAEA_PhspSource" &&
67  source->getObjectType()!="EGS_BeamSource") {
68  egsWarning("EGS_DynamicSource: source motion can only be synchronized with a source of type iaea_phsp_source or egs_beam_source.\n Will not synchronize.\n");
69  sync = false;
70  }
71 
72  //get control points
73  EGS_Input *dyninp = input->takeInputItem("motion");
74  if (dyninp) {
75  //get control points
76  ncpts=0;
77  vector<EGS_Float> point;
78  EGS_ControlPoint cpt;
79  stringstream itos;
80  int err;
81  int icpts=1;
82  itos << icpts;
83  string sstring = "control point " + itos.str();
84  while (!(err = dyninp->getInput(sstring,point))) {
85  if (point.size()!=8) {
86  egsWarning("EGS_DynamicSource: control point %i does not specify 8 values.\n",icpts);
87  valid = false;
88  }
89  else {
90  if (ncpts>0 && point[7] < cpts[ncpts-1].mu) {
91  egsWarning("EGS_DynamicSource: mu index of control point %i < mu index of control point %i\n",icpts,ncpts);
92  valid = false;
93  }
94  else if (point[7] < 0.) {
95  egsWarning("EGS_DynamicSource: mu index of control point %i < 0.0\n",icpts);
96  valid = false;
97  }
98  else {
99  ncpts++;
100  if (ncpts ==1 && point[7] > 0.0) {
101  egsWarning("EGS_DynamicSource: mu index of control point 1 > 0.0. This will generate many warning messages.\n");
102  }
103  //set cpt values
104  cpt.iso = EGS_Vector(point[0],point[1],point[2]);
105  cpt.dsource = point[3];
106  cpt.theta = point[4];
107  cpt.phi = point[5];
108  cpt.phicol = point[6];
109  cpt.mu = point[7];
110  //add it to the vector of control points
111  cpts.push_back(cpt);
112  icpts++;
113  itos.str("");
114  itos << icpts;
115  sstring = "control point " + itos.str();
116  }
117  }
118  }
119  if (ncpts<=1) {
120  egsWarning("EGS_DynamicSource: not enough or missing control points.\n");
121  valid = false;
122  }
123  if (cpts[ncpts-1].mu == 0.0) {
124  egsWarning("EGS_DynamicSource: mu index of last control point = 0. Something's wrong.\n");
125  valid = false;
126  }
127  else {
128  //normalize mu index to max. value
129  for (int i=0; i<ncpts-1; i++) {
130  cpts[i].mu /= cpts[ncpts-1].mu;
131  }
132  }
133  }
134  else {
135  egsWarning("EGS_DynamicSource: no control points input.\n");
136  valid = false;
137  }
138  setUp();
139 }
140 
141 void EGS_DynamicSource::setUp() {
142  //most setup done in constructor
143  otype="EGS_DynamicSource";
144  if (!isValid()) {
145  description = "Invalid dynamic source";
146  }
147  else {
148  description = "Dynamic source based on\n";
150  if (sync) {
151  description += "\n Source will be synched with mu values read in (if available).";
152  }
153  }
154 }
155 
156 //actually select the rotation coordinates for the incident particle
157 int EGS_DynamicSource::getCoord(EGS_Float rand, EGS_ControlPoint &ipt) {
158  int iindex=0;
159  int i;
160  for (i=0; i<ncpts; i++) {
161  if (rand < cpts[i].mu) {
162  iindex =i;
163  break;
164  }
165  }
166  if (i==ncpts) {
167  egsWarning("EGS_DynamicSource: could not locate control point.\n");
168  return 1;
169  }
170  EGS_Float factor = (rand-cpts[iindex-1].mu)/(cpts[iindex].mu-cpts[iindex-1].mu);
171  ipt.iso.x=cpts[iindex-1].iso.x+ (cpts[iindex].iso.x-cpts[iindex-1].iso.x)*factor;
172  ipt.iso.y=cpts[iindex-1].iso.y+ (cpts[iindex].iso.y-cpts[iindex-1].iso.y)*factor;
173  ipt.iso.z=cpts[iindex-1].iso.z+ (cpts[iindex].iso.z-cpts[iindex-1].iso.z)*factor;
174  ipt.dsource=cpts[iindex-1].dsource+ (cpts[iindex].dsource-cpts[iindex-1].dsource)*factor;
175  ipt.theta=cpts[iindex-1].theta+ (cpts[iindex].theta-cpts[iindex-1].theta)*factor;
176  ipt.phi=cpts[iindex-1].phi+ (cpts[iindex].phi-cpts[iindex-1].phi)*factor;
177  ipt.phicol=cpts[iindex-1].phicol+ (cpts[iindex].phicol-cpts[iindex-1].phicol)*factor;
178  return 0;
179 };
180 
181 extern "C" {
182 
183  EGS_DYNAMIC_SOURCE_EXPORT EGS_BaseSource *createSource(EGS_Input *input,
184  EGS_ObjectFactory *f) {
185  return
186  createSourceTemplate<EGS_DynamicSource>(input,f,"dynamic source");
187  }
188 
189 }
static EGS_BaseSource * createSource(EGS_Input *)
Create sources from the information pointed to by input.
A source with simulated time-varying rotations/translations.
EGS_Input class header file.
A class representing 3D vectors.
Definition: egs_vector.h:56
const char * getSourceDescription() const
Get a short description of this source.
EGS_DynamicSource(EGS_BaseSource *Source, vector< EGS_ControlPoint > cpts, const string &Name="", EGS_ObjectFactory *f=0)
Construct a dynamic source using Source as the source and cpts as the control points. Not sure if this will ever be used but here just in case.
An object factory.
static EGS_BaseSource * getSource(const string &Name)
Get a pointer to the source named Name.
const string & getObjectType() const
Get the object type.
EGS_BaseSource * source
The source being rotated.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
string otype
The object type.
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
Base source class. All particle sources must be derived from this class.
string description
A short source description.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.