| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
5 // * The Geant4 software is copyright of th << 7 //
6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4ImplicitEuler.cc,v 1.3 2000/11/01 15:15:53 gcosmo Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-03-01 $
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 // 10 //
26 // G4ImplicitEuler implementation <<
27 // 11 //
28 // Implicit Euler: 12 // Implicit Euler:
29 // 13 //
30 // x_1 = x_0 + h/2 * ( dx(t_0,x_0) + dx 14 // x_1 = x_0 + h/2 * ( dx(t_0,x_0) + dx(t_0+h,x_0+h*dx(t_0,x_0) ) )
31 // 15 //
32 // Second order solver. << 16 // second order solver
33 // Take the current derivative and add it to t 17 // Take the current derivative and add it to the current position.
34 // Take the output and its derivative. Add the 18 // Take the output and its derivative. Add the mean of both derivatives
35 // to form the final output. << 19 // to form the final output
>> 20 //
>> 21 // W.Wander <wwc@mit.edu> 12/09/97
>> 22 // 6.11.98 V.Grichine, new constructor, fNumberOfVariables
36 // 23 //
37 // Author: W. Wander <wwc@mit.edu>, 12.09.1997 <<
38 // ------------------------------------------- <<
39 24
40 #include "G4ImplicitEuler.hh" 25 #include "G4ImplicitEuler.hh"
41 #include "G4ThreeVector.hh" 26 #include "G4ThreeVector.hh"
42 27
43 ////////////////////////////////////////////// 28 /////////////////////////////////////////////////////////////////////////
44 // 29 //
45 // Constructor 30 // Constructor
46 31
47 G4ImplicitEuler::G4ImplicitEuler(G4EquationOfM << 32 G4ImplicitEuler::G4ImplicitEuler(G4Mag_EqRhs *EqRhs,
48 G4int numberO << 33 G4int numberOfVariables):
49 : G4MagErrorStepper(EqRhs, numberOfVariables << 34 G4MagErrorStepper(EqRhs, numberOfVariables),
>> 35 fNumberOfVariables(numberOfVariables)
50 { 36 {
51 unsigned int noVariables = std::max(numberOf <<
52 dydxTemp = new G4double[noVariables] ; <<
53 yTemp = new G4double[noVariables] ; <<
54 } 37 }
55 38
56 39
57 ////////////////////////////////////////////// 40 ////////////////////////////////////////////////////////////////////////
58 // 41 //
59 // Destructor 42 // Destructor
60 // << 43
61 G4ImplicitEuler::~G4ImplicitEuler() 44 G4ImplicitEuler::~G4ImplicitEuler()
62 { 45 {
63 delete [] dydxTemp; <<
64 delete [] yTemp; <<
65 } 46 }
66 47
67 ////////////////////////////////////////////// 48 //////////////////////////////////////////////////////////////////////
68 // 49 //
69 // DumbStepper <<
70 // 50 //
>> 51
71 void 52 void
72 G4ImplicitEuler::DumbStepper( const G4double y << 53 G4ImplicitEuler::DumbStepper( const G4double yIn[],
73 const G4double d << 54 const G4double dydx[],
74 G4double h << 55 G4double h,
75 G4double y << 56 G4double yOut[])
76 { 57 {
77 const G4int numberOfVariables = GetNumberOfV << 58 // const G4int nvar = 6 ;
>> 59 G4double* dydxTemp = new G4double[fNumberOfVariables] ;
>> 60 G4double* yTemp = new G4double[fNumberOfVariables] ;
78 61
79 // Initialise time to t0, needed when it is << 62 G4int i;
80 // <<
81 yTemp[7] = yOut[7] = yIn[7]; // Better to <<
82 63
83 for( G4int i = 0; i < numberOfVariables; ++i << 64 for( i = 0; i < fNumberOfVariables; i++ )
84 { 65 {
85 yTemp[i] = yIn[i] + h*dydx[i] ; 66 yTemp[i] = yIn[i] + h*dydx[i] ;
86 } 67 }
87 68
88 RightHandSide(yTemp,dydxTemp); 69 RightHandSide(yTemp,dydxTemp);
89 70
90 for( G4int i = 0; i < numberOfVariables; ++i << 71 for( i = 0; i < fNumberOfVariables; i++ )
91 { 72 {
92 yOut[i] = yIn[i] + 0.5 * h * ( dydx[i] + d 73 yOut[i] = yIn[i] + 0.5 * h * ( dydx[i] + dydxTemp[i] );
93 } 74 }
94 75
95 return; << 76 // NormaliseTangentVector( yOut );
>> 77
>> 78 return ;
96 } 79 }
97 80