From: Mathias Fr?h. <fro...@us...> - 2004-04-23 11:59:22
|
Update of /cvsroot/jsbsim/JSBSim In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv21252 Modified Files: Tag: FROHLICH FGDopri5.cpp FGDopri5.h FGExplicitEuler.cpp FGExplicitEuler.h FGImplicitEuler.cpp FGImplicitEuler.h FGMcFarland.cpp FGMcFarland.h FGPropagate.cpp FGPropagate.h FGTimestep.cpp FGTimestep.h Log Message: Move the current vehicle state and the current time value to FGPropagate. This fixes the known crash at exit/reset. Index: FGDopri5.cpp =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGDopri5.cpp,v retrieving revision 1.1.2.1 retrieving revision 1.1.2.2 diff -C2 -r1.1.2.1 -r1.1.2.2 *** FGDopri5.cpp 18 Apr 2004 09:52:48 -0000 1.1.2.1 --- FGDopri5.cpp 23 Apr 2004 11:59:10 -0000 1.1.2.2 *************** *** 100,136 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGDopri5::Integrate(Time tEnd) { // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd ) ) { // Cut the timestep if required. ! double curStepsize = ClipStepsizeToTEnd( tEnd ); // Compute the inner stages k1, ..., k7 ! FGVehicleState k1 = ComputeStateDeriv( mTime, mState ); ! FGVehicleState y2 = mState + (curStepsize*a21)*k1; ! FGVehicleState k2 = ComputeStateDeriv( mTime+c2*curStepsize, y2 ); ! FGVehicleState y3 = mState + (curStepsize*a31)*k1 + (curStepsize*a32)*k2; ! FGVehicleState k3 = ComputeStateDeriv( mTime+c3*curStepsize, y3 ); ! FGVehicleState y4 = mState + (curStepsize*a41)*k1 + (curStepsize*a42)*k2 + (curStepsize*a43)*k3; ! FGVehicleState k4 = ComputeStateDeriv( mTime+c4*curStepsize, y4 ); ! FGVehicleState y5 = mState + (curStepsize*a51)*k1 + (curStepsize*a52)*k2 + (curStepsize*a53)*k3 + (curStepsize*a54)*k4; ! FGVehicleState k5 = ComputeStateDeriv( mTime+c5*curStepsize, y5 ); ! FGVehicleState y6 = mState + (curStepsize*a61)*k1 + (curStepsize*a62)*k2 + (curStepsize*a63)*k3 + (curStepsize*a64)*k4 + (curStepsize*a65)*k5; ! FGVehicleState k6 = ComputeStateDeriv( mTime+curStepsize, y6 ); ! FGVehicleState y7 = mState + (curStepsize*a71)*k1 + (curStepsize*a73)*k3 + (curStepsize*a74)*k4 + (curStepsize*a75)*k5 + (curStepsize*a76)*k6; ! FGVehicleState k7 = ComputeStateDeriv( mTime+curStepsize, y7 ); // Dense output ... --- 100,136 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGDopri5::Integrate(Time tEnd, Time& t, FGVehicleState& state) { // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd, t ) ) { // Cut the timestep if required. ! double curStepsize = ClipStepsizeToTEnd( tEnd, t, GetStepsize() ); // Compute the inner stages k1, ..., k7 ! FGVehicleState k1 = ComputeStateDeriv( t, state ); ! FGVehicleState y2 = state + (curStepsize*a21)*k1; ! FGVehicleState k2 = ComputeStateDeriv( t+c2*curStepsize, y2 ); ! FGVehicleState y3 = state + (curStepsize*a31)*k1 + (curStepsize*a32)*k2; ! FGVehicleState k3 = ComputeStateDeriv( t+c3*curStepsize, y3 ); ! FGVehicleState y4 = state + (curStepsize*a41)*k1 + (curStepsize*a42)*k2 + (curStepsize*a43)*k3; ! FGVehicleState k4 = ComputeStateDeriv( t+c4*curStepsize, y4 ); ! FGVehicleState y5 = state + (curStepsize*a51)*k1 + (curStepsize*a52)*k2 + (curStepsize*a53)*k3 + (curStepsize*a54)*k4; ! FGVehicleState k5 = ComputeStateDeriv( t+c5*curStepsize, y5 ); ! FGVehicleState y6 = state + (curStepsize*a61)*k1 + (curStepsize*a62)*k2 + (curStepsize*a63)*k3 + (curStepsize*a64)*k4 + (curStepsize*a65)*k5; ! FGVehicleState k6 = ComputeStateDeriv( t+curStepsize, y6 ); ! FGVehicleState y7 = state + (curStepsize*a71)*k1 + (curStepsize*a73)*k3 + (curStepsize*a74)*k4 + (curStepsize*a75)*k5 + (curStepsize*a76)*k6; ! FGVehicleState k7 = ComputeStateDeriv( t+curStepsize, y7 ); // Dense output ... *************** *** 146,150 **** - double en = FGTimestep::ErrorNorm2( 1e-3, y7.GetLocation(), err.GetLocation() ); --- 146,149 ---- *************** *** 157,167 **** // We do unconditionally accept all steps, y7 is the new state value. ! mState = y7; // Increment the simulation time ... ! IncrementTime( curStepsize ); // Renormalize the quaternion and the location ! mState.Normalize(); } --- 156,166 ---- // We do unconditionally accept all steps, y7 is the new state value. ! state = y7; // Increment the simulation time ... ! t += curStepsize; // Renormalize the quaternion and the location ! state.Normalize(); } Index: FGDopri5.h =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGDopri5.h,v retrieving revision 1.1.2.1 retrieving revision 1.1.2.2 diff -C2 -r1.1.2.1 -r1.1.2.2 *** FGDopri5.h 18 Apr 2004 09:52:48 -0000 1.1.2.1 --- FGDopri5.h 23 Apr 2004 11:59:10 -0000 1.1.2.2 *************** *** 93,97 **** Perform the actual timesteps nothing exciting here ... */ ! virtual bool Integrate(Time tEnd); private: --- 93,97 ---- Perform the actual timesteps nothing exciting here ... */ ! virtual bool Integrate(Time tEnd, Time& t, FGVehicleState& state); private: Index: FGExplicitEuler.cpp =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGExplicitEuler.cpp,v retrieving revision 1.1.2.2 retrieving revision 1.1.2.3 diff -C2 -r1.1.2.2 -r1.1.2.3 *** FGExplicitEuler.cpp 18 Apr 2004 08:41:52 -0000 1.1.2.2 --- FGExplicitEuler.cpp 23 Apr 2004 11:59:10 -0000 1.1.2.3 *************** *** 59,78 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGExplicitEuler::Integrate(Time tEnd) { // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd ) ) { // Cut the timestep if required. ! double curStepsize = ClipStepsizeToTEnd( tEnd ); // Do one explicit euler step. ! mState += curStepsize * ComputeStateDeriv( mTime, mState ); // Increment the simulation time ... ! IncrementTime( curStepsize ); // Renormalize the quaternion and the location ! mState.Normalize(); } --- 59,78 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGExplicitEuler::Integrate(Time tEnd, Time& t, FGVehicleState& state) { // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd, t ) ) { // Cut the timestep if required. ! double curStepsize = ClipStepsizeToTEnd( tEnd, t, GetStepsize() ); // Do one explicit euler step. ! state += curStepsize * ComputeStateDeriv( t, state ); // Increment the simulation time ... ! t += curStepsize; // Renormalize the quaternion and the location ! state.Normalize(); } Index: FGExplicitEuler.h =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGExplicitEuler.h,v retrieving revision 1.1.2.1 retrieving revision 1.1.2.2 diff -C2 -r1.1.2.1 -r1.1.2.2 *** FGExplicitEuler.h 8 Apr 2004 19:37:51 -0000 1.1.2.1 --- FGExplicitEuler.h 23 Apr 2004 11:59:10 -0000 1.1.2.2 *************** *** 95,99 **** Perform the actual timesteps nothing exciting here ... */ ! virtual bool Integrate(Time tEnd); }; --- 95,99 ---- Perform the actual timesteps nothing exciting here ... */ ! virtual bool Integrate(Time tEnd, Time& t, FGVehicleState& state); }; Index: FGImplicitEuler.cpp =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGImplicitEuler.cpp,v retrieving revision 1.1.2.2 retrieving revision 1.1.2.3 diff -C2 -r1.1.2.2 -r1.1.2.3 *** FGImplicitEuler.cpp 20 Apr 2004 10:21:29 -0000 1.1.2.2 --- FGImplicitEuler.cpp 23 Apr 2004 11:59:10 -0000 1.1.2.3 *************** *** 67,71 **** { // FIXME: ! return mState + curStepsize*ComputeStateDeriv( mTime+curStepsize, x ) - x; } --- 67,71 ---- { // FIXME: ! return mCurrentState + curStepsize*ComputeStateDeriv( mCurrentTime+curStepsize, x ) - x; } *************** *** 75,79 **** { // compute a jacobian ... ! EOMJac = ComputeJacobian( mTime+curStepsize, x ); ImplEulJacStepsize = 0.0; // cout << EOMJac << endl; --- 75,79 ---- { // compute a jacobian ... ! EOMJac = ComputeJacobian( mCurrentTime+curStepsize, x ); ImplEulJacStepsize = 0.0; // cout << EOMJac << endl; *************** *** 116,120 **** { double en ! = FGTimestep::ErrorNorm2( 1e-3, ref.GetLocation(), diff.GetLocation() ); en += FGTimestep::ErrorNorm2( 1e-4, diff.GetOrientation() ); en += FGTimestep::ErrorNorm2( 1e-3, 1e-4, ref.GetUVW(), diff.GetUVW() ); --- 116,120 ---- { double en ! = FGTimestep::ErrorNorm2( 1e-4, ref.GetLocation(), diff.GetLocation() ); en += FGTimestep::ErrorNorm2( 1e-4, diff.GetOrientation() ); en += FGTimestep::ErrorNorm2( 1e-3, 1e-4, ref.GetUVW(), diff.GetUVW() ); *************** *** 126,134 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGImplicitEuler::Integrate(Time tEnd) { // Make shure we have a valid jacobian. ! if (ImplEulJac.Rows() != mState.Rows()) { ! ComputeNewJacobian( mState ); // set default stepsize curStepsize = GetStepsize(); --- 126,134 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGImplicitEuler::Integrate(Time tEnd, Time& t, FGVehicleState& state) { // Make shure we have a valid jacobian. ! if (ImplEulJac.Rows() != state.Rows()) { ! ComputeNewJacobian( state ); // set default stepsize curStepsize = GetStepsize(); *************** *** 138,147 **** if (mSteps%60 == 0) { cout << "NewStep: " ! << mState.GetLatitude()*radtodeg << ", " ! << mState.GetLongitude()*radtodeg << ", " ! << mFDMEx->GetAuxiliary()->Geth( mState ) << ", " ! << mState.GetUVW() << ", " ! << mState.GetEuler() << ", " ! << mState.GetPQR() << ", " << endl; } --- 138,147 ---- if (mSteps%60 == 0) { cout << "NewStep: " ! << state.GetLatitude()*radtodeg << ", " ! << state.GetLongitude()*radtodeg << ", " ! << mFDMEx->GetAuxiliary()->Geth( state ) << ", " ! << state.GetUVW() << ", " ! << state.GetEuler() << ", " ! << state.GetPQR() << ", " << endl; } *************** *** 150,157 **** // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd ) ) { // Cut the timestep if required. ! curStepsize = ClipStepsizeToTEnd( tEnd, 2.0*curStepsize ); int newtonmaxit = 6; --- 150,157 ---- // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd, t ) ) { // Cut the timestep if required. ! curStepsize = ClipStepsizeToTEnd( tEnd, t, 2.0*curStepsize ); int newtonmaxit = 6; *************** *** 159,174 **** int k = 0; // An explicit euler step for the initial guess ... // If someone has something better please cry!!!!! ! // FGVehicleState new_state = mState + (0.5*curStepsize)*ComputeStateDeriv( mTime, mState ); // .. ok, first try ... ! FGVehicleState new_state = mState; // ... not so good ... ! // ((FGVehicleVelocity&)new_state) += (0.2*curStepsize)*ComputeVelocityDeriv( mTime, mState ); // MAybe ... // Since it is stiff in the velocity space, predicting the position // with an explicit step might be a good idea ... ! // ((FGVehiclePosition&)new_state) += (0.5*curStepsize)*ComputePositionDeriv( mTime, mState ); ! // ((FGVehiclePosition&)new_state) += curStepsize*ComputePositionDeriv( mTime, mState ); bool converged; do { --- 159,178 ---- int k = 0; + // copy this to private storage, is accessed during the newton iteration. + mCurrentState = state; + mCurrentTime = t; + // An explicit euler step for the initial guess ... // If someone has something better please cry!!!!! ! // FGVehicleState new_state = state + (0.5*curStepsize)*ComputeStateDeriv( t, state ); // .. ok, first try ... ! FGVehicleState new_state = state; // ... not so good ... ! // ((FGVehicleVelocity&)new_state) += (0.2*curStepsize)*ComputeVelocityDeriv( t, state ); // MAybe ... // Since it is stiff in the velocity space, predicting the position // with an explicit step might be a good idea ... ! // ((FGVehiclePosition&)new_state) += (0.5*curStepsize)*ComputePositionDeriv( t, state ); ! // ((FGVehiclePosition&)new_state) += curStepsize*ComputePositionDeriv( t, state ); bool converged; do { *************** *** 186,190 **** curStepsize *= 0.5; new_state *= 0.2; ! new_state += 0.8*mState; ComputeNewJacobian( new_state ); newtonmaxit += 2; --- 190,194 ---- curStepsize *= 0.5; new_state *= 0.2; ! new_state += 0.8*state; ComputeNewJacobian( new_state ); newtonmaxit += 2; *************** *** 195,204 **** // cout << "Success " << endl; // successful step ... ! mState = new_state; // If it took tooo long to iterate in the newton method, // try with a new jacobian in the next step ... if ( maxnewtontry < k + 3 ) ! ComputeNewJacobian( mState ); } else { --- 199,208 ---- // cout << "Success " << endl; // successful step ... ! state = new_state; // If it took tooo long to iterate in the newton method, // try with a new jacobian in the next step ... if ( maxnewtontry < k + 3 ) ! ComputeNewJacobian( state ); } else { *************** *** 209,220 **** // curStepsize = ClipStepsizeToTEnd( tEnd ); ! mState += curStepsize*ComputeStateDeriv( mTime, mState ); } // Increment the simulation time ... ! IncrementTime( curStepsize ); // Renormalize the quaternion and the location ! mState.Normalize(); ++mSteps; --- 213,224 ---- // curStepsize = ClipStepsizeToTEnd( tEnd ); ! state += curStepsize*ComputeStateDeriv( t, state ); } // Increment the simulation time ... ! t += curStepsize; // Renormalize the quaternion and the location ! state.Normalize(); ++mSteps; Index: FGImplicitEuler.h =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGImplicitEuler.h,v retrieving revision 1.1.2.2 retrieving revision 1.1.2.3 diff -C2 -r1.1.2.2 -r1.1.2.3 *** FGImplicitEuler.h 20 Apr 2004 10:21:29 -0000 1.1.2.2 --- FGImplicitEuler.h 23 Apr 2004 11:59:10 -0000 1.1.2.3 *************** *** 91,95 **** Perform the actual timesteps nothing exciting here ... */ ! virtual bool Integrate(Time tEnd); private: --- 91,95 ---- Perform the actual timesteps nothing exciting here ... */ ! virtual bool Integrate(Time tEnd, Time& t, FGVehicleState& state); private: *************** *** 105,108 **** --- 105,110 ---- double ImplEulJacStepsize; + FGVehicleState mCurrentState; + double mCurrentTime; double curStepsize; Index: FGMcFarland.cpp =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGMcFarland.cpp,v retrieving revision 1.1.2.2 retrieving revision 1.1.2.3 diff -C2 -r1.1.2.2 -r1.1.2.3 *** FGMcFarland.cpp 18 Apr 2004 08:41:52 -0000 1.1.2.2 --- FGMcFarland.cpp 23 Apr 2004 11:59:10 -0000 1.1.2.3 *************** *** 62,73 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGMcFarland::Integrate(Time tEnd) { // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd ) ) { // Cut the timestep if required. ! double curStepsize = ClipStepsizeToTEnd( tEnd ); // Check if we have enough history for a regular step. --- 62,73 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGMcFarland::Integrate(Time tEnd, Time& t, FGVehicleState& state) { // Do the actual timesteps. // Cycle as long as required to reach tEnd. ! while ( !ReachedTEnd( tEnd, t ) ) { // Cut the timestep if required. ! double curStepsize = ClipStepsizeToTEnd( tEnd, t, GetStepsize() ); // Check if we have enough history for a regular step. *************** *** 78,90 **** // approach requires constant stepsizes this is perfectly ok. if ( mPastStepSize == curStepsize ) ! McFarlandStep( curStepsize ); else ! EulerStep( curStepsize ); // Increment the simulation time ... ! IncrementTime( curStepsize ); // Renormalize the quaternion and the location ! mState.Normalize(); } --- 78,90 ---- // approach requires constant stepsizes this is perfectly ok. if ( mPastStepSize == curStepsize ) ! McFarlandStep( curStepsize, t, state ); else ! EulerStep( curStepsize, t, state ); // Increment the simulation time ... ! t += curStepsize; // Renormalize the quaternion and the location ! state.Normalize(); } *************** *** 102,134 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! void FGMcFarland::EulerStep(Time stepSize) { // Do one explicit euler step at the beginning ... ! FGVehicleState stateDeriv = ComputeStateDeriv( mTime, mState ); ! mState += stepSize * stateDeriv; // ... and save the history. mPastStepSize = stepSize; mPastVelDeriv = stateDeriv; ! mPastPosDeriv = ComputePositionDeriv( mTime+stepSize, mState ); } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! void FGMcFarland::McFarlandStep(Time stepSize) { // Get the velocity derivatives at the current state. ! FGVehicleVelocity curVelDeriv = ComputeVelocityDeriv( mTime, mState ); // Get a reference to the subvectors. ! FGVehiclePosition& curPositions = mState; ! FGVehicleVelocity& curVelocities = mState; // Adams Barnsworth is // y_{n+1} = y_n + 0.5*h*(3*y'_n - y'_{n-1}) ! curVelocities += 1.5*stepSize*curVelDeriv; ! curVelocities -= 0.5*stepSize*mPastVelDeriv; ! // Save the old derivatives. mPastVelDeriv = curVelDeriv; --- 102,134 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! void FGMcFarland::EulerStep(Time stepSize, Time& t, FGVehicleState& state) { // Do one explicit euler step at the beginning ... ! FGVehicleState stateDeriv = ComputeStateDeriv( t, state ); ! state += stepSize * stateDeriv; // ... and save the history. mPastStepSize = stepSize; mPastVelDeriv = stateDeriv; ! mPastPosDeriv = ComputePositionDeriv( t+stepSize, state ); } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! void FGMcFarland::McFarlandStep(Time stepSize, Time& t, FGVehicleState& state) { // Get the velocity derivatives at the current state. ! FGVehicleVelocity curVelDeriv = ComputeVelocityDeriv( t, state ); // Get a reference to the subvectors. ! FGVehiclePosition& curPositions = state; ! FGVehicleVelocity& curVelocities = state; // Adams Barnsworth is // y_{n+1} = y_n + 0.5*h*(3*y'_n - y'_{n-1}) ! curVelocities += (1.5*stepSize)*curVelDeriv; ! curVelocities -= (0.5*stepSize)*mPastVelDeriv; ! // Save the current velocity derivatives as the old ones for the next step. mPastVelDeriv = curVelDeriv; *************** *** 141,150 **** // The trapezidual rule. // y_{n+1} = y_n + 0.5*h*(y'_{n+1} + y'_n) ! curPositions += 0.5*stepSize*mPastPosDeriv; mPastPosDeriv = ! ComputePositionDeriv( mTime+stepSize, ! FGVehicleState( mState, curVelocities ) ); ! curPositions += 0.5*stepSize*mPastPosDeriv; } --- 141,149 ---- // The trapezidual rule. // y_{n+1} = y_n + 0.5*h*(y'_{n+1} + y'_n) ! curPositions += (0.5*stepSize)*mPastPosDeriv; mPastPosDeriv = ! ComputePositionDeriv( t+stepSize, FGVehicleState( state, curVelocities ) ); ! curPositions += (0.5*stepSize)*mPastPosDeriv; } Index: FGMcFarland.h =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGMcFarland.h,v retrieving revision 1.1.2.1 retrieving revision 1.1.2.2 diff -C2 -r1.1.2.1 -r1.1.2.2 *** FGMcFarland.h 8 Apr 2004 19:37:59 -0000 1.1.2.1 --- FGMcFarland.h 23 Apr 2004 11:59:10 -0000 1.1.2.2 *************** *** 39,43 **** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ #include "FGTimestep.h" - #include "FGExplicitEuler.h" #include "FGVehicleState.h" --- 39,42 ---- *************** *** 92,96 **** /** Do timestepping. */ ! virtual bool Integrate(Time tEnd); /** Reset integration. --- 91,95 ---- /** Do timestepping. */ ! virtual bool Integrate(Time tEnd, Time& t, FGVehicleState& state); /** Reset integration. *************** *** 107,111 **** Also save the history required for the next step. */ ! void EulerStep(Time stepSize); /** McFarland step. --- 106,110 ---- Also save the history required for the next step. */ ! void EulerStep(Time stepSize, Time& t, FGVehicleState& state); /** McFarland step. *************** *** 116,120 **** Also save the history required for the next step. */ ! void McFarlandStep(Time stepSize); /** Stepsize of the last timestep we have history for. --- 115,119 ---- Also save the history required for the next step. */ ! void McFarlandStep(Time stepSize, Time& t, FGVehicleState& state); /** Stepsize of the last timestep we have history for. Index: FGPropagate.cpp =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/FGPropagate.cpp,v retrieving revision 1.1.2.5 retrieving revision 1.1.2.6 diff -C2 -r1.1.2.5 -r1.1.2.6 *** FGPropagate.cpp 20 Apr 2004 10:27:05 -0000 1.1.2.5 --- FGPropagate.cpp 23 Apr 2004 11:59:10 -0000 1.1.2.6 *************** *** 99,108 **** void FGPropagate::SetInitialState(const FGInitialCondition *FGIC) { - FGVehicleState state; // Set the position lat/lon/radius FGLocation newPos( FGIC->GetLongitudeRadIC(), FGIC->GetLatitudeRadIC(), FGIC->GetAltitudeFtIC() + FGIC->GetSeaLevelRadiusFtIC() ); ! state.SetLocation( newPos ); // Set the Orientation from the euler angles --- 99,107 ---- void FGPropagate::SetInitialState(const FGInitialCondition *FGIC) { // Set the position lat/lon/radius FGLocation newPos( FGIC->GetLongitudeRadIC(), FGIC->GetLatitudeRadIC(), FGIC->GetAltitudeFtIC() + FGIC->GetSeaLevelRadiusFtIC() ); ! mState.SetLocation( newPos ); // Set the Orientation from the euler angles *************** *** 110,114 **** FGIC->GetThetaRadIC(), FGIC->GetPsiRadIC() ); ! state.SetOrientation( newOrientation ); // Set the velocities in the instantaneus body frame --- 109,113 ---- FGIC->GetThetaRadIC(), FGIC->GetPsiRadIC() ); ! mState.SetOrientation( newOrientation ); // Set the velocities in the instantaneus body frame *************** *** 118,122 **** ic->GetVBodyFpsIC(), ic->GetWBodyFpsIC() ); ! state.SetUVW( newUVW ); // Set the angular velocities in the instantaneus body frame. --- 117,121 ---- ic->GetVBodyFpsIC(), ic->GetWBodyFpsIC() ); ! mState.SetUVW( newUVW ); // Set the angular velocities in the instantaneus body frame. *************** *** 124,131 **** FGIC->GetQRadpsIC(), FGIC->GetRRadpsIC() ); ! state.SetPQR( newPQR ); ! // Finally set it. ! mTimestep->SetState( state ); } --- 123,131 ---- FGIC->GetQRadpsIC(), FGIC->GetRRadpsIC() ); ! mState.SetPQR( newPQR ); ! // Tell the timestepping method, that it should throw away all history ! // which is invalid now. ! TimestepReset(); } *************** *** 138,141 **** --- 138,148 ---- return true; + // Some kind of assertation ... + if (!mTimestep) { + PutMessage( "Fatal error in FGPropagate: " + "reached FGPropagate::Run() without a timestepping method!" ); + return true; + } + // Successful return if propagateping is suspended. if (mSuspended) { *************** *** 143,148 **** // So compute the accelerations one time and exit ... // Will clean this up a later time. FIXME ! FGVehicleVelocity accel ! = mTimestep->ComputeVelocityDeriv( mTimestep->GetTime(), mTimestep->GetState() ); // HACK HACK HACK Auxiliary->GetUVWdot() = accel.GetUVW(); --- 150,154 ---- // So compute the accelerations one time and exit ... // Will clean this up a later time. FIXME ! FGVehicleVelocity accel = mTimestep->ComputeVelocityDeriv( mTime, mState ); // HACK HACK HACK Auxiliary->GetUVWdot() = accel.GetUVW(); *************** *** 153,163 **** // Do the actual timestepping. ! mTimestep->Integrate( mTimestep->GetTime() + rate*State->Getdt() ); // HACK HACK HACK // Finally update the acceleration to the realy current value // at interval end. ! FGVehicleVelocity accel ! = mTimestep->ComputeVelocityDeriv( mTimestep->GetTime(), mTimestep->GetState() ); Auxiliary->GetUVWdot() = accel.GetUVW(); Auxiliary->GetPQRdot() = accel.GetPQR(); --- 159,174 ---- // Do the actual timestepping. ! // Note that mTime and mState is modified during this call ! double tEnd = mTime + rate*State->Getdt(); ! if (!mTimestep->Integrate( tEnd, mTime, mState )) { ! PutMessage( "Fatal error in FGPropagate: " ! "Timestepping method did not return successful!" ); ! return true; ! } // HACK HACK HACK // Finally update the acceleration to the realy current value // at interval end. ! FGVehicleVelocity accel = mTimestep->ComputeVelocityDeriv( mTime, mState ); Auxiliary->GetUVWdot() = accel.GetUVW(); Auxiliary->GetPQRdot() = accel.GetPQR(); *************** *** 173,178 **** void FGPropagate::bind(void) { ! // FIXME really? ! ((FGVehicleState&)mTimestep->GetState()).bind( PropertyManager, "" ); PropertyManager->Tie("sim-time-sec", this, &FGPropagate::GetTime); } --- 184,190 ---- void FGPropagate::bind(void) { ! // Tie the current vehicle state into the property tree ! mState.bind( PropertyManager, "" ); ! // Tie the current time into the property tree PropertyManager->Tie("sim-time-sec", this, &FGPropagate::GetTime); } *************** *** 182,187 **** void FGPropagate::unbind(void) { ! // FIXME really? ! ((FGVehicleState&)mTimestep->GetState()).unbind( PropertyManager, "" ); PropertyManager->Untie("sim-time-sec"); } --- 194,199 ---- void FGPropagate::unbind(void) { ! // cleanup ... ! mState.unbind( PropertyManager, "" ); PropertyManager->Untie("sim-time-sec"); } Index: FGPropagate.h =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/FGPropagate.h,v retrieving revision 1.1.2.2 retrieving revision 1.1.2.3 diff -C2 -r1.1.2.2 -r1.1.2.3 *** FGPropagate.h 14 Apr 2004 20:11:26 -0000 1.1.2.2 --- FGPropagate.h 23 Apr 2004 11:59:10 -0000 1.1.2.3 *************** *** 106,114 **** /** Get the current simulation time. */ ! double GetTime(void) const { return mTimestep->GetTime(); } /** Set the current simulation time. */ ! void SetTime(double t) { mTimestep->SetTime(t); } /** Initializes the simulation state based on parameters from an --- 106,114 ---- /** Get the current simulation time. */ ! double GetTime(void) const { return mTime; } /** Set the current simulation time. */ ! void SetTime(double t) { mTime = t; TimestepReset(); } /** Initializes the simulation state based on parameters from an *************** *** 119,126 **** void SetInitialState(const FGInitialCondition *FGIC); ! const FGVehicleState& GetState(void) const { return mTimestep->GetState(); } private: void Debug(int from); --- 119,137 ---- void SetInitialState(const FGInitialCondition *FGIC); ! /** Current vehicle state. ! @return a const reference to the current vehicle state. ! Note that we never return a nonconst reference to the vehicle state ! since physics tells us that we can only change this state by applying ! a force. One exception: SetInitialState sets the complete vehicle state ! to the given initial argument. ! */ ! const FGVehicleState& GetState(void) const { return mState; } private: + /** Reset timestepping algorithm. + */ + void TimestepReset(void) { if (mTimestep) mTimestep->Reset(); } + void Debug(int from); *************** *** 128,134 **** --- 139,154 ---- void unbind(void); + // The current physical state of the aircraft. + FGVehicleState mState; + + // The current time. + double mTime; + // Is true if timestepping is suspended. bool mSuspended; + // The pointer to the currently used timestepping method. + // Can be zero if it is not yet configured. + // If we enter Run() it *must* be nonzerro ... FGTimestep* mTimestep; }; Index: FGTimestep.cpp =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGTimestep.cpp,v retrieving revision 1.1.2.6 retrieving revision 1.1.2.7 diff -C2 -r1.1.2.6 -r1.1.2.7 *** FGTimestep.cpp 20 Apr 2004 10:19:45 -0000 1.1.2.6 --- FGTimestep.cpp 23 Apr 2004 11:59:10 -0000 1.1.2.7 *************** *** 66,72 **** mFDMEx = fdmex; - // Reset Time - mTime = 0.0; - // Just for compatibility for now. mStepsize = 1.0/120.0; --- 66,69 ---- *************** *** 436,449 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGTimestep::ReachedTEnd( Time tEnd ) { ! // Check if mTime is equal to tEnd up to the order of roundoff. static double eps = 8.0*numeric_limits<double>::epsilon(); ! if ( eps*max(fabs(mTime), fabs(tEnd)) < fabs(mTime - tEnd) ) return false; ! // If it is equal up to roundoff, just set mTime to tEnd, so that others // do not get confused by this very small error. ! mTime = tEnd; return true; } --- 433,446 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! bool FGTimestep::ReachedTEnd( Time tEnd, Time& t ) { ! // Check if t is equal to tEnd up to the order of roundoff. static double eps = 8.0*numeric_limits<double>::epsilon(); ! if ( eps*max(fabs(t), fabs(tEnd)) < fabs(t - tEnd) ) return false; ! // If it is equal up to roundoff, just set t to tEnd, so that others // do not get confused by this very small error. ! t = tEnd; return true; } *************** *** 451,464 **** //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! double FGTimestep::ClipStepsizeToTEnd( Time tEnd, double stepsize ) const ! { ! return min( stepsize, tEnd - mTime ); ! } ! ! //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! double FGTimestep::ClipStepsizeToTEnd( Time tEnd ) const { ! return min( mStepsize, tEnd - mTime ); } --- 448,460 ---- //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! double FGTimestep::ClipStepsizeToTEnd( Time tEnd, Time t, double stepsize ) { ! // Clip the stepsize not to overshoot, ! // but don't bother with finite precision. ! static double eps = 4.0*numeric_limits<double>::epsilon(); ! if ( fabs(t + stepsize - tEnd) < eps*max(fabs(t+stepsize), fabs(tEnd)) ) ! return stepsize; ! else ! return min( stepsize, tEnd - t ); } Index: FGTimestep.h =================================================================== RCS file: /cvsroot/jsbsim/JSBSim/Attic/FGTimestep.h,v retrieving revision 1.1.2.7 retrieving revision 1.1.2.8 diff -C2 -r1.1.2.7 -r1.1.2.8 *** FGTimestep.h 18 Apr 2004 09:39:06 -0000 1.1.2.7 --- FGTimestep.h 23 Apr 2004 11:59:10 -0000 1.1.2.8 *************** *** 101,112 **** typedef double Time; - /** Get the current simulation time. - */ - Time GetTime(void) const { return mTime; } - - /** Set the current simulation time. - */ - void SetTime(Time t) { Reset(); mTime = t; } - /** Get the stepzize of the last timestep. */ --- 101,104 ---- *************** *** 117,132 **** void SetStepsize(double stepsize) { Reset(); mStepsize = stepsize; } - /** Set the statevalues to the given argument. - */ - void SetState(const FGVehicleState& s) { Reset(); mState = s; } - - /** Get the statevalues. - */ - const FGVehicleState& GetState(void) const { return mState; } - /** Do timestepping up to time tEnd. @param tEnd time to stop timestepping at. */ ! virtual bool Integrate(Time tEnd) = 0; /** Reset the integration scheme. */ --- 109,116 ---- void SetStepsize(double stepsize) { Reset(); mStepsize = stepsize; } /** Do timestepping up to time tEnd. @param tEnd time to stop timestepping at. */ ! virtual bool Integrate(Time tEnd, Time& t, FGVehicleState& state) = 0; /** Reset the integration scheme. */ *************** *** 150,176 **** protected: - // The current physical state of the aircraft. - FGVehicleState mState; - // Holds the actual stepsize of the timestepping method. // Since this is common to all timestepping methods this value is best // held here ... double mStepsize; - // ... the same goes for the actual time of integration. - Time mTime; /** Check if we reached tEnd. */ ! bool ReachedTEnd( Time tEnd ); /** Return the current stepsize clipped to tEnd. */ ! double ClipStepsizeToTEnd( Time tEnd, double stepsize ) const; ! double ClipStepsizeToTEnd( Time tEnd ) const; ! ! /** Increment the current time due to the given stepsize. ! */ ! void IncrementTime( double curStepsize ) { mTime += curStepsize; } ! /** The 'scaled norm' of x1-x2. --- 134,149 ---- protected: // Holds the actual stepsize of the timestepping method. // Since this is common to all timestepping methods this value is best // held here ... double mStepsize; /** Check if we reached tEnd. */ ! static bool ReachedTEnd( Time tEnd, Time& t ); /** Return the current stepsize clipped to tEnd. */ ! static double ClipStepsizeToTEnd( Time tEnd, Time t, double stepsize ); /** The 'scaled norm' of x1-x2. |