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 10.0.p4)


  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 66356 2012-12-18 09:02:32Z gcosmo $
                                                   >>  28 //
                                                   >>  29 //
                                                   >>  30 // - 18.09.2008 - J.Apostolakis, T.Nikitina - Created
 29 // -------------------------------------------     31 // -------------------------------------------------------------------
 30                                                    32 
 31 #include "G4ConstRK4.hh"                           33 #include "G4ConstRK4.hh"
 32 #include "G4ThreeVector.hh"                        34 #include "G4ThreeVector.hh"
 33 #include "G4LineSection.hh"                        35 #include "G4LineSection.hh"
 34                                                    36 
 35 //////////////////////////////////////////////     37 //////////////////////////////////////////////////////////////////
 36 //                                                 38 //
 37 // Constructor sets the number of *State* vari     39 // Constructor sets the number of *State* variables (default = 8)
 38 //   The number of variables integrated is alw     40 //   The number of variables integrated is always 6
 39 //                                             <<  41 
 40 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4i     42 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numStateVariables)
 41   : G4MagErrorStepper(EqRhs, 6, numStateVariab     43   : G4MagErrorStepper(EqRhs, 6, numStateVariables)
 42 {                                                  44 {
 43   // const G4int numberOfVariables= 6;             45   // const G4int numberOfVariables= 6;
 44   if( numStateVariables < 8 )                      46   if( numStateVariables < 8 ) 
 45   {                                                47   {
 46     std::ostringstream message;                    48     std::ostringstream message;
 47     message << "The number of State variables      49     message << "The number of State variables at least 8 " << G4endl
 48             << "Instead it is - numStateVariab     50             << "Instead it is - numStateVariables= " << numStateVariables;
 49     G4Exception("G4ConstRK4::G4ConstRK4()", "G     51     G4Exception("G4ConstRK4::G4ConstRK4()", "GeomField0002",
 50                 FatalException, message, "Use      52                 FatalException, message, "Use another Stepper!");
 51   }                                                53   }
 52                                                    54 
 53   fEq = EqRhs;                                     55   fEq = EqRhs;
 54   yMiddle  = new G4double[8];                      56   yMiddle  = new G4double[8];
 55   dydxMid  = new G4double[8];                      57   dydxMid  = new G4double[8];
 56   yInitial = new G4double[8];                      58   yInitial = new G4double[8];
 57   yOneStep = new G4double[8];                      59   yOneStep = new G4double[8];
 58                                                    60 
 59   dydxm = new G4double[8];                         61   dydxm = new G4double[8];
 60   dydxt = new G4double[8];                         62   dydxt = new G4double[8]; 
 61   yt    = new G4double[8];                         63   yt    = new G4double[8]; 
 62   Field[0]=0.; Field[1]=0.; Field[2]=0.;           64   Field[0]=0.; Field[1]=0.; Field[2]=0.;
 63 }                                                  65 }
 64                                                    66 
 65 //////////////////////////////////////////////     67 ////////////////////////////////////////////////////////////////
 66 //                                                 68 //
 67 // Destructor                                      69 // Destructor
 68                                                    70 
 69 G4ConstRK4::~G4ConstRK4()                          71 G4ConstRK4::~G4ConstRK4()
 70 {                                                  72 {
 71    delete [] yMiddle;                              73    delete [] yMiddle;
 72    delete [] dydxMid;                              74    delete [] dydxMid;
 73    delete [] yInitial;                             75    delete [] yInitial;
 74    delete [] yOneStep;                             76    delete [] yOneStep;
 75    delete [] dydxm;                                77    delete [] dydxm;
 76    delete [] dydxt;                                78    delete [] dydxt;
 77    delete [] yt;                                   79    delete [] yt;
 78 }                                                  80 }
 79                                                    81 
 80 //////////////////////////////////////////////     82 //////////////////////////////////////////////////////////////////////
 81 //                                                 83 //
 82 // Given values for the variables y[0,..,n-1]      84 // Given values for the variables y[0,..,n-1] and their derivatives
 83 // dydx[0,...,n-1] known at x, use the classic     85 // dydx[0,...,n-1] known at x, use the classical 4th Runge-Kutta
 84 // method to advance the solution over an inte     86 // method to advance the solution over an interval h and return the
 85 // incremented variables as yout[0,...,n-1], w     87 // incremented variables as yout[0,...,n-1], which is not a distinct
 86 // array from y. The user supplies the routine     88 // array from y. The user supplies the routine RightHandSide(x,y,dydx),
 87 // which returns derivatives dydx at x. The so     89 // which returns derivatives dydx at x. The source is routine rk4 from
 88 // NRC p. 712-713 .                                90 // NRC p. 712-713 .
 89 //                                             <<  91 
 90 void G4ConstRK4::DumbStepper( const G4double y <<  92 void G4ConstRK4::DumbStepper( const G4double  yIn[],
 91                               const G4double d <<  93                               const G4double  dydx[],
 92                                     G4double h <<  94                                     G4double  h,
 93                                     G4double y <<  95                                     G4double  yOut[])
 94 {                                                  96 {
 95    G4double  hh = h*0.5 , h6 = h/6.0  ;            97    G4double  hh = h*0.5 , h6 = h/6.0  ;
 96                                                    98    
 97    // 1st Step K1=h*dydx                           99    // 1st Step K1=h*dydx
 98    yt[5] = yIn[5] + hh*dydx[5] ;                  100    yt[5] = yIn[5] + hh*dydx[5] ;
 99    yt[4] = yIn[4] + hh*dydx[4] ;                  101    yt[4] = yIn[4] + hh*dydx[4] ;
100    yt[3] = yIn[3] + hh*dydx[3] ;                  102    yt[3] = yIn[3] + hh*dydx[3] ;
101    yt[2] = yIn[2] + hh*dydx[2] ;                  103    yt[2] = yIn[2] + hh*dydx[2] ;
102    yt[1] = yIn[1] + hh*dydx[1] ;                  104    yt[1] = yIn[1] + hh*dydx[1] ;
103    yt[0] = yIn[0] + hh*dydx[0] ;                  105    yt[0] = yIn[0] + hh*dydx[0] ;
104    RightHandSideConst(yt,dydxt) ;                 106    RightHandSideConst(yt,dydxt) ;        
105                                                   107 
106    // 2nd Step K2=h*dydxt                         108    // 2nd Step K2=h*dydxt
107    yt[5] = yIn[5] + hh*dydxt[5] ;                 109    yt[5] = yIn[5] + hh*dydxt[5] ;
108    yt[4] = yIn[4] + hh*dydxt[4] ;                 110    yt[4] = yIn[4] + hh*dydxt[4] ;
109    yt[3] = yIn[3] + hh*dydxt[3] ;                 111    yt[3] = yIn[3] + hh*dydxt[3] ;
110    yt[2] = yIn[2] + hh*dydxt[2] ;                 112    yt[2] = yIn[2] + hh*dydxt[2] ;
111    yt[1] = yIn[1] + hh*dydxt[1] ;                 113    yt[1] = yIn[1] + hh*dydxt[1] ;
112    yt[0] = yIn[0] + hh*dydxt[0] ;                 114    yt[0] = yIn[0] + hh*dydxt[0] ;
113    RightHandSideConst(yt,dydxm) ;                 115    RightHandSideConst(yt,dydxm) ;     
114                                                   116 
115    // 3rd Step K3=h*dydxm                         117    // 3rd Step K3=h*dydxm
116    // now dydxm=(K2+K3)/h                         118    // now dydxm=(K2+K3)/h
117    yt[5] = yIn[5] + h*dydxm[5] ;                  119    yt[5] = yIn[5] + h*dydxm[5] ;
118    dydxm[5] += dydxt[5] ;                         120    dydxm[5] += dydxt[5] ;  
119    yt[4] = yIn[4] + h*dydxm[4] ;                  121    yt[4] = yIn[4] + h*dydxm[4] ;
120    dydxm[4] += dydxt[4] ;                         122    dydxm[4] += dydxt[4] ;  
121    yt[3] = yIn[3] + h*dydxm[3] ;                  123    yt[3] = yIn[3] + h*dydxm[3] ;
122    dydxm[3] += dydxt[3] ;                         124    dydxm[3] += dydxt[3] ;  
123    yt[2] = yIn[2] + h*dydxm[2] ;                  125    yt[2] = yIn[2] + h*dydxm[2] ;
124    dydxm[2] += dydxt[2] ;                         126    dydxm[2] += dydxt[2] ;  
125    yt[1] = yIn[1] + h*dydxm[1] ;                  127    yt[1] = yIn[1] + h*dydxm[1] ;
126    dydxm[1] += dydxt[1] ;                         128    dydxm[1] += dydxt[1] ;  
127    yt[0] = yIn[0] + h*dydxm[0] ;                  129    yt[0] = yIn[0] + h*dydxm[0] ;
128    dydxm[0] += dydxt[0] ;                         130    dydxm[0] += dydxt[0] ;  
129    RightHandSideConst(yt,dydxt) ;                 131    RightHandSideConst(yt,dydxt) ;   
130                                                   132 
131    // 4th Step K4=h*dydxt                         133    // 4th Step K4=h*dydxt
132    yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*d    134    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    135    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    136    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    137    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    138    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    139    yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*dydxm[0]);
138                                                   140    
139 }  // end of DumbStepper .....................    141 }  // end of DumbStepper ....................................................
140                                                   142 
141 //////////////////////////////////////////////    143 ////////////////////////////////////////////////////////////////
142 //                                                144 //
143 // Stepper                                        145 // Stepper
144                                                   146 
145 void                                              147 void
146 G4ConstRK4::Stepper( const G4double yInput[],     148 G4ConstRK4::Stepper( const G4double yInput[],
147                      const G4double dydx[],       149                      const G4double dydx[],
148                            G4double hstep,        150                            G4double hstep,
149                            G4double yOutput[],    151                            G4double yOutput[],
150                            G4double yError []     152                            G4double yError [] )
151 {                                                 153 {
152    const G4int nvar = 6;  // number of variabl    154    const G4int nvar = 6;  // number of variables integrated
153    const G4int maxvar = GetNumberOfStateVariab << 155    const G4int maxvar= GetNumberOfStateVariables();
154                                                   156 
155    // Correction for Richardson extrapolation     157    // Correction for Richardson extrapolation
156    G4double  correction = 1. / ( (1 << Integra    158    G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
157                                                   159 
158    G4int i;                                       160    G4int i;
159                                                   161    
160    // Saving yInput because yInput and yOutput    162    // Saving yInput because yInput and yOutput can be aliases for same array
161    for (i=0;    i<maxvar; ++i) { yInitial[i]=  << 163    for (i=0;    i<maxvar; i++) { yInitial[i]= yInput[i]; }
162                                                   164  
163    // Must copy the part of the state *not* in    165    // Must copy the part of the state *not* integrated to the output
164    for (i=nvar; i<maxvar; ++i) { yOutput[i]=   << 166    for (i=nvar; i<maxvar; i++) { yOutput[i]=  yInput[i]; }
165                                                   167 
166    // yInitial[7]= yInput[7];  //  The time is    168    // yInitial[7]= yInput[7];  //  The time is typically needed
167    yMiddle[7]  = yInput[7];   // Copy the time    169    yMiddle[7]  = yInput[7];   // Copy the time from initial value 
168    yOneStep[7] = yInput[7];   // As it contrib    170    yOneStep[7] = yInput[7];   // As it contributes to final value of yOutput ?
169    // yOutput[7] = yInput[7];  // -> dumb step    171    // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
170    yError[7] = 0.0;                               172    yError[7] = 0.0;         
171                                                   173 
172    G4double halfStep = hstep * 0.5;               174    G4double halfStep = hstep * 0.5; 
173                                                   175 
174    // Do two half steps                           176    // Do two half steps
175    //                                             177    //
176    GetConstField(yInitial,Field);                 178    GetConstField(yInitial,Field);
177    DumbStepper  (yInitial,  dydx,   halfStep,     179    DumbStepper  (yInitial,  dydx,   halfStep, yMiddle);
178    RightHandSideConst(yMiddle, dydxMid);          180    RightHandSideConst(yMiddle, dydxMid);    
179    DumbStepper  (yMiddle, dydxMid, halfStep, y    181    DumbStepper  (yMiddle, dydxMid, halfStep, yOutput); 
180                                                   182 
181    // Store midpoint, chord calculation           183    // Store midpoint, chord calculation
182    //                                             184    //
183    fMidPoint = G4ThreeVector( yMiddle[0],  yMi    185    fMidPoint = G4ThreeVector( yMiddle[0],  yMiddle[1],  yMiddle[2]); 
184                                                   186 
185    // Do a full Step                              187    // Do a full Step
186    //                                             188    //
187    DumbStepper(yInitial, dydx, hstep, yOneStep    189    DumbStepper(yInitial, dydx, hstep, yOneStep);
188    for(i=0; i<nvar; ++i)                       << 190    for(i=0;i<nvar;i++)
189    {                                              191    {
190       yError [i] = yOutput[i] - yOneStep[i] ;     192       yError [i] = yOutput[i] - yOneStep[i] ;
191       yOutput[i] += yError[i]*correction ;        193       yOutput[i] += yError[i]*correction ;
192         // Provides accuracy increased by 1 or    194         // Provides accuracy increased by 1 order via the 
193         // Richardson extrapolation               195         // Richardson extrapolation  
194    }                                              196    }
195                                                   197 
196    fInitialPoint = G4ThreeVector( yInitial[0],    198    fInitialPoint = G4ThreeVector( yInitial[0], yInitial[1], yInitial[2]); 
197    fFinalPoint   = G4ThreeVector( yOutput[0],     199    fFinalPoint   = G4ThreeVector( yOutput[0],  yOutput[1],  yOutput[2]); 
198                                                   200 
199    return;                                        201    return;
200 }                                                 202 }
201                                                   203 
202 //////////////////////////////////////////////    204 ////////////////////////////////////////////////////////////////
203 //                                                205 //
204 // Estimate the maximum distance from the curv    206 // Estimate the maximum distance from the curve to the chord
205 //                                                207 //
206 // We estimate this using the distance of the     208 // We estimate this using the distance of the midpoint to chord.
207 // The method below is good only for angle dev    209 // The method below is good only for angle deviations < 2 pi;
208 // this restriction should not be a problem fo    210 // this restriction should not be a problem for the Runge Kutta methods, 
209 // which generally cannot integrate accurately    211 // which generally cannot integrate accurately for large angle deviations
210 //                                             << 212 
211 G4double G4ConstRK4::DistChord() const            213 G4double G4ConstRK4::DistChord() const 
212 {                                                 214 {
213   G4double distLine, distChord;                   215   G4double distLine, distChord; 
214                                                   216 
215   if (fInitialPoint != fFinalPoint)               217   if (fInitialPoint != fFinalPoint)
216   {                                               218   {
217      distLine= G4LineSection::Distline( fMidPo    219      distLine= G4LineSection::Distline( fMidPoint, fInitialPoint, fFinalPoint );
218        // This is a class method that gives di    220        // This is a class method that gives distance of Mid 
219        // from the Chord between the Initial a    221        // from the Chord between the Initial and Final points
220      distChord = distLine;                        222      distChord = distLine;
221   }                                               223   }
222   else                                            224   else
223   {                                               225   {
224      distChord = (fMidPoint-fInitialPoint).mag    226      distChord = (fMidPoint-fInitialPoint).mag();
225   }                                               227   }
226   return distChord;                               228   return distChord;
227 }                                                 229 }
228                                                   230