BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSArray4DCoords.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSArray4DCoords.hh"
20#include "BDSDebug.hh"
21#include "BDSException.hh"
22#include "BDSExtent.hh"
23#include "BDSFieldValue.hh"
24#include "BDSUtilities.hh"
25
26#include "CLHEP/Units/PhysicalConstants.h"
27
28#include <algorithm>
29#include <cmath>
30#include <limits>
31#include <ostream>
32#include <string>
33
34#include "globals.hh"
35
36BDSArray4DCoords::BDSArray4DCoords(G4int nXIn, G4int nYIn, G4int nZIn, G4int nTIn,
37 G4double xMinIn, G4double xMaxIn,
38 G4double yMinIn, G4double yMaxIn,
39 G4double zMinIn, G4double zMaxIn,
40 G4double tMinIn, G4double tMaxIn,
41 BDSDimensionType xDimensionIn,
42 BDSDimensionType yDimensionIn,
43 BDSDimensionType zDimensionIn,
44 BDSDimensionType tDimensionIn):
45 BDSArray4D(nXIn,nYIn,nZIn,nTIn),
46 xMin(xMinIn), xMax(xMaxIn),
47 yMin(yMinIn), yMax(yMaxIn),
48 zMin(zMinIn), zMax(zMaxIn),
49 tMin(tMinIn), tMax(tMaxIn),
50 smallestSpatialStep(std::numeric_limits<double>::max()),
51 xDimension(xDimensionIn),
52 yDimension(yDimensionIn),
53 zDimension(zDimensionIn),
54 tDimension(tDimensionIn),
55 dimensions{0,1,2,3},
56 timeVarying(true)
57{
58 // There are 1 fewer differences than the points.
59 if (nX > 1)
60 {
61 xStep = (xMax - xMin) / ((G4double)nX - 1);
62 CheckStep(xStep, "x");
63 smallestSpatialStep = std::min(smallestSpatialStep, xStep);
64 }
65 else
66 {xStep = 1;}
67 if (nY > 1)
68 {
69 yStep = (yMax - yMin) / ((G4double)nY - 1);
70 CheckStep(yStep, "y");
71 smallestSpatialStep = std::min(smallestSpatialStep, yStep);
72 }
73 else
74 {yStep = 1;}
75 if (nZ > 1)
76 {
77 zStep = (zMax - zMin) / ((G4double)nZ - 1);
78 CheckStep(zStep, "z");
79 smallestSpatialStep = std::min(smallestSpatialStep, zStep);
80 }
81 else
82 {zStep = 1;}
83 if (nT > 1)
84 {
85 tStep = (tMax - tMin) / ((G4double)nT - 1);
86 CheckStep(tStep, "t");
87 G4double lengthScale = CLHEP::c_light * tStep;
88 smallestSpatialStep = std::min(smallestSpatialStep, lengthScale);
89 }
90 else
91 {tStep = 1;}
93}
94
95void BDSArray4DCoords::CheckStep(G4double step, const G4String& name)
96{
97 if (!BDS::IsFinite(step))
98 {
99 throw BDSException(__METHOD_NAME__, "Invalid " + name + "min and " + name +
100 "max in array leading to 0 step size between points.");
101 }
102}
103
105 G4double y,
106 G4double z,
107 G4double t) const
108{
109 G4bool rx = x < xMin || x > xMax;
110 G4bool ry = y < yMin || y > yMax;
111 G4bool rz = z < zMin || z > zMax;
112 G4bool rt = t < tMin || t > tMax;
113 return rx || ry || rz || rt;
114}
115
117 G4double y,
118 G4double z,
119 G4double t) const
120{
121 if (OutsideCoords(x,y,z,t))
122 {
123 throw BDSException(__METHOD_NAME__, "(" +
124 std::to_string(x) + ", " +
125 std::to_string(y) + ", " +
126 std::to_string(z) + ", " +
127 std::to_string(t) + ") is outside array");
128 }
129}
130
132 G4double y,
133 G4double z,
134 G4double t,
135 BDSFieldValue (&localData)[2][2][2][2],
136 G4double& xFrac,
137 G4double& yFrac,
138 G4double& zFrac,
139 G4double& tFrac) const
140{
141 G4double xArrayCoords, yArrayCoords, zArrayCoords, tArrayCoords;
142 ArrayCoordsFromXYZT(x, xArrayCoords, y, yArrayCoords, z, zArrayCoords, t, tArrayCoords);
143 auto x1 = (G4int)std::floor(xArrayCoords);
144 auto y1 = (G4int)std::floor(yArrayCoords);
145 auto z1 = (G4int)std::floor(zArrayCoords);
146 auto t1 = (G4int)std::floor(tArrayCoords);
147
148 xFrac = xArrayCoords - x1;
149 yFrac = yArrayCoords - y1;
150 zFrac = zArrayCoords - z1;
151 tFrac = tArrayCoords - t1;
152
153 for (G4int i = 0; i < 2; i++)
154 {
155 for (G4int j = 0; j < 2; j++)
156 {
157 for (G4int k = 0; k < 2; k++)
158 {
159 for (G4int l = 0; l < 2; l++)
160 {localData[i][j][k][l] = GetConst(x1+i, y1+j, z1+k, t1+l);}
161 }
162 }
163 }
164}
165
167 G4double y,
168 G4double z,
169 G4double t,
170 BDSFieldValue (&localData)[4][4][4][4],
171 G4double& xFrac,
172 G4double& yFrac,
173 G4double& zFrac,
174 G4double& tFrac) const
175{
176 G4double xArrayCoords, yArrayCoords, zArrayCoords, tArrayCoords;
177 ArrayCoordsFromXYZT(x, xArrayCoords, y, yArrayCoords, z, zArrayCoords, t, tArrayCoords);
178 auto x1 = (G4int)std::floor(xArrayCoords);
179 auto y1 = (G4int)std::floor(yArrayCoords);
180 auto z1 = (G4int)std::floor(zArrayCoords);
181 auto t1 = (G4int)std::floor(tArrayCoords);
182
183 xFrac = xArrayCoords - x1;
184 yFrac = yArrayCoords - y1;
185 zFrac = zArrayCoords - z1;
186 tFrac = tArrayCoords - t1;
187
188 for (G4int i = 0; i < 4; i++)
189 {
190 for (G4int j = 0; j < 4; j++)
191 {
192 for (G4int k = 0; k < 4; k++)
193 {
194 for (G4int l = 0; l < 4; l++)
195 {localData[i][j][k][l] = GetConst(x1-1+i, y1-1+j, z1-1+k, t1-1+l);}
196 }
197 }
198 }
199}
200
202 G4double y,
203 G4double z,
204 G4double t) const
205{
206 G4int xind = NearestX(x);
207 G4int yind = NearestY(y);
208 G4int zind = NearestZ(z);
209 G4int tind = NearestT(t);
210 BDSFieldValue result = GetConst(xind, yind, zind, tind); // here we're constructing a copy on purpose
211 return result;
212}
213
214std::ostream& BDSArray4DCoords::Print(std::ostream& out) const
215{
216 out << "X: (" << xMin << ", " << xMax << ")" << G4endl;
217 out << "Y: (" << yMin << ", " << yMax << ")" << G4endl;
218 out << "Z: (" << zMin << ", " << zMax << ")" << G4endl;
219 out << "T: (" << tMin << ", " << tMax << ")" << G4endl;
220 return BDSArray4D::Print(out);
221}
222
223std::ostream& operator<< (std::ostream& out, BDSArray4DCoords const &a)
224{
225 return a.Print(out);
226}
227
229{
230 return BDSExtent(xMin, xMax, yMin, yMax, zMin, zMax);
231}
232
234{
235 dimensions[0] = FirstDimension().underlying();
236 dimensions[1] = SecondDimension().underlying();
237 dimensions[2] = ThirdDimension().underlying();
238 dimensions[3] = FourthDimension().underlying();
239}
240
242{
243 G4int result = -1;
244 G4int match = spatialDimension.underlying();
245 for (G4int i = 0; i < 4; i++)
246 {
247 if (dimensions[i] == match)
248 {
249 result = i;
250 break;
251 }
252 }
253 return result;
254}
Overlay of 4D array that provides uniform only spatial coordinate mapping.
G4double xMin
Dimension parameter - protected for derived class access.
virtual G4bool OutsideCoords(G4double x, G4double y, G4double z, G4double t) const
Whether the spatial coordinates lie outside the range of the array or not.
G4double yMin
Dimension parameter - protected for derived class access.
G4double yStep
Dimension parameter - protected for derived class access.
BDSDimensionType ThirdDimension() const
Accessor for each dimension label. e.g. array 'x' = spatial z.
virtual std::ostream & Print(std::ostream &out) const
virtual G4int NearestY(G4double y) const
Not much point in being both virtual and inline (in our use case) but has to be virtual.
virtual G4int NearestX(G4double x) const
Not much point in being both virtual and inline (in our use case) but has to be virtual.
BDSDimensionType SecondDimension() const
Accessor for each dimension label. e.g. array 'x' = spatial z.
virtual BDSExtent Extent() const
Return the SPATIAL (only) extent of this field without any offset. Ignores time.
G4double zMax
Dimension parameter - protected for derived class access.
virtual G4int NearestT(G4double t) const
Not much point in being both virtual and inline (in our use case) but has to be virtual.
G4double tMin
Dimension parameter - protected for derived class access.
BDSDimensionType FourthDimension() const
Accessor for each dimension label. e.g. array 'x' = spatial z.
G4double zMin
Dimension parameter - protected for derived class access.
virtual G4int NearestZ(G4double z) const
Not much point in being both virtual and inline (in our use case) but has to be virtual.
G4double xStep
Dimension parameter - protected for derived class access.
virtual void OutsideCoordsWarn(G4double x, G4double y, G4double z, G4double t) const
virtual void ExtractSection2x2x2x2(G4double x, G4double y, G4double z, G4double t, BDSFieldValue(&localData)[2][2][2][2], G4double &xFrac, G4double &yFrac, G4double &zFrac, G4double &tFrac) const
Extract 2x2x2x2 points lying around coordinate x.
G4int DimensionIndex(BDSDimensionType spatialDimension) const
G4double yMax
Dimension parameter - protected for derived class access.
G4double xMax
Dimension parameter - protected for derived class access.
virtual void ExtractSection4x4x4x4(G4double x, G4double y, G4double z, G4double t, BDSFieldValue(&localData)[4][4][4][4], G4double &xFrac, G4double &yFrac, G4double &zFrac, G4double &tFrac) const
Extract 4x4x4x4 points lying around coordinate x.
G4double tMax
Dimension parameter - protected for derived class access.
G4double zStep
Dimension parameter - protected for derived class access.
void ArrayCoordsFromXYZT(G4double &x, G4double &xArr, G4double &y, G4double &yArr, G4double &z, G4double &zArr, G4double &t, G4double &tArr) const
Utility version to forward to individual function.
BDSArray4DCoords()=delete
virtual BDSFieldValue ExtractNearest(G4double x, G4double y=0, G4double z=0, G4double t=0) const
Extract nearest field value from array.
G4double tStep
Dimension parameter - protected for derived class access.
BDSDimensionType FirstDimension() const
Accessor for each dimension label. e.g. array 'x' = spatial z.
void BuildDimensionIndex()
Build up an array of ints based on the order of dimensions stored in the array.
4D array and base class for 3,2 & 1D arrays.
Definition: BDSArray4D.hh:51
const G4int nZ
Dimension.
Definition: BDSArray4D.hh:122
const G4int nY
Dimension.
Definition: BDSArray4D.hh:121
virtual const BDSFieldValue & GetConst(G4int x, G4int y=0, G4int z=0, G4int t=0) const
Definition: BDSArray4D.cc:46
const G4int nX
Dimension.
Definition: BDSArray4D.hh:120
virtual std::ostream & Print(std::ostream &out) const
Definition: BDSArray4D.cc:91
const G4int nT
Dimension.
Definition: BDSArray4D.hh:123
General exception with possible name of object and message.
Definition: BDSException.hh:35
Holder for +- extents in 3 dimensions.
Definition: BDSExtent.hh:39
type underlying() const
return underlying value (can be used in switch statement)
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())