Copy the KDE 3.5 branch to branches/trinity for new KDE 3.5 features.

BUG:215923


git-svn-id: svn://anonsvn.kde.org/home/kde/branches/trinity/kdeartwork@1054174 283d02a7-25f6-0310-bc7c-ecb5cbfe19da

1 files changed, 252 insertions, 0 deletions

diff --git a/kscreensaver/kdesavers/rkodesolver.cpp b/kscreensaver/kdesavers/rkodesolver.cpp
new file mode 100644
index 00000000..2127a110
--- /dev/null
+++ b/kscreensaver/kdesavers/rkodesolver.cpp
@@ -0,0 +1,252 @@
+//============================================================================
+//
+// Ordinary differential equation solver using the Runge-Kutta method.
+// $Id$
+// Copyright (C) 2004 Georg Drenkhahn
+//
+// This file is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License version 2 as published by the
+// Free Software Foundation.
+//
+//============================================================================
+
+#include <kdebug.h>
+#include "rkodesolver.h"
+
+template<typename T>
+RkOdeSolver<T>::RkOdeSolver(const T &x,
+                            const std::valarray<T>& y,
+                            const T &dx,
+                            const T &eps)
+   : mX(x)
+{
+   Y(y);
+   dX(dx);
+   Eps(eps);
+}
+
+// virtual dtor
+template<typename T>
+RkOdeSolver<T>::~RkOdeSolver(void)
+{}
+
+// accessors
+
+template<typename T>
+void
+RkOdeSolver<T>::dX(const T &a)
+{
+   if (a < 0.0)
+   {
+      kdDebug() << "RkOdeSolver: dx was negative, made it positive" << endl;
+      mStep = -a;
+   }
+   else if (a == 0.0)
+   {
+      mStep = 0.001;            // a very arbitrary value
+      kdDebug() << "RkOdeSolver: dx == 0, set it to " << mStep << endl;
+   }
+   else
+   {
+      mStep = a;
+   }
+}
+
+template<typename T>
+void
+RkOdeSolver<T>::Eps(const T &a)
+{
+   if (a < 0.0)
+   {
+      kdDebug() << "RkOdeSolver: eps was negative, made it positive" << endl;
+      mEps = -a;
+   }
+   else if (a == 0.0)
+   {
+      mEps = 1e-5;              // a very arbitrary value
+      kdDebug() << "RkOdeSolver: eps == 0, set it to 1e-5" << endl;
+   }
+   else
+   {
+      mEps = a;
+   }
+}
+
+template<typename T>
+void
+RkOdeSolver<T>::Y(const std::valarray<T> &a)
+{
+   mY.resize(a.size());
+   mY = a;
+}
+
+
+// public member functions
+
+template<typename T>
+void
+RkOdeSolver<T>::integrate(const T &deltaX)
+{
+   if (deltaX == 0)
+   {
+      return;                   // nothing to integrate
+   }
+
+   // init dydx if uninitialised
+   if (mDydx.size() != mY.size())
+   {
+      mDydx.resize(mY.size());
+      mDydx = f(mX,mY);
+   }
+
+   static const unsigned int maxiter = 10000;
+   const T x2 = mX + deltaX;
+
+   unsigned int iter;
+   for (iter=0;
+        iter<maxiter && rkStepCheck(x2-mX) == false;
+        ++iter)
+   {}
+
+   if (iter>maxiter)
+   {
+      kdDebug() << "RkOdeSolver: More than " << maxiter
+                << " iterations in RkOdeSolver::integrate" << endl;
+      // TODO throw exeption here
+   }
+}
+
+
+// private member functions
+
+template<typename T>
+bool
+RkOdeSolver<T>::rkStepCheck(const T& dx_requested)
+{
+   static const T safety =  0.9;
+   static const T pshrnk = -0.25;
+   static const T pgrow  = -0.2;
+
+   // reduce step size by no more than a factor 10
+   static const T shrinkLimit = 0.1;
+   // enlarge step size by no more than a factor 5
+   static const T growthLimit = 5.0;
+   // errmax_sl = 6561.0
+   static const T errmax_sl = pow(shrinkLimit/safety, 1.0/pshrnk);
+   // errmax_gl = 1.89e-4
+   static const T errmax_gl = pow(growthLimit/safety, 1.0/pgrow);
+
+   static const unsigned int maxiter = 100;
+
+   if (dx_requested == 0)
+   {
+      return true;              // integration done
+   }
+
+   std::valarray<T> ytmp(mY.size());
+   std::valarray<T> yerr(mY.size());
+   std::valarray<T> t(mY.size());
+
+   bool stepSizeWasMaximal;
+   T dx;
+   if (std::abs(dx_requested) > mStep)
+   {
+      stepSizeWasMaximal = true;
+      dx = dx_requested>0 ? mStep : -mStep;
+   }
+   else
+   {
+      stepSizeWasMaximal = false;
+      dx = dx_requested;
+   }
+
+   // generic scaling factor
+   std::valarray<T> yscal = std::abs(mY) + std::abs(dx*mDydx) + 1e-15;
+
+   unsigned int iter = 0;
+   T errmax = 0;
+   do
+   {
+      if (errmax >= 1.0)
+      {
+         // reduce step size
+         dx *= errmax<errmax_sl ? safety * pow(errmax, pshrnk) : shrinkLimit;
+         stepSizeWasMaximal = true;
+         if (mX == mX + dx)
+         {
+            // stepsize below numerical resolution
+            kdDebug() << "RkOdeSolver: stepsize underflow in rkStepCheck"
+                      << endl;
+            // TODO throw exeption here
+            exit(0);
+         }
+         // new dx -> update scaling vector
+         yscal = std::abs(mY) + std::abs(dx*mDydx) + 1e-15;
+      }
+
+      ytmp   = rkStep(dx, yerr); // try to make a step forward
+      t      = std::abs(yerr/yscal); // calc the error vector
+      errmax = t.max()/mEps;    // calc the rel. maximal error
+      ++iter;
+   } while (iter < maxiter && errmax >= 1.0);
+
+   if (iter >= maxiter)
+   {
+      kdDebug() << "RkOdeSolver: too many iterations in rkStepCheck" << endl;
+      // TODO throw exeption here
+      exit(0);
+   }
+
+   if (stepSizeWasMaximal == true)
+   {
+      // estimate next step size if used step size was maximal
+      mStep  =
+         std::abs(dx)
+         * (errmax>errmax_gl ? safety * pow(errmax, pgrow) : growthLimit);
+   }
+   mX    += dx;                 // make step forward
+   mY     = ytmp;               // save new function values
+   mDydx  = f(mX,mY);           // and update derivatives
+
+   return std::abs(dx) < std::abs(dx_requested);
+}
+
+template<typename T>
+std::valarray<T>
+RkOdeSolver<T>::rkStep(const T& dx, std::valarray<T>& yerr) const
+{
+   static const T
+      a2=0.2, a3=0.3, a4=0.6, a5=1.0, a6=0.875,
+      b21=0.2,
+      b31=3.0/40.0,       b32=9.0/40.0,
+      b41=0.3,            b42=-0.9,        b43=1.2,
+      b51=-11.0/54.0,     b52=2.5,         b53=-70.0/27.0, b54=35.0/27.0,
+      b61=1631.0/55296.0, b62=175.0/512.0, b63=575.0/13824.0,
+      b64=44275.0/110592.0, b65=253.0/4096.0,
+      c1=37.0/378.0, c3=250.0/621.0, c4=125.0/594.0, c6=512.0/1771.0,
+      dc1=c1-2825.0/27648.0,  dc3=c3-18575.0/48384.0,
+      dc4=c4-13525.0/55296.0, dc5=-277.0/14336.0, dc6=c6-0.25;
+
+   std::valarray<T> ak2 = f(mX + a2*dx,
+                            mY + dx*b21*mDydx);             // 2. step
+   std::valarray<T> ak3 = f(mX + a3*dx,
+                            mY + dx*(b31*mDydx + b32*ak2)); // 3.step
+   std::valarray<T> ak4 = f(mX + a4*dx,
+                            mY + dx*(b41*mDydx + b42*ak2
+                                     + b43*ak3));           // 4.step
+   std::valarray<T> ak5 = f(mX + a5*dx,
+                            mY + dx*(b51*mDydx + b52*ak2
+                                     + b53*ak3 + b54*ak4)); // 5.step
+   std::valarray<T> ak6 = f(mX + a6*dx,
+                            mY + dx*(b61*mDydx + b62*ak2
+                                     + b63*ak3 + b64*ak4
+                                     + b65*ak5));           // 6.step
+   yerr      = dx*(dc1*mDydx + dc3*ak3 + dc4*ak4 + dc5*ak5 + dc6*ak6);
+   return mY + dx*( c1*mDydx +            c3*ak3 +  c4*ak4 +  c6*ak6);
+}
+
+
+// explicite instantiations
+//template RkOdeSolver<long double>;
+template class RkOdeSolver<double>;
+//template RkOdeSolver<float>;