Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
HadrontherapyElectricTabulatedField3D Class Reference

#include <HadrontherapyElectricTabulatedField3D.hh>

Inheritance diagram for HadrontherapyElectricTabulatedField3D:
Inheritance graph
[legend]
Collaboration diagram for HadrontherapyElectricTabulatedField3D:
Collaboration graph
[legend]

Public Member Functions

 HadrontherapyElectricTabulatedField3D (const char *filename, G4double exOffset, G4double eyOffset, G4double ezOffset)
 
void GetFieldValue (const G4double Epoint[4], G4double *Efield) const
 
- Public Member Functions inherited from G4ElectricField
 G4ElectricField ()
 
virtual ~G4ElectricField ()
 
 G4ElectricField (const G4ElectricField &r)
 
G4ElectricField & operator= (const G4ElectricField &p)
 
G4bool DoesFieldChangeEnergy () const
 
- Public Member Functions inherited from G4ElectroMagneticField
 G4ElectroMagneticField ()
 
virtual ~G4ElectroMagneticField ()
 
 G4ElectroMagneticField (const G4ElectroMagneticField &r)
 
G4ElectroMagneticField & operator= (const G4ElectroMagneticField &p)
 
- Public Member Functions inherited from G4Field
 G4Field (G4bool gravityOn=false)
 
 G4Field (const G4Field &)
 
virtual ~G4Field ()
 
G4Field & operator= (const G4Field &p)
 
G4bool IsGravityActive () const
 
void SetGravityActive (G4bool OnOffFlag)
 
virtual G4Field * Clone () const
 

Detailed Description

Definition at line 40 of file HadrontherapyElectricTabulatedField3D.hh.

Constructor & Destructor Documentation

HadrontherapyElectricTabulatedField3D::HadrontherapyElectricTabulatedField3D ( const char *  filename,
G4double  exOffset,
G4double  eyOffset,
G4double  ezOffset 
)

Definition at line 35 of file HadrontherapyElectricTabulatedField3D.cc.

36  :feXoffset(exOffset),feYoffset(eyOffset),feZoffset(ezOffset),einvertX(false),einvertY(false),einvertZ(false)
37 {
38  //The format file is: X Y Z Ex Ey Ez
39 
40  G4double ElenUnit= cm;
41  G4double EfieldUnit= volt/m;
42  G4cout << "\n-----------------------------------------------------------"
43  << "\n Electric field"
44  << "\n-----------------------------------------------------------";
45 
46  G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl;
47  G4AutoLock lock(&MyHadrontherapyLockEField);
48 
49  ifstream file( filename ); // Open the file for reading.
50 
51  // Ignore first blank line
52  char ebuffer[256];
53  file.getline(ebuffer,256);
54 
55  // Read table dimensions
56  file >> Enx >> Eny >> Enz; // Note dodgy order
57 
58  G4cout << " [ Number of values x,y,z: "
59  << Enx << " " << Eny << " " << Enz << " ] "
60  << endl;
61 
62  // Set up storage space for table
63  xEField.resize( Enx );
64  yEField.resize( Enx );
65  zEField.resize( Enx );
66  G4int ix, iy, iz;
67  for (ix=0; ix<Enx; ix++) {
68  xEField[ix].resize(Eny);
69  yEField[ix].resize(Eny);
70  zEField[ix].resize(Eny);
71  for (iy=0; iy<Eny; iy++) {
72  xEField[ix][iy].resize(Enz);
73  yEField[ix][iy].resize(Enz);
74  zEField[ix][iy].resize(Enz);
75  }
76  }
77 
78  // Read in the data
79  G4double Exval=0.;
80  G4double Eyval=0.;
81  G4double Ezval=0.;
82  G4double Ex=0.;
83  G4double Ey=0.;
84  G4double Ez=0.;
85  for (iz=0; iz<Enz; iz++) {
86  for (iy=0; iy<Eny; iy++) {
87  for (ix=0; ix<Enx; ix++) {
88  file >> Exval >> Eyval >> Ezval >> Ex >> Ey >> Ez;
89 
90  if ( ix==0 && iy==0 && iz==0 ) {
91  Eminx = Exval * ElenUnit;
92  Eminy = Eyval * ElenUnit;
93  Eminz = Ezval * ElenUnit;
94  }
95  xEField[ix][iy][iz] = Ex * EfieldUnit;
96  yEField[ix][iy][iz] = Ey * EfieldUnit;
97  zEField[ix][iy][iz] = Ez * EfieldUnit;
98  }
99  }
100  }
101  file.close();
102  lock.unlock();
103 
104  Emaxx = Exval * ElenUnit;
105  Emaxy = Eyval * ElenUnit;
106  Emaxz = Ezval * ElenUnit;
107 
108  G4cout << "\n ---> ... done reading " << endl;
109 
110  // G4cout << " Read values of field from file " << filename << endl;
111  G4cout << " ---> assumed the order: x, y, z, Ex, Ey, Ez "
112  << "\n ---> Min values x,y,z: "
113  << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
114  << "\n ---> Max values x,y,z: "
115  << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm "
116  << "\n ---> The field will be offset in x by " << exOffset/cm << " cm "
117  << "\n ---> The field will be offset in y by " << eyOffset/cm << " cm "
118  << "\n ---> The field will be offset in z by " << ezOffset/cm << " cm " << endl;
119 
120  // Should really check that the limits are not the wrong way around.
121  if (Emaxx < Eminx) {swap(Emaxx,Eminx); einvertX = true;}
122  if (Emaxy < Eminy) {swap(Emaxy,Eminy); einvertY = true;}
123  if (Emaxz < Eminz) {swap(Emaxz,Eminz); einvertZ = true;}
124  G4cout << "\nAfter reordering if neccesary"
125  << "\n ---> Min values x,y,z: "
126  << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
127  << " \n ---> Max values x,y,z: "
128  << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm ";
129 
130  dx1 = Emaxx - Eminx;
131  dy1 = Emaxy - Eminy;
132  dz1 = Emaxz - Eminz;
133  G4cout << "\n ---> Dif values x,y,z (range): "
134  << dx1/cm << " " << dy1/cm << " " << dz1/cm << " cm "
135  << "\n-----------------------------------------------------------" << endl;
136 }
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
m
static constexpr double m
Definition: G4SIunits.hh:129
cm
static constexpr double cm
Definition: G4SIunits.hh:119
volt
static constexpr double volt
Definition: G4SIunits.hh:244
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Member Function Documentation

