Geant4 Cross Reference |
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 // G4TsitourasRK45 implementation << 27 // 26 // 28 // Author: Somnath Banerjee, Google Summer of << 27 // Tsitouras - 5(4) RK steppers ( non-FSAL version ) 29 // Supervision: John Apostolakis, CERN << 28 // >> 29 // Implements RK tableau from 'Table 1' of >> 30 // C. Tsitouras, “Runge–Kutta pairs of order 5(4) satisfying only >> 31 // the first column simplifying assumption,” >> 32 // Computers & Mathematics with Applications, >> 33 // vol. 62, no. 2, pp. 770–775, 2011. >> 34 // >> 35 // Adaptation / Geant4 implementation by Somnath Banerjee >> 36 // Supervision / code review: John Apostolakis >> 37 // >> 38 // Sponsored by Google in Google Summer of Code 2015. >> 39 // >> 40 // First version: 12 June 2015 >> 41 // 30 // ------------------------------------------- 42 // ------------------------------------------------------------------- 31 43 32 #include "G4TsitourasRK45.hh" 44 #include "G4TsitourasRK45.hh" 33 #include "G4LineSection.hh" 45 #include "G4LineSection.hh" 34 46 35 ////////////////////////////////////////////// 47 ///////////////////////////////////////////////////////////////////// 36 // 48 // 37 // Constructor 49 // Constructor 38 // << 50 39 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfM 51 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfMotion *EqRhs, 40 G4int noInteg << 52 G4int noIntegrationVariables, 41 G4bool primar << 53 G4bool primary) 42 : G4MagIntegratorStepper(EqRhs, noIntegratio << 54 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables), >> 55 fLastStepLength(0.), fAuxStepper(0) 43 { 56 { 44 const G4int numberOfVariables = noIntegratio 57 const G4int numberOfVariables = noIntegrationVariables; >> 58 G4cout << "G4TsitourasRK45 constructor called." << G4endl; 45 59 46 ak2 = new G4double[numberOfVariables] ; 60 ak2 = new G4double[numberOfVariables] ; 47 ak3 = new G4double[numberOfVariables] ; 61 ak3 = new G4double[numberOfVariables] ; 48 ak4 = new G4double[numberOfVariables] ; 62 ak4 = new G4double[numberOfVariables] ; 49 ak5 = new G4double[numberOfVariables] ; 63 ak5 = new G4double[numberOfVariables] ; 50 ak6 = new G4double[numberOfVariables] ; 64 ak6 = new G4double[numberOfVariables] ; 51 ak7 = new G4double[numberOfVariables] ; 65 ak7 = new G4double[numberOfVariables] ; 52 ak8 = new G4double[numberOfVariables] ; 66 ak8 = new G4double[numberOfVariables] ; 53 67 >> 68 54 // Must ensure space extra 'state' variables 69 // Must ensure space extra 'state' variables exists - i.e. yIn[7] 55 // << 56 const G4int numStateMax = std::max(GetNumbe 70 const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8); 57 const G4int numStateVars = std::max(noIntegr 71 const G4int numStateVars = std::max(noIntegrationVariables, 58 numState 72 numStateMax ); >> 73 // GetNumberOfStateVariables() ); 59 74 60 yTemp = new G4double[numStateVars] ; 75 yTemp = new G4double[numStateVars] ; 61 yIn = new G4double[numStateVars] ; 76 yIn = new G4double[numStateVars] ; 62 77 63 fLastInitialVector = new G4double[numberOfVa 78 fLastInitialVector = new G4double[numberOfVariables] ; 64 fLastFinalVector = new G4double[numberOfVari 79 fLastFinalVector = new G4double[numberOfVariables] ; 65 80 66 fLastDyDx = new G4double[numberOfVariables]; 81 fLastDyDx = new G4double[numberOfVariables]; 67 82 68 fMidVector = new G4double[numberOfVariables] 83 fMidVector = new G4double[numberOfVariables]; 69 fMidError = new G4double[numberOfVariables] 84 fMidError = new G4double[numberOfVariables]; 70 << 71 if( primary ) 85 if( primary ) 72 { 86 { 73 fAuxStepper = new G4TsitourasRK45(EqRhs, n 87 fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary); 74 } 88 } 75 } 89 } 76 90 77 ////////////////////////////////////////////// 91 ///////////////////////////////////////////////////////////////////// 78 // 92 // 79 // Destructor 93 // Destructor 80 // << 94 81 G4TsitourasRK45::~G4TsitourasRK45() 95 G4TsitourasRK45::~G4TsitourasRK45() 82 { 96 { 83 delete [] ak2; << 97 delete[] ak2; 84 delete [] ak3; << 98 delete[] ak3; 85 delete [] ak4; << 99 delete[] ak4; 86 delete [] ak5; << 100 delete[] ak5; 87 delete [] ak6; << 101 delete[] ak6; 88 delete [] ak7; << 102 delete[] ak7; 89 delete [] ak8; << 103 delete[] ak8; 90 << 104 91 delete [] yTemp; << 105 delete[] yTemp; 92 delete [] yIn; << 106 delete[] yIn; 93 << 107 94 delete [] fLastInitialVector; << 108 delete[] fLastInitialVector; 95 delete [] fLastFinalVector; << 109 delete[] fLastFinalVector; 96 delete [] fLastDyDx; << 110 delete[] fLastDyDx; 97 delete [] fMidVector; << 111 delete[] fMidVector; 98 delete [] fMidError; << 112 delete[] fMidError; 99 113 100 delete fAuxStepper; 114 delete fAuxStepper; 101 } 115 } 102 116 103 // The following coefficients have been obtain << 117 //The following coefficients have been obtained from 104 // Table 1: The Coefficients of the new pair 118 // Table 1: The Coefficients of the new pair 105 // << 119 //---Ref--- 106 // C. Tsitouras, "Runge–Kutta pairs of order << 120 // C. Tsitouras, “Runge–Kutta pairs of order 5(4) satisfying only 107 // the first column simplifying << 121 // the first column simplifying assumption,” 108 // Computers & Mathematics with Applications, << 122 // Computers & Mathematics with Applications, 109 // << 123 // vol. 62, no. 2, pp. 770–775, 2011. 110 // A corresponding matlab code was also found << 124 //----------------------------------- 111 // http://users.ntua.gr/tsitoura/new54.m << 125 // A corresponding matlab code was also found @ http://users.ntua.gr/tsitoura/new54.m 112 // << 126 113 // Doing a step 127 // Doing a step 114 // << 115 void 128 void 116 G4TsitourasRK45::Stepper( const G4double yInpu << 129 G4TsitourasRK45::Stepper( const G4double yInput[], 117 const G4double dydx[ << 130 const G4double dydx[], 118 G4double Step, << 131 G4double Step, 119 G4double yOut[ << 132 G4double yOut[], 120 G4double yErr[ << 133 G4double yErr[]) 121 { 134 { 122 const G4double b21 = 0.161 , << 135 G4int i; 123 b31 = -0.008480655492356988 << 136 124 b32 = 0.335480655492356989 << 137 const G4double 125 << 138 b21 = 0.161 , 126 b41 = 2.89715305710549343 << 139 127 b42 = -6.35944848997507484 << 140 b31 = -0.00848065549235698854 , 128 b43 = 4.36229543286958141 , << 141 b32 = 0.335480655492356989 , 129 << 130 b51 = 5.325864828439257, << 131 b52 = -11.748883564062828, << 132 b53 = 7.49553934288983621 , << 133 b54 = -0.09249506636175525, << 134 << 135 b61 = 5.8614554429464200, << 136 b62 = -12.9209693178471093 << 137 b63 = 8.1593678985761586 , << 138 b64 = -0.071584973281400997 << 139 b65 = -0.028269050394068382 << 140 << 141 b71 = 0.0964607668180652295 << 142 b72 = 0.01, << 143 b73 = 0.479889650414499575, << 144 b74 = 1.37900857410374189, << 145 b75 = -3.2900695154360807, << 146 b76 = 2.32471052409977398, << 147 << 148 // c1 = 0.001780011052226 , << 149 // c2 = 0.000816434459657 , << 150 // c3 = -0.007880878010262 , << 151 // c4 = 0.144711007173263 , << 152 // c5 = -0.582357165452555 , << 153 // c6 = 0.458082105929187 , << 154 // c7 = 1.0/66.0 ; << 155 << 156 dc1 = 0.0935237485818927066 << 157 dc2 = 0.0086528831415663676 << 158 dc3 = 0.492893099131431868 << 159 dc4 = 1.14023541226785810 - << 160 dc5 = - 2.3291801924393646 << 161 dc6 = 1.56887504931661552 - << 162 dc7 = 0.025; //- 1.0/66.0 ; << 163 << 164 // dc1 = -3.0/1280.0, << 165 // dc2 = 0.0, << 166 // dc3 = 6561.0/632320.0, << 167 // dc4 = -343.0/20800.0, << 168 // dc5 = 243.0/12800.0, << 169 // dc6 = -1.0/95.0, << 170 // dc7 = 0.0 ; << 171 142 172 const G4int numberOfVariables = GetNumberO << 143 b41 = 2.89715305710549343 , >> 144 b42 = -6.35944848997507484 , >> 145 b43 = 4.36229543286958141 , >> 146 >> 147 b51 = 5.325864828439257, >> 148 b52 = -11.748883564062828, >> 149 b53 = 7.49553934288983621 , >> 150 b54 = -0.09249506636175525, >> 151 >> 152 b61 = 5.8614554429464200, >> 153 b62 = -12.9209693178471093 , >> 154 b63 = 8.1593678985761586 , >> 155 b64 = -0.071584973281400997, >> 156 b65 = -0.0282690503940683829, >> 157 >> 158 >> 159 b71 = 0.0964607668180652295 , >> 160 b72 = 0.01, >> 161 b73 = 0.479889650414499575, >> 162 b74 = 1.37900857410374189, >> 163 b75 = -3.2900695154360807, >> 164 b76 = 2.32471052409977398, >> 165 >> 166 // c1 = 0.001780011052226 , >> 167 // c2 = 0.000816434459657 , >> 168 // c3 = -0.007880878010262 , >> 169 // c4 = 0.144711007173263 , >> 170 // c5 = -0.582357165452555 , >> 171 // c6 = 0.458082105929187 , >> 172 // c7 = 1.0/66.0 ; >> 173 >> 174 dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226, >> 175 dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657, >> 176 dc3 = 0.492893099131431868 - b73 , // + 0.007880878010262, >> 177 dc4 = 1.14023541226785810 - b74 , // 0.144711007173263, >> 178 dc5 = - 2.3291801924393646 - b75, // + 0.582357165452555, >> 179 dc6 = 1.56887504931661552 - b76 , // - 0.458082105929187, >> 180 dc7 = 0.025; //- 1.0/66.0 ; >> 181 >> 182 // dc1 = -3.0/1280.0, >> 183 // dc2 = 0.0, >> 184 // dc3 = 6561.0/632320.0, >> 185 // dc4 = -343.0/20800.0, >> 186 // dc5 = 243.0/12800.0, >> 187 // dc6 = -1.0/95.0, >> 188 // dc7 = 0.0 ; >> 189 >> 190 const G4int numberOfVariables= this->GetNumberOfVariables(); 173 191 174 // The number of variables to be integrate 192 // The number of variables to be integrated over 175 // << 193 yOut[7] = yTemp[7] = yIn[7]; 176 yOut[7] = yTemp[7] = yIn[7] = yInput[7]; << 177 << 178 // Saving yInput because yInput and yOut 194 // Saving yInput because yInput and yOut can be aliases for same array 179 // << 195 180 for(G4int i=0; i<numberOfVariables; ++i) << 196 for(i=0;i<numberOfVariables;i++) 181 { 197 { 182 yIn[i]=yInput[i]; 198 yIn[i]=yInput[i]; 183 } 199 } 184 // RightHandSide(yIn, dydx) ; // 1st Ste << 200 >> 201 // RightHandSide(yIn, dydx) ; >> 202 // 1st Step - Not doing, getting passed 185 203 186 for(G4int i=0; i<numberOfVariables; ++i) << 204 for(i=0;i<numberOfVariables;i++) 187 { 205 { 188 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; 206 yTemp[i] = yIn[i] + b21*Step*dydx[i] ; 189 } 207 } 190 RightHandSide(yTemp, ak2) ; / 208 RightHandSide(yTemp, ak2) ; // 2nd Stage 191 209 192 for(G4int i=0; i<numberOfVariables; ++i) << 210 for(i=0;i<numberOfVariables;i++) 193 { 211 { 194 yTemp[i] = yIn[i] + Step*(b31*dydx[i] 212 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ; 195 } 213 } 196 RightHandSide(yTemp, ak3) ; / 214 RightHandSide(yTemp, ak3) ; // 3rd Stage 197 215 198 for(G4int i=0; i<numberOfVariables; ++i) << 216 for(i=0;i<numberOfVariables;i++) 199 { 217 { 200 yTemp[i] = yIn[i] + Step*(b41*dydx[i] 218 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ; 201 } 219 } 202 RightHandSide(yTemp, ak4) ; / 220 RightHandSide(yTemp, ak4) ; // 4th Stage 203 221 204 for(G4int i=0; i<numberOfVariables; ++i) << 222 for(i=0;i<numberOfVariables;i++) 205 { 223 { 206 yTemp[i] = yIn[i] + Step*(b51*dydx[i] 224 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] + 207 b54*ak4[i]) 225 b54*ak4[i]) ; 208 } 226 } 209 RightHandSide(yTemp, ak5) ; / 227 RightHandSide(yTemp, ak5) ; // 5th Stage 210 228 211 for(G4int i=0; i<numberOfVariables; ++i) << 229 for(i=0;i<numberOfVariables;i++) 212 { 230 { 213 yTemp[i] = yIn[i] + Step*(b61*dydx[i] 231 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] + 214 b64*ak4[i] + 232 b64*ak4[i] + b65*ak5[i]) ; 215 } 233 } 216 RightHandSide(yTemp, ak6) ; / 234 RightHandSide(yTemp, ak6) ; // 6th Stage 217 235 218 for(G4int i=0; i<numberOfVariables; ++i) << 236 for(i=0;i<numberOfVariables;i++) 219 { 237 { 220 yOut[i] = yIn[i] + Step*(b71*dydx[i] + 238 yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] + 221 b74*ak4[i] + 239 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]); 222 } 240 } 223 RightHandSide(yOut, ak7); / << 241 RightHandSide(yOut, ak7); //7th Stage 224 242 225 //Calculate the error in the step: 243 //Calculate the error in the step: 226 for(G4int i=0; i<numberOfVariables; ++i) << 244 for(i=0;i<numberOfVariables;i++) 227 { 245 { 228 yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[ 246 yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] + 229 dc5*ak5[i] + dc6*ak6[ 247 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ; 230 248 231 // Store Input and Final values, for p 249 // Store Input and Final values, for possible use in calculating chord 232 // << 233 fLastInitialVector[i] = yIn[i] ; 250 fLastInitialVector[i] = yIn[i] ; 234 fLastFinalVector[i] = yOut[i]; 251 fLastFinalVector[i] = yOut[i]; 235 fLastDyDx[i] = dydx[i]; 252 fLastDyDx[i] = dydx[i]; 236 } 253 } 237 254 238 fLastStepLength = Step; 255 fLastStepLength = Step; 239 256 240 return ; 257 return ; 241 } 258 } 242 259 243 void G4TsitourasRK45::SetupInterpolation() << 260 void G4TsitourasRK45::SetupInterpolation() // (const G4double *yInput, const G4double *dydx, const G4double Step) 244 // (const G4double *yInput, const G4double * << 245 { 261 { 246 // Nothing to be done << 262 //Nothing to be done 247 } 263 } 248 264 249 void G4TsitourasRK45::Interpolate(const G4doub << 265 250 const G4doub << 266 void G4TsitourasRK45::Interpolate(const G4double *yInput, const G4double *dydx, const G4double Step, G4double *yOut, G4double tau){ 251 const G4doub << 267 252 G4doub << 268 253 G4doub << 254 { << 255 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7 269 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7; 256 // Coefficients for all the seven stages << 270 // Coefficients for all the seven stages. 257 271 258 const G4int numberOfVariables = GetNumberO << 272 const G4int numberOfVariables= this->GetNumberOfVariables(); 259 273 260 G4double tau0 = tau; 274 G4double tau0 = tau; 261 275 262 for(G4int i=0; i<numberOfVariables; ++i) << 276 for(int i=0;i<numberOfVariables;i++) 263 { 277 { 264 yIn[i] = yInput[i]; << 278 yIn[i]=yInput[i]; 265 } 279 } 266 280 267 G4double tau_2 = tau0*tau0 ; << 281 G4double 268 // tau_3 = tau0*tau_2, << 282 tau_2 = tau0*tau0 ; 269 // tau_4 = tau_2*tau_2; << 283 // tau_3 = tau0*tau_2, >> 284 // tau_4 = tau_2*tau_2; 270 285 271 bf1 = -1.0530884977290216*tau*(tau - 1.329 286 bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 - 272 1.4364028541716351*tau + 0.713 << 287 1.4364028541716351*tau + 0.7139816917074209), 273 bf2 = 0.1017*tau_2*(tau_2 - 2.196656833824 288 bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau + 274 1.2949852507374631); << 289 1.2949852507374631), 275 bf3 = 2.490627285651252793*tau_2*(tau_2 - 290 bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau 276 + 1.5780 << 291 + 1.57803468208092486) , 277 bf4 = -16.54810288924490272*(tau - 1.21712 292 bf4 = -16.54810288924490272*(tau - 1.21712927295533244)* 278 (tau - 0.6 << 293 (tau - 0.61620406037800089)*tau_2, 279 bf5 = 47.37952196281928122*(tau - 1.203071 294 bf5 = 47.37952196281928122*(tau - 1.203071208372362603)* 280 (tau - 0.6 << 295 (tau - 0.658047292653547382)*tau_2, 281 bf6 = -34.87065786149660974*(tau - 1.2)*(t 296 bf6 = -34.87065786149660974*(tau - 1.2)*(tau - 282 0.666666666666 << 297 0.666666666666666667)*tau_2, 283 bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2; 298 bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2; 284 299 285 // Putting together the coefficients calcu << 300 //Putting together the coefficients calculated as the respective stage coefficients 286 // stage coefficients << 301 for( int i=0; i<numberOfVariables; i++){ 287 // << 302 yOut[i] = yIn[i] + Step*( bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i] 288 for(G4int i=0; i<numberOfVariables; ++i) << 303 + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ; 289 { << 290 yOut[i] = yIn[i] + Step*( bf1*dydx[i] << 291 + bf4*ak4[i] << 292 + bf7*ak7[i] << 293 } 304 } 294 } 305 } 295 306 >> 307 /////////////////////////////////////////////////////////////////////////////// >> 308 296 G4double G4TsitourasRK45::DistChord() const 309 G4double G4TsitourasRK45::DistChord() const 297 { 310 { 298 G4double distLine, distChord; 311 G4double distLine, distChord; 299 G4ThreeVector initialPoint, finalPoint, midP 312 G4ThreeVector initialPoint, finalPoint, midPoint; 300 313 301 // Store last initial and final points (they << 314 // Store last initial and final points (they will be overwritten in self-Stepper call!) 302 // overwritten in self-Stepper call!) << 303 // << 304 initialPoint = G4ThreeVector( fLastInitialVe 315 initialPoint = G4ThreeVector( fLastInitialVector[0], 305 fLastInitialVe 316 fLastInitialVector[1], fLastInitialVector[2]); 306 finalPoint = G4ThreeVector( fLastFinalVect 317 finalPoint = G4ThreeVector( fLastFinalVector[0], 307 fLastFinalVect 318 fLastFinalVector[1], fLastFinalVector[2]); 308 319 309 // Do half a step using StepNoErr 320 // Do half a step using StepNoErr 310 // << 321 311 fAuxStepper->Stepper( fLastInitialVector, fL 322 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 312 fMidVector, fMidError ); 323 fMidVector, fMidError ); 313 324 314 midPoint = G4ThreeVector( fMidVector[0], fMi 325 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]); 315 326 316 // Use stored values of Initial and Endpoint 327 // Use stored values of Initial and Endpoint + new Midpoint to evaluate 317 // distance of Chord 328 // distance of Chord 318 // << 329 >> 330 319 if (initialPoint != finalPoint) 331 if (initialPoint != finalPoint) 320 { 332 { 321 distLine = G4LineSection::Distline( midP 333 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint ); 322 distChord = distLine; 334 distChord = distLine; 323 } 335 } 324 else 336 else 325 { 337 { 326 distChord = (midPoint-initialPoint).mag() 338 distChord = (midPoint-initialPoint).mag(); 327 } 339 } 328 return distChord; 340 return distChord; 329 } 341 } 330 342