| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4CashKarpRKF45 implementation << 23 //
>> 24 // $Id: G4CashKarpRKF45.cc,v 1.8.2.1 2001/06/28 19:08:18 gunter Exp $
>> 25 // GEANT4 tag $Name: $
27 // 26 //
28 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 meth 27 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
29 // order method (giving fifth-order accuracy) << 28 // order method (giving fifth-order accuracy) for the solution
30 // Two different fourth order estimates are ca << 29 // of an ODE. Two different fourth order estimates are calculated;
31 // gives an error estimate. [ref. Numerical Re << 30 // their difference gives an error estimate.
32 // It is used to integrate the equations of th << 31 // (We use it to integrate the equations of the motion of a particle
33 // in a magnetic field. << 32 // in a magnetic field. )
34 // 33 //
35 // [ref. Numerical Recipes in C, 2nd Edition] << 34 // Similar to Numerical Recipes, .... put REFerence here!
36 // 35 //
37 // Authors: J.Apostolakis, V.Grichine - 30.01. <<
38 // ------------------------------------------- <<
39 36
40 #include "G4CashKarpRKF45.hh" 37 #include "G4CashKarpRKF45.hh"
41 #include "G4LineSection.hh" 38 #include "G4LineSection.hh"
42 39
43 ////////////////////////////////////////////// 40 /////////////////////////////////////////////////////////////////////
44 // 41 //
45 // Constructor 42 // Constructor
46 // << 43
47 G4CashKarpRKF45::G4CashKarpRKF45(G4EquationOfM << 44 G4CashKarpRKF45::G4CashKarpRKF45(G4EquationOfMotion *EqRhs, G4int numberOfVariables, G4bool primary)
48 G4int noInteg << 45 : G4MagIntegratorStepper(EqRhs, numberOfVariables)
49 G4bool primar <<
50 : G4MagIntegratorStepper(EqRhs, noIntegratio <<
51 { 46 {
52 const G4int numberOfVariables = << 47 fNumberOfVariables = numberOfVariables ;
53 std::max( noIntegrationVariables, << 48
54 ( ( (noIntegrationVariables-1)/ << 49 ak2 = new G4double[fNumberOfVariables] ;
55 // For better alignment with cache-line << 50 ak3 = new G4double[fNumberOfVariables] ;
56 << 51 ak4 = new G4double[fNumberOfVariables] ;
57 ak2 = new G4double[numberOfVariables] ; << 52 ak5 = new G4double[fNumberOfVariables] ;
58 ak3 = new G4double[numberOfVariables] ; << 53 ak6 = new G4double[fNumberOfVariables] ;
59 ak4 = new G4double[numberOfVariables] ; << 54 yTemp = new G4double[fNumberOfVariables] ;
60 ak5 = new G4double[numberOfVariables] ; << 55 yIn = new G4double[fNumberOfVariables] ;
61 ak6 = new G4double[numberOfVariables] ; << 56
62 // ak7 = 0; << 57 fLastInitialVector = new G4double[fNumberOfVariables] ;
63 << 58 fLastFinalVector = new G4double[fNumberOfVariables] ;
64 // Must ensure space extra 'state' variables << 59 fLastDyDx = new G4double[fNumberOfVariables];
65 const G4int numStateMax = std::max(GetNumbe << 60
66 const G4int numStateVars = std::max(noIntegr << 61 fMidVector = new G4double[fNumberOfVariables];
67 numState << 62 fMidError = new G4double[fNumberOfVariables];
68 // GetNumbe << 63 fAuxStepper = 0;
69 << 64 if( primary )
70 yTemp = new G4double[numStateVars] ; << 65 fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
71 yIn = new G4double[numStateVars] ; << 66
72 <<
73 fLastInitialVector = new G4double[numStateVa <<
74 fLastFinalVector = new G4double[numStateVars <<
75 fLastDyDx = new G4double[numberOfVariables]; <<
76 <<
77 fMidVector = new G4double[numStateVars]; <<
78 fMidError = new G4double[numStateVars]; <<
79 if( primary ) <<
80 { <<
81 fAuxStepper = new G4CashKarpRKF45(EqRhs, n <<
82 } <<
83 } 67 }
84 68
85 ////////////////////////////////////////////// 69 /////////////////////////////////////////////////////////////////////
86 // 70 //
87 // Destructor 71 // Destructor
88 // << 72
89 G4CashKarpRKF45::~G4CashKarpRKF45() 73 G4CashKarpRKF45::~G4CashKarpRKF45()
90 { 74 {
91 delete [] ak2; << 75 delete[] ak2;
92 delete [] ak3; << 76 delete[] ak3;
93 delete [] ak4; << 77 delete[] ak4;
94 delete [] ak5; << 78 delete[] ak5;
95 delete [] ak6; << 79 delete[] ak6;
96 // delete [] ak7; << 80 delete[] ak7;
97 delete [] yTemp; << 81 delete[] yTemp;
98 delete [] yIn; << 82 delete[] yIn;
99 << 83
100 delete [] fLastInitialVector; << 84 delete[] fLastInitialVector;
101 delete [] fLastFinalVector; << 85 delete[] fLastFinalVector;
102 delete [] fLastDyDx; << 86 delete[] fLastDyDx;
103 delete [] fMidVector; << 87 delete[] fMidVector;
104 delete [] fMidError; << 88 delete[] fMidError;
105 89
106 delete fAuxStepper; 90 delete fAuxStepper;
107 } 91 }
108 92
109 ////////////////////////////////////////////// 93 //////////////////////////////////////////////////////////////////////
110 // 94 //
111 // Given values for n = 6 variables yIn[0,..., 95 // Given values for n = 6 variables yIn[0,...,n-1]
112 // known at x, use the fifth-order Cash-Karp 96 // known at x, use the fifth-order Cash-Karp Runge-
113 // Kutta-Fehlberg-4-5 method to advance the so 97 // Kutta-Fehlberg-4-5 method to advance the solution over an interval
114 // Step and return the incremented variables a 98 // Step and return the incremented variables as yOut[0,...,n-1]. Also
115 // return an estimate of the local truncation 99 // return an estimate of the local truncation error yErr[] using the
116 // embedded 4th-order method. The user supplie 100 // embedded 4th-order method. The user supplies routine
117 // RightHandSide(y,dydx), which returns deriva 101 // RightHandSide(y,dydx), which returns derivatives dydx for y .
118 // << 102
119 void 103 void
120 G4CashKarpRKF45::Stepper(const G4double yInput 104 G4CashKarpRKF45::Stepper(const G4double yInput[],
121 const G4double dydx[] << 105 const G4double dydx[],
122 G4double Step, << 106 G4double Step,
123 G4double yOut[] << 107 G4double yOut[],
124 G4double yErr[] << 108 G4double yErr[])
125 { 109 {
126 // const G4int nvar = 6 ; 110 // const G4int nvar = 6 ;
127 // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 111 // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
128 G4int i; << 112 G4int i;
129 <<
130 const G4double b21 = 0.2 , <<
131 b31 = 3.0/40.0 , b32 = 9.0/4 <<
132 b41 = 0.3 , b42 = -0.9 , b43 <<
133 113
134 b51 = -11.0/54.0 , b52 = 2.5 << 114 const G4double b21 = 0.2 ,
135 b54 = 35.0/27.0 , << 115 b31 = 3.0/40.0 , b32 = 9.0/40.0 ,
>> 116 b41 = 0.3 , b42 = -0.9 , b43 = 1.2 ,
136 117
137 b61 = 1631.0/55296.0 , b62 = << 118 b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 ,
138 b63 = 575.0/13824.0 , b64 = << 119 b54 = 35.0/27.0 ,
139 b65 = 253.0/4096.0 , <<
140 120
141 c1 = 37.0/378.0 , c3 = 250.0 << 121 b61 = 1631.0/55296.0 , b62 = 175.0/512.0 ,
142 c6 = 512.0/1771.0 , dc5 = -2 << 122 b63 = 575.0/13824.0 , b64 = 44275.0/110592.0 ,
>> 123 b65 = 253.0/4096.0 ,
143 124
144 const G4double dc1 = c1 - 2825.0/27648.0 , << 125 c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 ,
145 dc4 = c4 - 13525.0/55296.0 , << 126 c6 = 512.0/1771.0 ,
>> 127 dc5 = -277.0/14336.0 ;
146 128
147 // Initialise time to t0, needed when it is << 129 const G4double dc1 = c1 - 2825.0/27648.0 , dc3 = c3 - 18575.0/48384.0 ,
148 // [ Note: Only for time dependent fi << 130 dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ;
149 // is it neccessary to inte <<
150 yOut[7] = yTemp[7] = yIn[7] = yInput[7]; <<
151 131
152 const G4int numberOfVariables= this->GetNumb <<
153 // The number of variables to be integrate <<
154 132
155 // Saving yInput because yInput and yOut ca << 133 // Saving yInput because yInput and yOut can be aliases for same array
156 134
157 for(i=0; i<numberOfVariables; ++i) << 135 for(i=0;i<fNumberOfVariables;i++)
158 { << 136 {
159 yIn[i]=yInput[i]; << 137 yIn[i]=yInput[i];
160 } << 138 }
161 // RightHandSide(yIn, dydx) ; // << 139 // RightHandSide(yIn, dydx) ; // 1st Step
162 140
163 for(i=0; i<numberOfVariables; ++i) << 141 for(i=0;i<fNumberOfVariables;i++)
164 { << 142 {
165 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; << 143 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
166 } << 144 }
167 RightHandSide(yTemp, ak2) ; // << 145 RightHandSide(yTemp, ak2) ; // 2nd Step
168 146
169 for(i=0; i<numberOfVariables; ++i) << 147 for(i=0;i<fNumberOfVariables;i++)
170 { << 148 {
171 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b3 149 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
172 } << 150 }
173 RightHandSide(yTemp, ak3) ; // << 151 RightHandSide(yTemp, ak3) ; // 3rd Step
174 152
175 for(i=0; i<numberOfVariables; ++i) << 153 for(i=0;i<fNumberOfVariables;i++)
176 { << 154 {
177 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b4 155 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
178 } << 156 }
179 RightHandSide(yTemp, ak4) ; // << 157 RightHandSide(yTemp, ak4) ; // 4th Step
180 <<
181 for(i=0; i<numberOfVariables; ++i) <<
182 { <<
183 yTemp[i] = yIn[i] + Step*(b51*dydx[i] <<
184 + b52*ak2[i] + b53*ak3[i <<
185 } <<
186 RightHandSide(yTemp, ak5) ; // <<
187 158
188 for(i=0; i<numberOfVariables; ++i) << 159 for(i=0;i<fNumberOfVariables;i++)
189 { << 160 {
190 yTemp[i] = yIn[i] + Step*(b61*dydx[i] << 161 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
191 + b62*ak2[i] + b63*ak3[i << 162 b54*ak4[i]) ;
192 } << 163 }
193 RightHandSide(yTemp, ak6) ; // << 164 RightHandSide(yTemp, ak5) ; // 5th Step
>> 165
>> 166 for(i=0;i<fNumberOfVariables;i++)
>> 167 {
>> 168 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
>> 169 b64*ak4[i] + b65*ak5[i]) ;
>> 170 }
>> 171 RightHandSide(yTemp, ak6) ; // 6th Step
194 172
195 for(i=0; i<numberOfVariables; ++i) << 173 for(i=0;i<fNumberOfVariables;i++)
196 { << 174 {
197 // Accumulate increments with proper weigh 175 // Accumulate increments with proper weights
198 // << 176
199 yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*a 177 yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ;
200 178
201 // Estimate error as difference between 4t << 179 // Estimate error as difference between 4th and
202 // << 180 // 5th order methods
203 yErr[i] = Step*(dc1*dydx[i] << 181
204 + dc3*ak3[i] + dc4*ak4[i] + dc5*ak << 182 yErr[i] = Step*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] +
>> 183 dc5*ak5[i] + dc6*ak6[i]) ;
205 184
206 // Store Input and Final values, for possi 185 // Store Input and Final values, for possible use in calculating chord
207 // <<
208 fLastInitialVector[i] = yIn[i] ; 186 fLastInitialVector[i] = yIn[i] ;
209 fLastFinalVector[i] = yOut[i]; 187 fLastFinalVector[i] = yOut[i];
210 fLastDyDx[i] = dydx[i]; 188 fLastDyDx[i] = dydx[i];
211 } << 189 }
212 // NormaliseTangentVector( yOut ); // Not wa << 190 // NormaliseTangentVector( yOut ); // Not wanted
>> 191
>> 192 fLastStepLength =Step;
213 193
214 fLastStepLength = Step; << 194 return ;
215 195
216 return; <<
217 } 196 }
218 197
219 ////////////////////////////////////////////// 198 ///////////////////////////////////////////////////////////////////////////////
220 // << 199
221 void 200 void
222 G4CashKarpRKF45::StepWithEst( const G4double*, << 201 G4CashKarpRKF45::StepWithEst(const G4double yInput[],
223 const G4double*, << 202 const G4double dydx[],
224 G4double, << 203 G4double Step,
225 G4double*, << 204 G4double yOut[],
226 G4double&, << 205 G4double& alpha2,
227 G4double&, << 206 G4double& beta2,
228 const G4double*, << 207 const G4double B1[],
229 G4double* << 208 G4double B2[] )
230 { 209 {
231 G4Exception("G4CashKarpRKF45::StepWithEst()" << 210
232 FatalException, "Method no longe << 211 G4Exception("G4CashKarpRKF45::StepWithEst ERROR: This Method is no longer used.");
233 return ; << 212
>> 213 #if 0
>> 214 // const G4int nvar = 6 ;
>> 215 G4int i;
>> 216
>> 217 // Saving yInput because yInput and yOut can be aliases for same array
>> 218
>> 219 for(i=0;i<fNumberOfVariables;i++)
>> 220 {
>> 221 yIn[i]=yInput[i];
>> 222 }
>> 223
>> 224 const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875 ,
>> 225
>> 226 b21 = 0.2 ,
>> 227 b31 = 3.0/40.0 , b32 = 9.0/40.0 ,
>> 228 b41 = 0.3 , b42 = -0.9 , b43 = 1.2 ,
>> 229
>> 230 b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 ,
>> 231 b54 = 35.0/27.0 ,
>> 232
>> 233 b61 = 1631.0/55296.0 , b62 = 175.0/512.0 ,
>> 234 b63 = 575.0/13824.0 , b64 = 44275.0/110592.0 ,
>> 235 b65 = 253.0/4096.0 ,
>> 236
>> 237 c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 ,
>> 238 c6 = 512.0/1771.0 ,
>> 239 dc5 = -277.0/14336.0 ;
>> 240
>> 241 const G4double dc1 = c1 - 2825.0/27648.0 , dc3 = c3 - 18575.0/48384.0 ,
>> 242 dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ;
>> 243
>> 244 alpha2 = 0 ;
>> 245 beta2 = 0 ;
>> 246
>> 247 // RightHandSide(yIn, dydx) ; // 1st Step
>> 248
>> 249 for(i=0;i<3;i++)
>> 250 {
>> 251 alpha2 += B1[i]*B1[i] ;
>> 252 beta2 += dydx[i+3]*dydx[i+3] ;
>> 253 }
>> 254
>> 255 for(i=0;i<fNumberOfVariables;i++)
>> 256 {
>> 257 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
>> 258 }
>> 259
>> 260 GetEquationOfMotion()->EvaluateRhsReturnB(yTemp,ak2,B2) ; // Calculates yderive &
>> 261 // returns B too!
>> 262 for(i=0;i<3;i++)
>> 263 {
>> 264 alpha2 += B2[i]*B2[i] ;
>> 265 beta2 += ak2[i+3]*ak2[i+3] ;
>> 266 }
>> 267
>> 268 for(i=0;i<fNumberOfVariables;i++)
>> 269 {
>> 270 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
>> 271 }
>> 272 GetEquationOfMotion()->EvaluateRhsReturnB(yTemp,ak3,B2) ;
>> 273 for(i=0;i<3;i++)
>> 274 {
>> 275 alpha2 += B2[i]*B2[i] ;
>> 276 beta2 += ak3[i+3]*ak3[i+3] ;
>> 277 }
>> 278
>> 279 for(i=0;i<fNumberOfVariables;i++)
>> 280 {
>> 281 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
>> 282 }
>> 283 GetEquationOfMotion()->EvaluateRhsReturnB(yTemp,ak4,B2) ;
>> 284
>> 285 for(i=0;i<3;i++)
>> 286 {
>> 287 alpha2 += B2[i]*B2[i] ;
>> 288 beta2 += ak4[i+3]*ak4[i+3] ;
>> 289 }
>> 290
>> 291 for(i=0;i<fNumberOfVariables;i++)
>> 292 {
>> 293 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
>> 294 b54*ak4[i]) ;
>> 295 }
>> 296 GetEquationOfMotion()->EvaluateRhsReturnB(yTemp,ak5,B2) ;
>> 297
>> 298 for(i=0;i<3;i++)
>> 299 {
>> 300 alpha2 += B2[i]*B2[i] ;
>> 301 beta2 += ak5[i+3]*ak5[i+3] ;
>> 302 }
>> 303
>> 304 for(i=0;i<fNumberOfVariables;i++)
>> 305 {
>> 306 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
>> 307 b64*ak4[i] + b65*ak5[i]) ;
>> 308 }
>> 309 GetEquationOfMotion()->EvaluateRhsReturnB(yTemp,ak6,B2) ;
>> 310
>> 311 for(i=0;i<3;i++)
>> 312 {
>> 313 alpha2 += B2[i]*B2[i] ;
>> 314 beta2 += ak6[i+3]*ak6[i+3] ;
>> 315 }
>> 316
>> 317 for(i=0;i<fNumberOfVariables;i++)
>> 318 {
>> 319 // Accumulate increments with proper weights
>> 320
>> 321 yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ;
>> 322 }
>> 323
>> 324 alpha2 *= sqr(GetEquationOfMotion()->FCof()*Step)/6.0 ;
>> 325 beta2 *= (Step*Step)/6.0 ;
>> 326 // NormaliseTangentVector( yOut );
>> 327 #endif
>> 328
>> 329 return ;
>> 330
234 } 331 }
235 332
236 ////////////////////////////////////////////// 333 /////////////////////////////////////////////////////////////////
237 // << 334
238 G4double G4CashKarpRKF45::DistChord() const 335 G4double G4CashKarpRKF45::DistChord() const
239 { 336 {
240 G4double distLine, distChord; 337 G4double distLine, distChord;
241 G4ThreeVector initialPoint, finalPoint, midP 338 G4ThreeVector initialPoint, finalPoint, midPoint;
242 339
243 // Store last initial and final points << 340 // Store last initial and final points (they will be overwritten in self-Stepper call!)
244 // (they will be overwritten in self-Stepper <<
245 // <<
246 initialPoint = G4ThreeVector( fLastInitialVe 341 initialPoint = G4ThreeVector( fLastInitialVector[0],
247 fLastInitialVe 342 fLastInitialVector[1], fLastInitialVector[2]);
248 finalPoint = G4ThreeVector( fLastFinalVect 343 finalPoint = G4ThreeVector( fLastFinalVector[0],
249 fLastFinalVect 344 fLastFinalVector[1], fLastFinalVector[2]);
250 345
251 // Do half a step using StepNoErr 346 // Do half a step using StepNoErr
252 // << 347
253 fAuxStepper->Stepper( fLastInitialVector, fL << 348 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
254 0.5 * fLastStepLength, << 349 fMidVector, fMidError );
255 350
256 midPoint = G4ThreeVector( fMidVector[0], fMi 351 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
257 352
258 // Use stored values of Initial and Endpoint 353 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
259 // distance of Chord << 354 // distance of Chord
260 // << 355
>> 356
261 if (initialPoint != finalPoint) 357 if (initialPoint != finalPoint)
262 { 358 {
263 distLine = G4LineSection::Distline( midP 359 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
264 distChord = distLine; 360 distChord = distLine;
265 } 361 }
266 else 362 else
267 { 363 {
268 distChord = (midPoint-initialPoint).mag() 364 distChord = (midPoint-initialPoint).mag();
269 } 365 }
270 return distChord; 366 return distChord;
271 } 367 }
>> 368
>> 369
272 370