103{
105
106
107
109 b31 = 7.0/72.0 ,
110 b32 = 7.0/24.0 ,
111
112 b41 = 3043.0/3528.0 ,
113 b42 = -3757.0/1176.0 ,
114 b43 = 1445.0/441.0,
115
116 b51 = 17617.0/11662.0 ,
117 b52 = -4023.0/686.0 ,
118 b53 = 9372.0/1715.0 ,
119 b54 = -66.0/595.0 ,
120
121 b61 = 29.0/238.0 ,
122 b62 = 0.0 ,
123 b63 = 216.0/385.0 ,
124 b64 = 54.0/85.0 ,
125 b65 = -7.0/22.0 ,
126
127 dc1 = 363.0/2975.0 - b61 ,
128 dc2 = 0.0 - b62 ,
129 dc3 = 981.0/1750.0 - b63,
130 dc4 = 2709.0/4250.0 - b64 ,
131 dc5 = -3.0/10.0 - b65 ,
132 dc6 = -1.0/50.0 ;
133
135
136
137
138 yOut[7] = yTemp[7] = yIn[7];
139
140
141
142 for(i=0; i<numberOfVariables; ++i)
143 {
144 yIn[i]=yInput[i];
145 }
146
147
148
149 for(i=0; i<numberOfVariables; ++i)
150 {
151 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
152 }
154
155 for(i=0; i<numberOfVariables; ++i)
156 {
157 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
158 }
160
161 for(i=0; i<numberOfVariables; ++i)
162 {
163 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
164 }
166
167 for(i=0; i<numberOfVariables; ++i)
168 {
169 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i]
170 + b53*ak3[i] + b54*ak4[i]) ;
171 }
173
174 for(i=0; i<numberOfVariables; ++i)
175 {
176 yOut[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i]
177 + b64*ak4[i] + b65*ak5[i]) ;
178 }
180
181 for(i=0; i<numberOfVariables; ++i)
182 {
183
184 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
185 + dc5*ak5[i] + dc6*ak6[i] ) ;
186
187
188
189 fLastInitialVector[i] = yIn[i] ;
190 fLastFinalVector[i] = yOut[i];
191 fLastDyDx[i] = DyDx[i];
192 }
193
194 fLastStepLength = Step;
195
196 return ;
197}
void RightHandSide(const G4double y[], G4double dydx[]) const