void HadrontherapyElectricTabulatedField3D::GetFieldValue ( const G4double  Epoint[4],
G4double *  Efield 
) const
virtual

Implements G4ElectricField.

Definition at line 138 of file HadrontherapyElectricTabulatedField3D.cc.

140 {
141  G4double x1 = Epoint[0] + feXoffset;
142  G4double y1 = Epoint[1] + feYoffset;
143  G4double z1 = Epoint[2] + feZoffset;
144 
145  // Position of given point within region, normalized to the range
146  // [0,1]
147  G4double Exfraction = (x1 - Eminx) / dx1;
148  G4double Eyfraction = (y1 - Eminy) / dy1;
149  G4double Ezfraction = (z1 - Eminz) / dz1;
150 
151  if (einvertX) { Exfraction = 1 - Exfraction;}
152  if (einvertY) { Eyfraction = 1 - Eyfraction;}
153  if (einvertZ) { Ezfraction = 1 - Ezfraction;}
154 
155  // Need addresses of these to pass to modf below.
156  // modf uses its second argument as an OUTPUT argument.
157  G4double exdindex, eydindex, ezdindex;
158 
159  // Position of the point within the cuboid defined by the
160  // nearest surrounding tabulated points
161  G4double exlocal = ( std::modf(Exfraction*(Enx-1), &exdindex));
162  G4double eylocal = ( std::modf(Eyfraction*(Eny-1), &eydindex));
163  G4double ezlocal = ( std::modf(Ezfraction*(Enz-1), &ezdindex));
164 
165  // The indices of the nearest tabulated point whose coordinates
166  // are all less than those of the given point
167  G4int exindex = static_cast<G4int>(std::floor(exdindex));
168  G4int eyindex = static_cast<G4int>(std::floor(eydindex));
169  G4int ezindex = static_cast<G4int>(std::floor(ezdindex));
170 
171  if ((exindex < 0) || (exindex >= Enx - 1) ||
172  (eyindex < 0) || (eyindex >= Eny - 1) ||
173  (ezindex < 0) || (ezindex >= Enz - 1))
174  {
175  Efield[0] = 0.0;
176  Efield[1] = 0.0;
177  Efield[2] = 0.0;
178  Efield[3] = 0.0;
179  Efield[4] = 0.0;
180  Efield[5] = 0.0;
181  }
182  else
183  {
184 
185 /*
186 #ifdef DEBUG_G4intERPOLATING_FIELD
187  G4cout << "Local x,y,z: " << exlocal << " " << eylocal << " " << ezlocal << endl;
188  G4cout << "Index x,y,z: " << exindex << " " << eyindex << " " << ezindex << endl;
189  G4double valx0z0, mulx0z0, valx1z0, mulx1z0;
190  G4double valx0z1, mulx0z1, valx1z1, mulx1z1;
191  valx0z0= table[exindex ][0][ezindex]; mulx0z0= (1-exlocal) * (1-ezlocal);
192  valx1z0= table[exindex+1][0][ezindex]; mulx1z0= exlocal * (1-ezlocal);
193  valx0z1= table[exindex ][0][ezindex+1]; mulx0z1= (1-exlocal) * ezlocal;
194  valx1z1= table[exindex+1][0][ezindex+1]; mulx1z1= exlocal * ezlocal;
195 #endif
196 */
197  // Full 3-dimensional version
198 
199  Efield[0] = 0.0;
200  Efield[1] = 0.0;
201  Efield[2] = 0.0;
202 
203  Efield[3] =
204  xEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
205  xEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
206  xEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
207  xEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
208  xEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
209  xEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
210  xEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
211  xEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
212  Efield[4] =
213  yEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
214  yEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
215  yEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
216  yEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
217  yEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
218  yEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
219  yEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
220  yEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
221  Efield[5] =
222  zEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
223  zEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
224  zEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
225  zEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
226  zEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
227  zEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
228  zEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
229  zEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
230  }
231 //G4cout << "Getting electric field " << Efield[3]/(volt/m) << " " << Efield[4]/(volt/m) << " " << Efield[5]/(volt/m) << endl;
232 //G4cout << "For coordinates: " << Epoint[0] << " " << Epoint[1] << " " << Epoint[2] << G4endl;
233 
234 /*std::ofstream WriteDataIn("ElectricFieldFC.out", std::ios::app);
235  WriteDataIn << Epoint[0] << '\t' << " "
236  << Epoint[1] << '\t' << " "
237  << Epoint[2] << '\t' << " "
238  << Efield[3]/(volt/m) << '\t' << " "
239  << Efield[4]/(volt/m) << '\t' << " "
240  << Efield[5]/(volt/m) << '\t' << " "
241  << G4endl; */
242 }
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
The documentation for this class was generated from the following files: