BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSCavityFactoryElliptical.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 "BDSCavity.hh"
20#include "BDSCavityFactoryElliptical.hh"
21#include "BDSCavityInfo.hh"
22
23#include "globals.hh"
24#include "G4LogicalVolume.hh"
25#include "G4Polycone.hh"
26#include "G4SubtractionSolid.hh"
27#include "G4Tubs.hh"
28#include "G4VisAttributes.hh"
29#include "G4VSolid.hh"
30
31#include "CLHEP/Units/SystemOfUnits.h"
32
33#include <cmath>
34#include <vector>
35
36BDSCavityFactoryElliptical::BDSCavityFactoryElliptical()
37{;}
38
39BDSCavityFactoryElliptical::~BDSCavityFactoryElliptical()
40{;}
41
43 G4double totalChordLength,
44 const BDSCavityInfo* info)
45{
46 G4double chordLength = totalChordLength;
47 // Elliptical Cavity Parameters
48 // irisRSemiAxis --> Semi-axis of the iris ellipse perpendicular to the length of the cavity.
49 // irisZSemiAxis --> Semi-axis of the iris ellipse along the length of the cavity.
50 // equatorRSemiAxis --> Semi-axis of the equator ellipse perpendicular to the length of the cavity.
51 // equatorZSemiAxis --> Semi-axis of the equator ellipse along the length of the cavity.
52 // tangentAngle --> Angle between the line vertical to the z-axis and the line tangent to both ellipses.
53 // irisRadius --> Radius of the iris.
54 // noPoints --> Number of points used to define the shape. This is divided up between
55 // the 3 constituent shapes, 1/4 in each iris ellipse and 1/2 in central dome.
56 // length --> Length of the cell.
57 // equatorRadius --> Radius at the equator.
58 G4double equatorRadius = info->equatorRadius;
59 G4double irisRSemiAxis = info->irisVerticalAxis;
60 G4double irisZSemiAxis = info->irisHorizontalAxis; //iris ellipse horizontal semiaxis
61 G4double equatorRSemiAxis = info->equatorVerticalAxis ; //equator ellipse vertical semiaxis
62 G4double equatorZSemiAxis = info->equatorHorizontalAxis; //equator ellipse horizontal semiaxis
63 G4double tangentAngle = info->tangentLineAngle;
64 G4double irisRadius = info->irisRadius;
65 G4double thickness = info->thickness;
66 unsigned int noPoints = info->numberOfPoints;
67
68 // calculate cartesian coordinates (z, r) from parameters.
69 // 2D cylindrical coordinates, z along the beamline:
70 G4double zi = chordLength * 0.5; // z coord of iris ellipse centre
71 G4double ri = irisRadius + irisRSemiAxis; // r coord of iris ellipse centre
72 G4double ze = 0.0; // z coord of equator ellipse centre
73 G4double re = equatorRadius - equatorRSemiAxis; // r coord of equator ellipse centre.
74 // gradient of line connecting the ellipses.
75 // Add a pi/2 because angle is defined from the vertical, clockwise.
76 G4double m = std::tan(tangentAngle + CLHEP::halfpi);
77
78 // Gradient from tangentAngle. Find the derivative of ellipses. Equate
79 // and solve for the parameter.
80 // atan finds solution in the first quadrant.
81 G4double equatorParameterTangentPoint = std::atan(-equatorRSemiAxis/(m*equatorZSemiAxis));
82 // Add pi to get desired solution (third quadrant)
83 G4double irisParameterTangentPoint = std::atan(-irisRSemiAxis/(m*irisZSemiAxis)) + CLHEP::pi;
84
85 // Rounding down to a multiple of 4. This is so that number of points are
86 // share equally and consistently between the constituent ellipses.
87 noPoints = noPoints - (noPoints % 4);
88
89 // Vector Definitions:
90 // Parametric equation parameters:
91 // equatorParameter --> the values of a parameter for generating the central ellipses
92 // irisParameter --> the values of a parameter for generating the iris ellipse
93
94 //The Cartesian coordinates to be calculated from above parameters:
95 //rInnerCoord --> the values of radius for the inner geometry
96 //rOuterCoord --> the values of radius for the outer geometry
97 //zInnerCoord --> the values of z for the inner geometry
98 //zOuterCoord --> the values of z for the outer goemetry
99
100 //Vector declaration:
101 std::vector<G4double> equatorParameter;
102 std::vector<G4double> irisParameter;
103 std::vector<G4double> rInnerCoord;
104 std::vector<G4double> rOuterCoord;
105 std::vector<G4double> zInnerCoord;
106 std::vector<G4double> zOuterCoord;
107
108 // Generate values for the parameters used to define the ellipse shapes
109
110 // Vector of parameter values for generating the iris ellipse. This can then be
111 // reflected to find the iris ellipse on the other side.
112 for (unsigned int i = 0; i < noPoints/4; i++) //iris allocated one quarter of points.
113 {irisParameter.push_back(1.5*CLHEP::pi + (1.5*CLHEP::pi - irisParameterTangentPoint )*i/((noPoints/4)-1));}
114
115 //Vector of parameter values for generating the equator ellipse
116 for (unsigned int i = 0; i < noPoints/2; i++) //equator ellipse allocated half of the points
117 {
118 G4double term1 = CLHEP::pi - equatorParameterTangentPoint;
119 G4double term2 = CLHEP::pi - 2*equatorParameterTangentPoint;
120 G4double term3 = (noPoints*0.5) - 1;
121 G4double value = term1 - term2 * i / term3;
122 equatorParameter.push_back(value);
123 }
124
125 // Calculate the coordinates of the inner cavity shape, starting from -z through
126 // to +z, in 3 sections: -z iris ellipse -> equator ellipse -> +z iris ellipse.
127 // Geant requires the coordinates to be done in this (negative to positive)
128
129 // -z iris ellipse:
130 // noPoints/4 because of the total noPoints, the iris ellipse in -z is allocated a quarter of points.
131 for (unsigned int i = 0; i < noPoints/4; i++)
132 {
133 // An extra point at the very start for unambigious boolean subtraction later on.
134 // Will not actually be part of the final geometry.
135 if (i == 0)
136 {
137 //-10mm so boundaries don't lay upon each other when using boolean operation
138 zInnerCoord.push_back(- zi + irisZSemiAxis * std::cos(irisParameter[i]) - 10*CLHEP::mm);
139 rInnerCoord.push_back(ri + irisRSemiAxis * std::sin(irisParameter[i]));
140 }
141 zInnerCoord.push_back(- zi + irisZSemiAxis * std::cos(irisParameter[i]));
142 rInnerCoord.push_back(ri + irisRSemiAxis * std::sin(irisParameter[i]));
143 }
144
145 //equator ellipse (middle):
146 for (unsigned int i = 0; i < noPoints/2; i++) //central ellipse is allocated half the total points.
147 {
148 zInnerCoord.push_back(ze + equatorZSemiAxis * std::cos(equatorParameter[i]));
149 rInnerCoord.push_back(re + equatorRSemiAxis * std::sin(equatorParameter[i]));
150 };
151
152 //+z iris ellipse:
153 for (unsigned int i = 0; i < noPoints/4; i++) //the +z iris ellipse is allocated noPoints/4.
154 {
155 zInnerCoord.push_back(zi + irisZSemiAxis * std::cos(CLHEP::pi - irisParameter[noPoints/4 - 1 - i]));
156 rInnerCoord.push_back(ri + irisRSemiAxis * std::sin(CLHEP::pi - irisParameter[noPoints/4 - 1 - i]));
157 //If last point, add an extra point purely for unambigious boolean subtraction later.
158 if (i == (noPoints/4 - 1))
159 {
160 //+10mm so boundaries don't lay upon each other
161 zInnerCoord.push_back(zi + irisZSemiAxis * std::cos(CLHEP::pi - irisParameter[noPoints/4 - 1 - i]) + 10*CLHEP::mm);
162 rInnerCoord.push_back(ri + irisRSemiAxis * std::sin(CLHEP::pi - irisParameter[noPoints/4 - 1 - i]));
163 }
164 }
165
166 // Generating points for the entire outer cavity shape by using the inner cavity
167 // vector of coordinates--- This works by calculating the gradient at each point on the curve.
168 // It then finds a unit vector pointing perpendicularly outwards from this point and multiplies
169 // by the thickness. This creates a second vector of points for the outer geometry.
170 //starting from 1 and ending at n-1 so that the coordinates used for the boolean subtraction are excluded.
171 for (unsigned int i = 1; i < zInnerCoord.size() - 1; i++)
172 {
173 //calculate the gradient at point i
174 G4double gradientAt_i = (rInnerCoord[i-1]-rInnerCoord[i+1])/(zInnerCoord[i-1]-zInnerCoord[i+1]);
175
176 // Since the calculation relies on using the two points either side, and the
177 // ends have only one coordinate either side, take the ends as having gradient = 0.
178 // So that the vector points vertically upwards. Otherwise just do the calculation described above.
179 if (i == 1 || i == (zInnerCoord.size() - 2))
180 {
181 zOuterCoord.push_back(zInnerCoord[i]);
182 rOuterCoord.push_back(rInnerCoord[i] + thickness);
183 }
184 else
185 {
186 zOuterCoord.push_back(- gradientAt_i*thickness * std::pow( (gradientAt_i * gradientAt_i + 1), -0.5) + zInnerCoord[i]);
187 rOuterCoord.push_back(+ thickness * std::pow( (gradientAt_i*gradientAt_i + 1), -0.5) + rInnerCoord[i]);
188 }
189 }
190
191#ifdef BDSDEBUG
192 G4cout << "Now printing the values of (zInnerCoord,rInnerCoord):" << G4endl;
193 G4cout << "Length of zInnerCoord = " << zInnerCoord.size() << G4endl;
194 G4cout << "Length of rInnerCoord = " << rInnerCoord.size() << G4endl;
195 for (G4int i = 0; i < (G4int)zInnerCoord.size(); i++)
196 {G4cout << "(" << zInnerCoord[i] << "," << rInnerCoord[i] << ")" << G4endl;};
197
198 G4cout << "Now printing the values of (zOuterCoord,rOuterCoord):" << G4endl;
199 G4cout << "Length of zOuterCoord = " << zOuterCoord.size() << G4endl;
200 G4cout << "Length of rOuterCoord = " << rOuterCoord.size() << G4endl;
201 for (G4int i = 0; i < (G4int)zOuterCoord.size(); i++)
202 {G4cout << "(" << zOuterCoord[i] << "," << rOuterCoord[i] << ")" << G4endl;};
203#endif
204
205 // Array of inner r coordinates. zeroes ensures the polycone will be solid.
206 std::vector<G4double> solidArrayOuter(noPoints, 0.0);
207 // Array of inner r coordinates. zeroes ensures the polycone will be solid. 1 extra
208 // point either side for unambiguous boolean subtraction.
209 std::vector<G4double> solidArrayInner(noPoints+2, 0.0);
210
211 // Define the inner solid which is to be subtracted from the outer and also used to
212 // define the vacuum.
213 G4VSolid* innerSolid = new G4Polycone(name + "_inner_solid", //name
214 0.0, //start angle
215 CLHEP::twopi, //sweep angle
216 zInnerCoord.size(),
217 zInnerCoord.data(),
218 solidArrayInner.data(),
219 rInnerCoord.data()); //r coordinates
220
221 // Define the outer solid from which the inner is subtracted.
222 G4VSolid* outerSolid = new G4Polycone(name + "_outer_solid", // name
223 0.0, // start angle
224 CLHEP::twopi, // sweep angle
225 zOuterCoord.size(), // number of corners
226 zOuterCoord.data(), // zOuterCoord array of points
227 solidArrayOuter.data(),// array of zeroes for solid shape
228 rOuterCoord.data()); // rOuterCoord array of points.
229 allSolids.insert(innerSolid);
230 allSolids.insert(outerSolid);
231
232 // Do the subtraction
233 cavitySolid = new G4SubtractionSolid(name + "_cavity_solid", // name
234 outerSolid, // solid1
235 innerSolid); // minus solid2
236
237 // Delete first and last elements of the InnerCoord Vectors, as these entries.
238 // The reason we need to do this is the same vector is also used for making the
239 // vacuum filling the cavity, but without the extra points for subtraction.
240 rInnerCoord.erase (rInnerCoord.begin());
241 rInnerCoord.pop_back();
242 zInnerCoord.erase (zInnerCoord.begin());
243 zInnerCoord.pop_back();
244
245 //Take lengthsafety from all r points and the end z points, to avoid overlap with the cavity.
246 for (unsigned int i = 0; i < noPoints; i++)
247 {rInnerCoord[i] = rInnerCoord[i] - lengthSafetyLarge;}
248
249 zInnerCoord[0] = zInnerCoord[0] + 2*lengthSafety;
250 G4int lastInd = (G4int)zInnerCoord.size() - 1;
251 zInnerCoord[lastInd] = zInnerCoord[lastInd] - 2*lengthSafety; // last element
252
253 vacuumSolid = new G4Polycone(name + "_inner_solid", //name
254 0.0, //start angle
255 CLHEP::twopi, //sweep angle
256 zInnerCoord.size(),
257 zInnerCoord.data(),
258 solidArrayInner.data(),
259 rInnerCoord.data()); //r coordinates
260 allSolids.insert(vacuumSolid);
261
262 G4double containerRadius = equatorRadius + thickness + lengthSafetyLarge;
263 containerSolid = new G4Tubs(name + "_container_solid", // name
264 0.0, // innerRadius
265 containerRadius, // outerRadius
266 chordLength*0.5, // half length
267 0.0, // starting angle
268 CLHEP::twopi); // sweep angle
269 return containerRadius;
270}
271
273{
275}
virtual void SetVisAttributes(G4String colourName="rfcavity")
Set vis attributes for cavityLV, vacuumLV and containerLV.
G4VSolid * vacuumSolid
Cache of particular solid or lv for common functionality in this class.
G4VSolid * containerSolid
Cache of particular solid or lv for common functionality in this class.
G4VSolid * cavitySolid
Cache of particular solid or lv for common functionality in this class.
virtual G4double CreateSolids(G4String name, G4double totalChordLength, const BDSCavityInfo *info)
Create vacuumSolid and cavitySolid. Must return the container radius.
virtual void SetVisAttributes(G4String colourName)
Overload base class method, but use base class method with different colour.
Holder for all Geometrical information required to create an RF cavity.
G4double equatorRadius
Equator radius - half width of widest part.
G4double equatorVerticalAxis
Equator ellipse vertical semi-axis.
G4double tangentLineAngle
Tangent angle.
G4double irisRadius
Iris radius - half width of narrowest part.
G4double irisVerticalAxis
Iris ellipse vertical semi-axis.
G4int numberOfPoints
G4double irisHorizontalAxis
Iris ellipse horizontal semi-axis.
G4double equatorHorizontalAxis
Equator ellipse horizontal semi-axis.
G4double thickness
Thickness of wall material.
G4double lengthSafety
Cache of global constants variable.
G4double lengthSafetyLarge
Cache of global constants variable.