00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id: G4MagErrorStepper.cc 69786 2013-05-15 09:38:51Z gcosmo $ 00028 // 00029 // -------------------------------------------------------------------- 00030 00031 #include "G4MagErrorStepper.hh" 00032 #include "G4LineSection.hh" 00033 00034 G4MagErrorStepper::~G4MagErrorStepper() 00035 { 00036 delete[] yMiddle; 00037 delete[] dydxMid; 00038 delete[] yInitial; 00039 delete[] yOneStep; 00040 } 00041 00042 void 00043 G4MagErrorStepper::Stepper( const G4double yInput[], 00044 const G4double dydx[], 00045 G4double hstep, 00046 G4double yOutput[], 00047 G4double yError [] ) 00048 { 00049 const G4int nvar = this->GetNumberOfVariables() ; 00050 const G4int maxvar= GetNumberOfStateVariables(); 00051 00052 G4int i; 00053 // correction for Richardson Extrapolation. 00054 G4double correction = 1. / ( (1 << IntegratorOrder()) -1 ); 00055 00056 // Saving yInput because yInput and yOutput can be aliases for same array 00057 00058 for(i=0;i<nvar;i++) yInitial[i]=yInput[i]; 00059 yInitial[7]= yInput[7]; // Copy the time in case ... even if not really needed 00060 yMiddle[7] = yInput[7]; // Copy the time from initial value 00061 yOneStep[7] = yInput[7]; // As it contributes to final value of yOutput ? 00062 // yOutput[7] = yInput[7]; // -> dumb stepper does it too for RK4 00063 for(i=nvar;i<maxvar;i++) yOutput[i]=yInput[i]; 00064 // yError[7] = 0.0; 00065 00066 G4double halfStep = hstep * 0.5; 00067 00068 // Do two half steps 00069 00070 DumbStepper (yInitial, dydx, halfStep, yMiddle); 00071 RightHandSide(yMiddle, dydxMid); 00072 DumbStepper (yMiddle, dydxMid, halfStep, yOutput); 00073 00074 // Store midpoint, chord calculation 00075 00076 fMidPoint = G4ThreeVector( yMiddle[0], yMiddle[1], yMiddle[2]); 00077 00078 // Do a full Step 00079 DumbStepper(yInitial, dydx, hstep, yOneStep); 00080 for(i=0;i<nvar;i++) { 00081 yError [i] = yOutput[i] - yOneStep[i] ; 00082 yOutput[i] += yError[i]*correction ; // Provides accuracy increased 00083 // by 1 order via the 00084 // Richardson Extrapolation 00085 } 00086 00087 fInitialPoint = G4ThreeVector( yInitial[0], yInitial[1], yInitial[2]); 00088 fFinalPoint = G4ThreeVector( yOutput[0], yOutput[1], yOutput[2]); 00089 00090 return ; 00091 } 00092 00093 00094 00095 G4double 00096 G4MagErrorStepper::DistChord() const 00097 { 00098 // Estimate the maximum distance from the curve to the chord 00099 // 00100 // We estimate this using the distance of the midpoint to 00101 // chord (the line between 00102 // 00103 // Method below is good only for angle deviations < 2 pi, 00104 // This restriction should not a problem for the Runge cutta methods, 00105 // which generally cannot integrate accurately for large angle deviations. 00106 G4double distLine, distChord; 00107 00108 if (fInitialPoint != fFinalPoint) { 00109 distLine= G4LineSection::Distline( fMidPoint, fInitialPoint, fFinalPoint ); 00110 // This is a class method that gives distance of Mid 00111 // from the Chord between the Initial and Final points. 00112 00113 distChord = distLine; 00114 }else{ 00115 distChord = (fMidPoint-fInitialPoint).mag(); 00116 } 00117 00118 return distChord; 00119 } 00120