Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4ConstRK4.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 /geometry/magneticfield/src/G4ConstRK4.cc (Version 11.3.0) and /geometry/magneticfield/src/G4ConstRK4.cc (Version 9.4.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4ConstRK4 implementation                   << 
 27 //                                                 26 //
 28 // Created: J.Apostolakis, T.Nikitina - 18.09. <<  27 // $Id: G4ConstRK4.cc,v 1.5 2010-09-10 15:51:10 japost Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-04-patch-01 $
                                                   >>  29 //
                                                   >>  30 //
                                                   >>  31 // - 18.09.2008 - J.Apostolakis, T.Nikitina - Created
 29 // -------------------------------------------     32 // -------------------------------------------------------------------
 30                                                    33 
 31 #include "G4ConstRK4.hh"                           34 #include "G4ConstRK4.hh"
 32 #include "G4ThreeVector.hh"                        35 #include "G4ThreeVector.hh"
 33 #include "G4LineSection.hh"                        36 #include "G4LineSection.hh"
 34                                                    37 
 35 //////////////////////////////////////////////     38 //////////////////////////////////////////////////////////////////
 36 //                                                 39 //
 37 // Constructor sets the number of *State* vari     40 // Constructor sets the number of *State* variables (default = 8)
 38 //   The number of variables integrated is alw     41 //   The number of variables integrated is always 6
 39 //                                             <<  42 
 40 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4i     43 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numStateVariables)
 41   : G4MagErrorStepper(EqRhs, 6, numStateVariab     44   : G4MagErrorStepper(EqRhs, 6, numStateVariables)
 42 {                                                  45 {
 43   // const G4int numberOfVariables= 6;             46   // const G4int numberOfVariables= 6;
 44   if( numStateVariables < 8 )                      47   if( numStateVariables < 8 ) 
 45   {                                                48   {
 46     std::ostringstream message;                <<  49     G4cerr << "ERROR in G4ConstRK4::G4ConstRK4 " 
 47     message << "The number of State variables  <<  50      << "   The number of State variables at least 8 " << G4endl;
 48             << "Instead it is - numStateVariab <<  51     G4cerr << "   Instead it is  numStateVariables= " << numStateVariables << G4endl; 
 49     G4Exception("G4ConstRK4::G4ConstRK4()", "G <<  52     G4Exception("G4ConstRK4::G4ConstRK4()", "InvalidSetup", FatalException,
 50                 FatalException, message, "Use  <<  53   "Valid only for number of state variables of 8 or more. Use another Stepper!");
 51   }                                                54   }
 52                                                    55 
 53   fEq = EqRhs;                                     56   fEq = EqRhs;
 54   yMiddle  = new G4double[8];                      57   yMiddle  = new G4double[8];
 55   dydxMid  = new G4double[8];                      58   dydxMid  = new G4double[8];
 56   yInitial = new G4double[8];                      59   yInitial = new G4double[8];
 57   yOneStep = new G4double[8];                      60   yOneStep = new G4double[8];
 58                                                    61 
 59   dydxm = new G4double[8];                         62   dydxm = new G4double[8];
 60   dydxt = new G4double[8];                         63   dydxt = new G4double[8]; 
 61   yt    = new G4double[8];                         64   yt    = new G4double[8]; 
 62   Field[0]=0.; Field[1]=0.; Field[2]=0.;           65   Field[0]=0.; Field[1]=0.; Field[2]=0.;
 63 }                                                  66 }
 64                                                    67 
 65 //////////////////////////////////////////////     68 ////////////////////////////////////////////////////////////////
 66 //                                                 69 //
 67 // Destructor                                      70 // Destructor
 68                                                    71 
 69 G4ConstRK4::~G4ConstRK4()                          72 G4ConstRK4::~G4ConstRK4()
 70 {                                                  73 {
 71    delete [] yMiddle;                              74    delete [] yMiddle;
 72    delete [] dydxMid;                              75    delete [] dydxMid;
 73    delete [] yInitial;                             76    delete [] yInitial;
 74    delete [] yOneStep;                             77    delete [] yOneStep;
 75    delete [] dydxm;                                78    delete [] dydxm;
 76    delete [] dydxt;                                79    delete [] dydxt;
 77    delete [] yt;                                   80    delete [] yt;
 78 }                                                  81 }
 79                                                    82 
 80 //////////////////////////////////////////////     83 //////////////////////////////////////////////////////////////////////
 81 //                                                 84 //
 82 // Given values for the variables y[0,..,n-1]      85 // Given values for the variables y[0,..,n-1] and their derivatives
 83 // dydx[0,...,n-1] known at x, use the classic     86 // dydx[0,...,n-1] known at x, use the classical 4th Runge-Kutta
 84 // method to advance the solution over an inte     87 // method to advance the solution over an interval h and return the
 85 // incremented variables as yout[0,...,n-1], w     88 // incremented variables as yout[0,...,n-1], which is not a distinct
 86 // array from y. The user supplies the routine     89 // array from y. The user supplies the routine RightHandSide(x,y,dydx),
 87 // which returns derivatives dydx at x. The so     90 // which returns derivatives dydx at x. The source is routine rk4 from
 88 // NRC p. 712-713 .                                91 // NRC p. 712-713 .
 89 //                                             <<  92 
 90 void G4ConstRK4::DumbStepper( const G4double y <<  93 void G4ConstRK4::DumbStepper( const G4double  yIn[],
 91                               const G4double d <<  94                               const G4double  dydx[],
 92                                     G4double h <<  95                                     G4double  h,
 93                                     G4double y <<  96                                     G4double  yOut[])
 94 {                                                  97 {
 95    G4double  hh = h*0.5 , h6 = h/6.0  ;            98    G4double  hh = h*0.5 , h6 = h/6.0  ;
 96                                                    99    
 97    // 1st Step K1=h*dydx                          100    // 1st Step K1=h*dydx
 98    yt[5] = yIn[5] + hh*dydx[5] ;                  101    yt[5] = yIn[5] + hh*dydx[5] ;
 99    yt[4] = yIn[4] + hh*dydx[4] ;                  102    yt[4] = yIn[4] + hh*dydx[4] ;
100    yt[3] = yIn[3] + hh*dydx[3] ;                  103    yt[3] = yIn[3] + hh*dydx[3] ;
101    yt[2] = yIn[2] + hh*dydx[2] ;                  104    yt[2] = yIn[2] + hh*dydx[2] ;
102    yt[1] = yIn[1] + hh*dydx[1] ;                  105    yt[1] = yIn[1] + hh*dydx[1] ;
103    yt[0] = yIn[0] + hh*dydx[0] ;                  106    yt[0] = yIn[0] + hh*dydx[0] ;
104    RightHandSideConst(yt,dydxt) ;                 107    RightHandSideConst(yt,dydxt) ;        
105                                                   108 
106    // 2nd Step K2=h*dydxt                         109    // 2nd Step K2=h*dydxt
107    yt[5] = yIn[5] + hh*dydxt[5] ;                 110    yt[5] = yIn[5] + hh*dydxt[5] ;
108    yt[4] = yIn[4] + hh*dydxt[4] ;                 111    yt[4] = yIn[4] + hh*dydxt[4] ;
109    yt[3] = yIn[3] + hh*dydxt[3] ;                 112    yt[3] = yIn[3] + hh*dydxt[3] ;
110    yt[2] = yIn[2] + hh*dydxt[2] ;                 113    yt[2] = yIn[2] + hh*dydxt[2] ;
111    yt[1] = yIn[1] + hh*dydxt[1] ;                 114    yt[1] = yIn[1] + hh*dydxt[1] ;
112    yt[0] = yIn[0] + hh*dydxt[0] ;                 115    yt[0] = yIn[0] + hh*dydxt[0] ;
113    RightHandSideConst(yt,dydxm) ;                 116    RightHandSideConst(yt,dydxm) ;     
114                                                   117 
115    // 3rd Step K3=h*dydxm                         118    // 3rd Step K3=h*dydxm
116    // now dydxm=(K2+K3)/h                         119    // now dydxm=(K2+K3)/h
117    yt[5] = yIn[5] + h*dydxm[5] ;                  120    yt[5] = yIn[5] + h*dydxm[5] ;
118    dydxm[5] += dydxt[5] ;                         121    dydxm[5] += dydxt[5] ;  
119    yt[4] = yIn[4] + h*dydxm[4] ;                  122    yt[4] = yIn[4] + h*dydxm[4] ;
120    dydxm[4] += dydxt[4] ;                         123    dydxm[4] += dydxt[4] ;  
121    yt[3] = yIn[3] + h*dydxm[3] ;                  124    yt[3] = yIn[3] + h*dydxm[3] ;
122    dydxm[3] += dydxt[3] ;                         125    dydxm[3] += dydxt[3] ;  
123    yt[2] = yIn[2] + h*dydxm[2] ;                  126    yt[2] = yIn[2] + h*dydxm[2] ;
124    dydxm[2] += dydxt[2] ;                         127    dydxm[2] += dydxt[2] ;  
125    yt[1] = yIn[1] + h*dydxm[1] ;                  128    yt[1] = yIn[1] + h*dydxm[1] ;
126    dydxm[1] += dydxt[1] ;                         129    dydxm[1] += dydxt[1] ;  
127    yt[0] = yIn[0] + h*dydxm[0] ;                  130    yt[0] = yIn[0] + h*dydxm[0] ;
128    dydxm[0] += dydxt[0] ;                         131    dydxm[0] += dydxt[0] ;  
129    RightHandSideConst(yt,dydxt) ;                 132    RightHandSideConst(yt,dydxt) ;   
130                                                   133 
131    // 4th Step K4=h*dydxt                         134    // 4th Step K4=h*dydxt
132    yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*d    135    yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*dydxm[5]);
133    yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*d    136    yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*dydxm[4]);
134    yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*d    137    yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*dydxm[3]);
135    yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*d    138    yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*dydxm[2]);
136    yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*d    139    yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*dydxm[1]);
137    yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*d    140    yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*dydxm[0]);
138                                                   141    
139 }  // end of DumbStepper .....................    142 }  // end of DumbStepper ....................................................
140                                                   143 
141 //////////////////////////////////////////////    144 ////////////////////////////////////////////////////////////////
142 //                                                145 //
143 // Stepper                                        146 // Stepper
144                                                   147 
145 void                                              148 void
146 G4ConstRK4::Stepper( const G4double yInput[],     149 G4ConstRK4::Stepper( const G4double yInput[],
147                      const G4double dydx[],       150                      const G4double dydx[],
148                            G4double hstep,        151                            G4double hstep,
149                            G4double yOutput[],    152                            G4double yOutput[],
150                            G4double yError []     153                            G4double yError [] )
151 {                                                 154 {
152    const G4int nvar = 6;  // number of variabl    155    const G4int nvar = 6;  // number of variables integrated
153    const G4int maxvar = GetNumberOfStateVariab << 156    const G4int maxvar= GetNumberOfStateVariables();
154                                                   157 
155    // Correction for Richardson extrapolation     158    // Correction for Richardson extrapolation
156    G4double  correction = 1. / ( (1 << Integra    159    G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
157                                                   160 
158    G4int i;                                       161    G4int i;
159                                                   162    
160    // Saving yInput because yInput and yOutput    163    // Saving yInput because yInput and yOutput can be aliases for same array
161    for (i=0;    i<maxvar; ++i) { yInitial[i]=  << 164    for (i=0;    i<maxvar; i++) { yInitial[i]= yInput[i]; }
162                                                   165  
163    // Must copy the part of the state *not* in    166    // Must copy the part of the state *not* integrated to the output
164    for (i=nvar; i<maxvar; ++i) { yOutput[i]=   << 167    for (i=nvar; i<maxvar; i++) { yOutput[i]=  yInput[i]; }
165                                                   168 
166    // yInitial[7]= yInput[7];  //  The time is    169    // yInitial[7]= yInput[7];  //  The time is typically needed
167    yMiddle[7]  = yInput[7];   // Copy the time    170    yMiddle[7]  = yInput[7];   // Copy the time from initial value 
168    yOneStep[7] = yInput[7];   // As it contrib    171    yOneStep[7] = yInput[7];   // As it contributes to final value of yOutput ?
169    // yOutput[7] = yInput[7];  // -> dumb step    172    // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
170    yError[7] = 0.0;                               173    yError[7] = 0.0;         
171                                                   174 
172    G4double halfStep = hstep * 0.5;               175    G4double halfStep = hstep * 0.5; 
173                                                   176 
174    // Do two half steps                           177    // Do two half steps
175    //                                             178    //
176    GetConstField(yInitial,Field);                 179    GetConstField(yInitial,Field);
177    DumbStepper  (yInitial,  dydx,   halfStep,     180    DumbStepper  (yInitial,  dydx,   halfStep, yMiddle);
178    RightHandSideConst(yMiddle, dydxMid);          181    RightHandSideConst(yMiddle, dydxMid);    
179    DumbStepper  (yMiddle, dydxMid, halfStep, y    182    DumbStepper  (yMiddle, dydxMid, halfStep, yOutput); 
180                                                   183 
181    // Store midpoint, chord calculation           184    // Store midpoint, chord calculation
182    //                                             185    //
183    fMidPoint = G4ThreeVector( yMiddle[0],  yMi    186    fMidPoint = G4ThreeVector( yMiddle[0],  yMiddle[1],  yMiddle[2]); 
184                                                   187 
185    // Do a full Step                              188    // Do a full Step
186    //                                             189    //
187    DumbStepper(yInitial, dydx, hstep, yOneStep    190    DumbStepper(yInitial, dydx, hstep, yOneStep);
188    for(i=0; i<nvar; ++i)                       << 191    for(i=0;i<nvar;i++)
189    {                                              192    {
190       yError [i] = yOutput[i] - yOneStep[i] ;     193       yError [i] = yOutput[i] - yOneStep[i] ;
191       yOutput[i] += yError[i]*correction ;        194       yOutput[i] += yError[i]*correction ;
192         // Provides accuracy increased by 1 or    195         // Provides accuracy increased by 1 order via the 
193         // Richardson extrapolation               196         // Richardson extrapolation  
194    }                                              197    }
195                                                   198 
196    fInitialPoint = G4ThreeVector( yInitial[0],    199    fInitialPoint = G4ThreeVector( yInitial[0], yInitial[1], yInitial[2]); 
197    fFinalPoint   = G4ThreeVector( yOutput[0],     200    fFinalPoint   = G4ThreeVector( yOutput[0],  yOutput[1],  yOutput[2]); 
198                                                   201 
199    return;                                        202    return;
200 }                                                 203 }
201                                                   204 
202 //////////////////////////////////////////////    205 ////////////////////////////////////////////////////////////////
203 //                                                206 //
204 // Estimate the maximum distance from the curv    207 // Estimate the maximum distance from the curve to the chord
205 //                                                208 //
206 // We estimate this using the distance of the     209 // We estimate this using the distance of the midpoint to chord.
207 // The method below is good only for angle dev    210 // The method below is good only for angle deviations < 2 pi;
208 // this restriction should not be a problem fo    211 // this restriction should not be a problem for the Runge Kutta methods, 
209 // which generally cannot integrate accurately    212 // which generally cannot integrate accurately for large angle deviations
210 //                                             << 213 
211 G4double G4ConstRK4::DistChord() const            214 G4double G4ConstRK4::DistChord() const 
212 {                                                 215 {
213   G4double distLine, distChord;                   216   G4double distLine, distChord; 
214                                                   217 
215   if (fInitialPoint != fFinalPoint)               218   if (fInitialPoint != fFinalPoint)
216   {                                               219   {
217      distLine= G4LineSection::Distline( fMidPo    220      distLine= G4LineSection::Distline( fMidPoint, fInitialPoint, fFinalPoint );
218        // This is a class method that gives di    221        // This is a class method that gives distance of Mid 
219        // from the Chord between the Initial a    222        // from the Chord between the Initial and Final points
220      distChord = distLine;                        223      distChord = distLine;
221   }                                               224   }
222   else                                            225   else
223   {                                               226   {
224      distChord = (fMidPoint-fInitialPoint).mag    227      distChord = (fMidPoint-fInitialPoint).mag();
225   }                                               228   }
226   return distChord;                               229   return distChord;
227 }                                                 230 }
228                                                   231