Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/purging_magnet/src/PurgMagTabulatedField3D.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/purging_magnet/src/PurgMagTabulatedField3D.cc (Version 11.3.0) and /examples/advanced/purging_magnet/src/PurgMagTabulatedField3D.cc (Version 4.0.p2)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // Code developed by:                             
 27 //  S.Larsson and J. Generowicz.                  
 28 //                                                
 29 //    *************************************       
 30 //    *                                   *       
 31 //    *    PurgMagTabulatedField3D.cc     *       
 32 //    *                                   *       
 33 //    *************************************       
 34 //                                                
 35 //                                                
 36                                                   
 37 #include "PurgMagTabulatedField3D.hh"             
 38 #include "G4SystemOfUnits.hh"                     
 39 #include "G4AutoLock.hh"                          
 40                                                   
 41 namespace{                                        
 42   G4Mutex myPurgMagTabulatedField3DLock = G4MU    
 43 }                                                 
 44                                                   
 45 using namespace std;                              
 46                                                   
 47 PurgMagTabulatedField3D::PurgMagTabulatedField    
 48              double zOffset )                     
 49   :fZoffset(zOffset),invertX(false),invertY(fa    
 50 {                                                 
 51                                                   
 52   double lenUnit= meter;                          
 53   double fieldUnit= tesla;                        
 54   G4cout << "\n-------------------------------    
 55    << "\n      Magnetic field"                    
 56    << "\n-------------------------------------    
 57                                                   
 58   G4cout << "\n ---> " "Reading the field grid    
 59                                                   
 60   //                                              
 61   //This is a thread-local class and we have t    
 62   //file at the same time                         
 63   G4AutoLock lock(&myPurgMagTabulatedField3DLo    
 64                                                   
 65   ifstream file( filename ); // Open the file     
 66                                                   
 67   if (!file.is_open())                            
 68     {                                             
 69       G4ExceptionDescription ed;                  
 70       ed << "Could not open input file " << fi    
 71       G4Exception("PurgMagTabulatedField3D::Pu    
 72       "pugmag001",FatalException,ed);             
 73     }                                             
 74                                                   
 75   // Ignore first blank line                      
 76   char buffer[256];                               
 77   file.getline(buffer,256);                       
 78                                                   
 79   // Read table dimensions                        
 80   file >> nx >> ny >> nz; // Note dodgy order     
 81                                                   
 82   G4cout << "  [ Number of values x,y,z: "        
 83    << nx << " " << ny << " " << nz << " ] "       
 84    << G4endl;                                     
 85                                                   
 86   // Set up storage space for table               
 87   xField.resize( nx );                            
 88   yField.resize( nx );                            
 89   zField.resize( nx );                            
 90   int ix, iy, iz;                                 
 91   for (ix=0; ix<nx; ix++) {                       
 92     xField[ix].resize(ny);                        
 93     yField[ix].resize(ny);                        
 94     zField[ix].resize(ny);                        
 95     for (iy=0; iy<ny; iy++) {                     
 96       xField[ix][iy].resize(nz);                  
 97       yField[ix][iy].resize(nz);                  
 98       zField[ix][iy].resize(nz);                  
 99     }                                             
100   }                                               
101                                                   
102   // Ignore other header information              
103   // The first line whose second character is     
104   // be the last line of the header.              
105   do {                                            
106     file.getline(buffer,256);                     
107   } while ( buffer[1]!='0');                      
108                                                   
109   // Read in the data                             
110   double xval,yval,zval,bx,by,bz;                 
111   double permeability; // Not used in this exa    
112   for (ix=0; ix<nx; ix++) {                       
113     for (iy=0; iy<ny; iy++) {                     
114       for (iz=0; iz<nz; iz++) {                   
115         file >> xval >> yval >> zval >> bx >>     
116         if ( ix==0 && iy==0 && iz==0 ) {          
117           minx = xval * lenUnit;                  
118           miny = yval * lenUnit;                  
119           minz = zval * lenUnit;                  
120         }                                         
121         xField[ix][iy][iz] = bx * fieldUnit;      
122         yField[ix][iy][iz] = by * fieldUnit;      
123         zField[ix][iy][iz] = bz * fieldUnit;      
124       }                                           
125     }                                             
126   }                                               
127   file.close();                                   
128                                                   
129   lock.unlock();                                  
130                                                   
131   maxx = xval * lenUnit;                          
132   maxy = yval * lenUnit;                          
133   maxz = zval * lenUnit;                          
134                                                   
135   G4cout << "\n ---> ... done reading " << G4e    
136                                                   
137   // G4cout << " Read values of field from fil    
138   G4cout << " ---> assumed the order:  x, y, z    
139    << "\n ---> Min values x,y,z: "                
140    << minx/cm << " " << miny/cm << " " << minz    
141    << "\n ---> Max values x,y,z: "                
142    << maxx/cm << " " << maxy/cm << " " << maxz    
143    << "\n ---> The field will be offset by " <    
144                                                   
145   // Should really check that the limits are n    
146   if (maxx < minx) {swap(maxx,minx); invertX =    
147   if (maxy < miny) {swap(maxy,miny); invertY =    
148   if (maxz < minz) {swap(maxz,minz); invertZ =    
149   G4cout << "\nAfter reordering if neccesary"     
150    << "\n ---> Min values x,y,z: "                
151    << minx/cm << " " << miny/cm << " " << minz    
152    << " \n ---> Max values x,y,z: "               
153    << maxx/cm << " " << maxy/cm << " " << maxz    
154                                                   
155   dx = maxx - minx;                               
156   dy = maxy - miny;                               
157   dz = maxz - minz;                               
158   G4cout << "\n ---> Dif values x,y,z (range):    
159    << dx/cm << " " << dy/cm << " " << dz/cm <<    
160    << "\n-------------------------------------    
161 }                                                 
162                                                   
163 void PurgMagTabulatedField3D::GetFieldValue(co    
164               double *Bfield ) const              
165 {                                                 
166                                                   
167   double x = point[0];                            
168   double y = point[1];                            
169   double z = point[2] + fZoffset;                 
170                                                   
171   // Check that the point is within the define    
172   if ( x>=minx && x<=maxx &&                      
173        y>=miny && y<=maxy &&                      
174        z>=minz && z<=maxz ) {                     
175                                                   
176     // Position of given point within region,     
177     // [0,1]                                      
178     double xfraction = (x - minx) / dx;           
179     double yfraction = (y - miny) / dy;           
180     double zfraction = (z - minz) / dz;           
181                                                   
182     if (invertX) { xfraction = 1 - xfraction;}    
183     if (invertY) { yfraction = 1 - yfraction;}    
184     if (invertZ) { zfraction = 1 - zfraction;}    
185                                                   
186     // Need addresses of these to pass to modf    
187     // modf uses its second argument as an OUT    
188     double xdindex, ydindex, zdindex;             
189                                                   
190     // Position of the point within the cuboid    
191     // nearest surrounding tabulated points       
192     double xlocal = ( std::modf(xfraction*(nx-    
193     double ylocal = ( std::modf(yfraction*(ny-    
194     double zlocal = ( std::modf(zfraction*(nz-    
195                                                   
196     // The indices of the nearest tabulated po    
197     // are all less than those of the given po    
198     int xindex = static_cast<int>(xdindex);       
199     int yindex = static_cast<int>(ydindex);       
200     int zindex = static_cast<int>(zdindex);       
201                                                   
202                                                   
203 #ifdef DEBUG_INTERPOLATING_FIELD                  
204     G4cout << "Local x,y,z: " << xlocal << " "    
205     G4cout << "Index x,y,z: " << xindex << " "    
206     double valx0z0, mulx0z0, valx1z0, mulx1z0;    
207     double valx0z1, mulx0z1, valx1z1, mulx1z1;    
208     valx0z0= table[xindex  ][0][zindex];  mulx    
209     valx1z0= table[xindex+1][0][zindex];  mulx    
210     valx0z1= table[xindex  ][0][zindex+1]; mul    
211     valx1z1= table[xindex+1][0][zindex+1]; mul    
212 #endif                                            
213                                                   
214         // Full 3-dimensional version             
215     Bfield[0] =                                   
216       xField[xindex  ][yindex  ][zindex  ] * (    
217       xField[xindex  ][yindex  ][zindex+1] * (    
218       xField[xindex  ][yindex+1][zindex  ] * (    
219       xField[xindex  ][yindex+1][zindex+1] * (    
220       xField[xindex+1][yindex  ][zindex  ] *      
221       xField[xindex+1][yindex  ][zindex+1] *      
222       xField[xindex+1][yindex+1][zindex  ] *      
223       xField[xindex+1][yindex+1][zindex+1] *      
224     Bfield[1] =                                   
225       yField[xindex  ][yindex  ][zindex  ] * (    
226       yField[xindex  ][yindex  ][zindex+1] * (    
227       yField[xindex  ][yindex+1][zindex  ] * (    
228       yField[xindex  ][yindex+1][zindex+1] * (    
229       yField[xindex+1][yindex  ][zindex  ] *      
230       yField[xindex+1][yindex  ][zindex+1] *      
231       yField[xindex+1][yindex+1][zindex  ] *      
232       yField[xindex+1][yindex+1][zindex+1] *      
233     Bfield[2] =                                   
234       zField[xindex  ][yindex  ][zindex  ] * (    
235       zField[xindex  ][yindex  ][zindex+1] * (    
236       zField[xindex  ][yindex+1][zindex  ] * (    
237       zField[xindex  ][yindex+1][zindex+1] * (    
238       zField[xindex+1][yindex  ][zindex  ] *      
239       zField[xindex+1][yindex  ][zindex+1] *      
240       zField[xindex+1][yindex+1][zindex  ] *      
241       zField[xindex+1][yindex+1][zindex+1] *      
242                                                   
243   } else {                                        
244     Bfield[0] = 0.0;                              
245     Bfield[1] = 0.0;                              
246     Bfield[2] = 0.0;                              
247   }                                               
248 }                                                 
249                                                   
250