//? tut2.cpp
//? C+- by Ulrich Mutze. Status of work 2022-05-06.
//? Copyright (c) 2022 Ulrich Mutze
//? contact: see contact-info at www.ulrichmutze.de
//?
//? This program is free software: you can redistribute it and/or
//? modify it under the terms of the GNU General Public License as
//? published by the Free Software Foundation, either version 3 of
//? the License, or (at your option) any later version.
//?
//? This program is distributed in the hope that it will be useful,
//? but WITHOUT ANY WARRANTY; without even the implied warranty of
//? MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//? GNU General Public License <http://www.gnu.org/licenses/> for
//? more details.


/*******************************************************************************
   Studying several integrators applied to the Kepler oscillator
*******************************************************************************/


#include <cpmbas1.h>
   // Actually, this is a collecion of header files.
#include <cpmvectormore.h>
#include <cpmrmatrix.h>
#include <cpmcmatrix.h>
#include <cpmrmateig.h>
#include <cpmvectorgraphics.h>
#include <type_traits>

using namespace std;
using namespace std::filesystem;

using namespace CpmRoot;
using namespace CpmRootX;
using namespace CpmSystem;
using namespace CpmArrays;
using namespace CpmGraphics;
using namespace CpmImaging;

using CpmLinAlg::R2;
using CpmGeo::Iv;
using CpmAlgorithms::solKepEqu;
using CpmAlgorithms::CyclicAlarm;
using CpmImaging::RED;
using CpmImaging::BROWN;
using CpmImaging::ORANGE;
using CpmImaging::YELLOW;
using CpmImaging::GREEN;
using CpmImaging::LIGHTGRAY;
using CpmImaging::LIGHTGREEN;
using CpmImaging::CYAN;
using CpmImaging::MAGENTA;
using CpmImaging::darkSky;
using CpmImaging::twiSky;
using CpmImaging::morSky;
using CpmImaging::spectral10;
using CpmImaging::spectral10DarkGround;
using CpmImaging::WRITECOLOR;
//using CpmDim2::Camera;
//using CpmDim2::Wrt;
//using CpmDim2::Path;
//using CpmDim2::Vec;
//using CpmDim2::AxVec;
//using CpmDim2::Spc;

namespace CpmKepler{

class KepOsc{ // 'Kepler oscillator'
// My name for the system defined as the radial part of the Kepler
// motion arround a fixed central mass. Physical units are used in which
// the moving mass, the gravitational potential at unit distance, and the
// angular momentum have value 1.
// see www.ulrichmutze.de/articles/99-271.pdf, Section 4
// This class implements the dynamics in a naive manner which is meant
// only as a pedagogical device. Classes that implement useful second order
// algorithms for dynamics will be derived from KepOsc.

protected:
// data

   R t_{0_R},x_,v_,dt_; // time, position, velocity, time step

// static functions
   static R f(R const& x){R xi=cpminv(x); return (xi-1_R)*xi*xi;}
      // force law

// member functions
   void init_( R eps, R nPerRev, R E, R t);
      // assumes 0<=eps<1, E arbitrary
      // initializes the common data

 public:
   virtual Word nameOf()const{ return "KepOsc";}
   KepOsc(){init_(0.,32.,0.,0.);}
   explicit KepOsc(R eps, R nPerRev=64., R E=0., R t=0.)
      {init_(eps,nPerRev,E,t);}
   virtual KepOsc* clone()const{ return new KepOsc(*this);}
   R t()const{ return t_;}
   R dt()const{ return dt_;}
   R x()const{ return x_;}
   R v()const{ return v_;}
   R eKin()const{ return v_*v_*0.5_R;}
      // kinetic energy
   R ePot()const{ R xi=cpminv(x_); return xi*(xi*0.5_R-1_R);}
      // potential energy
   static R pot(R const& x){ R xi=cpminv(x); return xi*(xi*0.5_R-1_R);}
   R eTot()const{ return eKin()+ePot();}
      //: e total

   R force()const{ return f(x_);}
      //: force
      // The superficially correct name 'for' is not correct since it is
      // a C++ key word.
   R_Vector orbDes()const;
      //: orbit descriptors
      // Returns a list of the usual orbit elements
      // a: major semi-axis, eps: numerical excentricity,
      // n: mean motion, E: excentric anomaly, M: mean anomaly,
      // tP: orbital period (revolution time)

   R2 ephem(R t, R acc)const;

   R_Vector ephem5(R t, R acc)const;

   R_Vector ephem7(R t, R acc)const;

   virtual void rev_(){ dt_=-dt_;}
      //: reversed
      // Notice name ending in '_' since it is a non-constant function
   void extEvl_(R t, R acc=1e-8);
      //: exact evolution as mutating operation

   KepOsc& eval_(R t, R acc=1e-8);

   void rev1_(){ v_=-v_;}
   void rev2_(){ dt_=-dt_;}
   void rev3_(){ dt_=-dt_;v_=-v_;}
      // reversion methods for experimentation

   virtual void step_()
      //: (evolution) step
      // Here only a implementation of the explicit
      // Euler rule. Notice mass is 1 so f(x_) is both force and
      // acceleration. It is meant as a definition of the law of motion
      // which in mathematical style would be written as the
      // first order ODE dx/dt=v, dv/dt=f(x)
   {
      R xOld=x_;
      x_+=v_*dt_;
      v_+= f(xOld)*dt_;
      t_+= dt_;
   }

   virtual void shf_(R dx, R dv){ x_+=dx; v_+=dv; }
      //: shift
      // virtual since in derived classes which use dependent state
      // variables one has to adjust to the new values of x_ and v_.

   void set_dt_(R dt){ dt_=dt;}

   X3<R,Iv,Iv> orbSize()const;
      //: orbit size
      // Returns the triplet tP (period), x_-range , and v_-range

   static V<R2> orb(R const& a, R const& eps, Z np=400_Z);
      //: orbit
      // Assume V<R2> vr2=KepOsc::orb(a,eps,np). For generic i write
      // vr2[i]=:R2(x,y). Then x and y are Cartesian coordinates of a point
      // which are taken with respect to the center point of the orbit ellipse
      // (which thus is taken as origin).


   static V<Color> colList(Graph const& g){
      if(g.isDark()) return spectral10DarkGround;
      else return spectral10;
   }

   virtual Color col(Graph const& g)const{
      return KepOsc::colList(g)[1];
   }
      //: color for visualization

   void mark(Graph& g, Z size=1)const
    //: mark
   {
      if (size==1){
         g.mark(R2(x_,v_),col(g));
         return;
      }
      g.mark(R2(x_,v_),size,col(g));
   }

   R2 phsPos()const{ return R2(x_,v_);}
      //: phase (space) position

   R per()const{ return orbSize().c1();}
      //:period

   virtual Word met()const{ return "Euler";}

   void diagnosis()const{
      cpmdata<<"Data from diagnosis:"<<endl;
      cpmdata<<"t_ = "<<t_<<endl;
      cpmdata<<"x_ = "<<x_<<endl;
      cpmdata<<"v_ = "<<v_<<endl;
      Word method=met();
      cpmdata<<"method = "<<-method<<endl;
   }
 };

void KepOsc::init_( R eps, R nPerRev, R E, R t)
// initializing x_ anv v_ from numerical excentricity and excentric
// anomaly. Further, initialization of the dependent quantities.
// Since we have E given, we have not to solve Kepler's equation and thus
// don't need acc
{
   Word loc("KepOsc::init_(R,R,R,R)");
   cpmassert(eps>=0,loc);
   cpmassert(eps<1,loc);
   t_=t;
   R ai=(1.-eps)*(1.+eps);
   R a=1./ai;
   x_=a*(1.-eps*cpmcos(E));
   v_=eps*cpmsqrt(a)*cpmsin(E)/x_;
   R tP=cpm2pi*cpmpow(a,1.5);
   dt_=tP/nPerRev;
   cpmdata<<"Data from KepOsc::init_()"<<endl;
   cpmdata<<"a = "<<a<<endl;
   cpmdata<<"E = "<<E<<endl;
   cpmdata<<"eps = "<<eps<<endl;
   cpmdata<<"nPerRev = "<<nPerRev<<endl;
   cpmdata<<"x_ = "<<x_<<endl;
   cpmdata<<"v_ = "<<v_<<endl;
   cpmdata<<"t_ = "<<t_<<endl;
   cpmdata<<"dt_ = "<<dt_<<endl;
}

R_Vector KepOsc::orbDes()const
{
   Word loc("KepOsc::orbDes()");
   Z mL=2;
   CPM_MA
   R e=eTot();
   if (e>=0){
      cpmcerr<<"We have e>=0, un-bound state. x_="<<x_<<" v_="<<v_<<endl;
      cpmdebug(ePot());
      cpmdebug(eKin());
   }
   cpmassert(e<0,loc); // we need a bound state
      // then also x>0.5
   R a=-0.5*cpminv(e);
   R eps=cpmsqrt(1.-1./a);
   R n=cpmpow(a,-1.5);
   R tP=cpm2pi/n;
   R z1=1-x_/a;
   R z2=x_*v_/cpmsqrt(a);
   C z(z1,z2);
   R E=z.arg();
   R M=E-eps*cpmsin(E);
   CPM_MZ
   return R_Vector{a,eps,n,E,M,tP};
}

R2 KepOsc::ephem(R t, R acc)const
// 'ephemeris generator'. For a given instance ko of class KepOsc and a
// time t we return the exact phase position (x,v) which the Kepler oscillator
// ko shows at time t (assuming 'exact dynamics' as implemented by solving Kepler's
// equation with accuracy acc). Notice that ko remains constant; it neither
// attains time t nor the quantities x and v which the function returns.
{
   Word loc("KepOsc::ephem(R,R)");
   Z mL=2;
   CPM_MA
   R a,eps,n,M,eaan;
   R_Vector y=orbDes();
   a=y[1];
   eps=y[2];
   n=y[3];
   M=y[5];
   eaan=eps*a*a*n;
   R Mt=M+n*(t-t_);
   R E=solKepEqu(Mt,eps,acc);
   R x=a*(1.-eps*cpmcos(E));
   cpmassert(x>=0.5,loc); // for safety, not needed
   R v=eaan*cpmsin(E)/x;
   CPM_MZ
   return R2(x,v);
}

R_Vector KepOsc::ephem5(R t, R acc)const
// 'ephemeris generator'. For a given instance ko of class KepOsc and a
// time t we return the exact phase position (x,v) which the Kepler oscillator
// ko shows at time t (assuming 'exact dynamics' as implemented by solving Kepler's
// equation with accuracy acc). Notice that ko remains constant; it neither
// attains time t nor the quantities x and v which the function returns.
{
   Word loc("KepOsc::ephem(R,R)");
   Z mL=2;
   CPM_MA
   R a,eps,n,M,eaan;
   R_Vector y=orbDes();
   a=y[1];
   eps=y[2];
   n=y[3];
   M=y[5];
   eaan=eps*a*a*n;
   R Mt=M+n*(t-t_);
   R E=solKepEqu(Mt,eps,acc);
   R x=a*(1.-eps*cpmcos(E));
   cpmassert(x>=0.5,loc); // for safety, not needed
   R v=eaan*cpmsin(E)/x;
   CPM_MZ
   return R_Vector{x,v,a,eps,E};
}

R_Vector KepOsc::ephem7(R t, R acc)const
// 'ephemeris generator'. For a given instance ko of class KepOsc and a
// time t we return the exact phase position (x,v) which the Kepler oscillator
// ko shows at time t (assuming 'exact dynamics' as implemented by solving Kepler's
// equation with accuracy acc). The point (x,y) is the exact position on
// the Kepler orbit at time t.Notice that ko remains constant; it neither
// attains time t nor the quantities x and v which the function returns.
{
   Word loc("KepOsc::ephem(R,R)");
   Z mL=2;
   CPM_MA
   R a,b,eps,n,M,eaan;
   R_Vector od=orbDes();
   a=od[1];
   eps=od[2];
   b=cpmsqrt(1_R-eps*eps);
   n=od[3];
   M=od[5];
   eaan=eps*a*a*n;
   R Mt=M+n*(t-t_);
   R E=solKepEqu(Mt,eps,acc);
   R x=a*(1.-eps*cpmcos(E));
   R y=b*cpmsin(E);
   cpmassert(x>=0.5,loc); // for safety, not needed
   R v=eaan*cpmsin(E)/x;
   CPM_MZ
   return R_Vector{x,y,v,a,b,eps,E};
}



X3<R,Iv,Iv> KepOsc::orbSize()const
{
   Z mL=2;
   Word loc("KepOsc::orbSize()");
   CPM_MA
   R_Vector y=orbDes();
   R a=y[1],eps=y[2],n=y[3],tP=y[6];
   R xMin=a*(1_R-eps);
   R xMax=a*(1_R+eps);
   R vMax=eps*a*a*n;
   CPM_MZ
   return X3<R,Iv,Iv>(tP,Iv(xMin,xMax),Iv(-vMax,vMax));
}

V<R2> KepOsc::orb(R const& a, R const& eps, Z np)
{
   R b=a*cpmsqrt(1_R-eps*eps);
   V<R2> res(np);
   R_Vector maL(np); // mean anomaly list
   maL[1]=0_R; // not really needed, sice all components are initialized as 0.
   maL[np]=cpm2pi;
   R xSun=a*eps;
   R delta=cpm2pi/(np-1_R);
   for (Z i=2;i<np;++i) maL[i]=maL[i-1]+delta;
   R acc= 1.e-12_R;
   for (Z i=1;i<=np;++i){
      R Mi=maL[i];
      R E=solKepEqu(Mi,eps,acc);
      res[i]=R2(xSun+a*(cpmcos(E)-eps),b*cpmsin(E));
   }
   return res;
}

// A version of KepOsc which does the exact evolution more economicly
// by transfering some work into the constructor and storing the result
// in suitable data elements.

class KepExc:public KepOsc{ // Exc: exact
   //Kepler Oscillator with exact dynamics (up to
   // numerical noise, of course).

// variables R t_,x_,v_,dt_; (time, position, velocity, time step)
// come from the base KepOsc.
// Additional variables in KepExt:
   R a_,eps_,n_,M_,eaan_,acc_{1e-12_R};
   void ini_() // class KepOsc brings a function 'init_' with it
   {
      R_Vector y=orbDes();
      a_=y[1];eps_=y[2];n_=y[3];M_=y[5];
      eaan_=eps_*a_*a_*n_; // do this once
      cpmcerr<<"a_ = "<<a_<<endl;
      cpmcerr<<"eps_ = "<<eps_<<endl;
      cpmcerr<<"n_ = "<<n_<<endl;
      cpmcerr<<"eaan_ = "<<eaan_<<endl;
      cpmcerr<<"M_ = "<<M_<<endl;
   }
public:
   KepExc():KepOsc(){ini_();}

   KepExc(R eps, R nPerRev, R E=0., R t=0.):
      KepOsc(eps,nPerRev,E,t) // by not calling ini_() here we spare a redundant
         // call to orbDes()
      {
         eps_=eps; M_=E-eps_*cpmsin(E);
         R nom=(1_R-eps_)*(1_R+eps_);
         a_=1_R/nom;
         n_=cpmpow(a_,-1.5_R); eaan_=eps_*a_*a_*n_;
      }

//    KepExc(R eps, R nPerRev, R E=0., R t=0.):
//      KepOsc(eps,nPerRev,E,t){ ini_();}  // ok, but more work than above!

   explicit KepExc(KepOsc const& ko):KepOsc(ko){ini_();}

   KepExc& setAcc_(R acc){acc_=acc;return *this;}
   R getAcc()const{ return acc_;}

   KepOsc* clone()const{ return new KepExc(*this);}

   void step_(){
      M_+=n_*dt_;
      R E=solKepEqu(M_,eps_,acc_);
      x_=a_*(1.-eps_*cpmcos(E));
      v_=eaan_*cpmsin(E)/x_;
      t_+=dt_; // a_,n_,eps_,eaan_, acc_, dt_ remain unchanged
   }

   void rev_(){ dt_=-dt_;}
      //: reversed

   void shf_(R dx, R dv){ x_+=dx; v_+=dv; ini_(); }

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[10];
   }
   Word met()const{ return "Exc";}
};

// Now the other integration methods are defined a derived classes with
// properly ovewritten virtual functions

// Runge Kutta integrator (2nd order)
// needs only the common data.

class KepRK2 :public KepOsc{ // RK2: Runge Kutta second order
   public:
   Word nameOf()const{ return "KepRK2";}
   KepRK2():KepOsc(){}
   explicit KepRK2(R eps, R nPerRev=64, R E=0):KepOsc(eps,nPerRev,E){}
   KepOsc* clone()const{ return new KepRK2(*this);}
      // all derived classes have to define clone in order to make my
      // polymorphic containers Pp<> and Vp<> work
   void step_()
   {
      R tau=dt_*0.5;
      R a=f(x_);
      R am=f(x_+v_*tau);
      R dv=am*dt_;
      R dx=(v_+a*tau)*dt_;
      x_+=dx;
      v_+=dv;
      t_+=dt_;
   }
  Color col(Graph const& g)const{
      return KepOsc::colList(g)[2];
  }
  Word met()const{ return "RK2";}
};

// Direct midpoint integrator (2nd order)
// needs only the common data.

class KepDMI :public KepOsc{ // DMI: direct midpoint integrator
   // also known as Stoermer Verlet integrator
   public:
   Word nameOf()const{ return "KepDMI";}
   KepDMI():KepOsc(){}
   explicit KepDMI(R eps, R nPerRev=64, R E=0):KepOsc(eps,nPerRev,E){}
   KepOsc* clone()const{ return new KepDMI(*this);}
      // all derived classes have to define clone in order to make my
      // polymorphic containers Pp<> and Vp<> work

   void step_()
   {
      R tau=dt_*0.5;
      t_+=tau;
      x_+=(v_*tau);
      v_+=f(x_)*dt_;
      x_+=(v_*tau);
      t_+=tau;
   }

    Color col(Graph const& g)const{
      return KepOsc::colList(g)[3];
   }
   Word met()const{ return "DMI";}
};

class KepALF :public KepOsc{ // ALF: asynchronous leap frog
   R w_,a_;
      // Asynchronous leap frog methods hold two states psi and phi. Here these
      // are psi=R2(x_,v_) and phi=R2(w_,a_).
public:
   Word nameOf()const{ return "KepALF";}
   KepALF():KepOsc(),w_(v_),a_(f(x_)){}
   explicit KepALF(R eps, R nPerRev=64, R E=0):
      KepOsc(eps,nPerRev,E),w_(v_),a_(f(x_)){}
   KepOsc* clone()const{ return new KepALF(*this);}

   void shf_(R dx, R dv)
   { x_+=dx; v_+=dv; w_=v_; a_=f(x_);}

   void step_()
   {
      R tau=dt_*0.5;
      t_ += tau;
      x_ += tau*w_;
      v_ += tau*a_;
      w_ = v_*2.-w_;
      a_ = f(x_)*2.-a_;
      v_ += tau*a_;
      x_ += tau*w_;
      t_ += tau;
   }

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[4];
   }
   Word met()const{ return "ALF";}
};

class KepDALF :public KepOsc{ // DALF: densified asynchronous leap frog
   R w_,a_;
   // Asynchronous leap frog methods hold two states psi and phi. Here these
   // are psi=R2(x_,v_) and phi=R2(w_,a_).
public:
   Word nameOf()const{ return "KepDALF";}
   KepDALF():KepOsc(),w_(v_),a_(f(x_)){}
   explicit KepDALF(R eps, R nPerRev=64, R E=0):
      KepOsc(eps,nPerRev,E),w_(v_),a_(f(x_)){}
   KepOsc* clone()const{ return new KepDALF(*this);}

   void shf_(R dx, R dv)
   { x_+=dx; v_+=dv; w_=v_; a_=f(x_);}

   void step_() // two ALF steps of length dt_/2
   // the third substep of the first step and the first substep of the
   // second step are merged to a single substep.
   {
      R tau=dt_*0.25_R;
      R tau2=dt_*0.5_R;
      t_ += tau;
      x_ += tau*w_;
      v_ += tau*a_;
      w_ = v_*2_R-w_;
      a_ = f(x_)*2_R-a_;
      t_ += tau2;
      x_ += tau2*w_;
      v_ += tau2*a_;
      w_ = v_*2_R-w_;
      a_ = f(x_)*2_R-a_;
      v_ += tau*a_;
      x_ += tau*w_;
      t_ += tau;
   }

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[5];
   }
   Word met()const{ return "DALF";}
};


class KepLF :public KepOsc{
// LF: leap frog
   R tp_,xp_,vp_; // p for previous
      // Leap frog methods hold two states, here these are R2(x_,v_) and
      // R2(xp_,vp_)
   void initLF_()
      // does the additional initializations
      // simplified version, we need not be more accurate than second
      // order
   {
      dt_*=2.;
      R h= -dt_*0.5;
      tp_=t_+h;
      R a=f(x_);
      R dv=a*h*0.5;
      xp_=x_+(v_+dv)*h; // here xp_ is previous x_
      vp_=v_+(a+f(x_))*h*0.5;
   }
public:
   Word nameOf()const{ return "KepLF";}
   KepLF():KepOsc(){initLF_();}
   explicit KepLF(R eps, R nPerRev=64, R E=0):
      KepOsc(eps,nPerRev,E){initLF_();}
   KepOsc* clone()const{ return new KepLF(*this);}

   void step_()
   {
      R tn=tp_+dt_;
      R vn=vp_+f(x_)*dt_;
      R xn=xp_+v_*dt_;
      tp_=t_;
      vp_=v_;
      xp_=x_;
      t_=tn;
      v_=vn;
      x_=xn;
   }

   void rev_(){
      R xMem=x_;
      R vMem=v_;
      R tMem=t_;
      step_();
      R tn=t();
      R xn=x();
      R vn=v();
      xp_=xn;
      vp_=vn;
      tp_=tn;
      x_=xMem;
      t_=tMem;
      v_=vMem;
      dt_=-dt_;
   }

   void shf_(R dx, R dv)
   { x_+=dx; v_+=dv; initLF_();}

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[6];
   }
   Word met()const{ return "LF";}
};

class KepADALF :public KepOsc{
// ADALF: averaged densified leap frog
   R w_,a_;
   // Asynchronous leap frog methods hold two states psi and phi. Here these
   // are psi=R2(x_,v_) and phi=R2(w_,a_).
public:
   Word nameOf()const{ return "KepADALF";}
   KepADALF():KepOsc(),w_(v_),a_(f(x_)){}
   explicit KepADALF(R eps, R nPerRev=64, R E=0):
      KepOsc(eps,nPerRev,E),w_(v_),a_(f(x_)){}
   KepOsc* clone()const{ return new KepADALF(*this);}

   void shf_(R dx, R dv)
   { x_+=dx; v_+=dv; w_=v_; a_=f(x_);}

   void step_()
   {
      R tau=dt_*0.25;
      R tau2=dt_*0.5;
      t_ += tau;
      x_ += tau*w_;
      v_ += tau*a_;
      w_ = v_*2.-w_;
      a_ = f(x_)*2.-a_;
      R wm=w_;
      R am=a_;
      t_ += tau2;
      x_ += tau2*w_;
      v_ += tau2*a_;
      w_ = v_*2.-w_;
      a_ = f(x_)*2.-a_;
      v_ += tau*a_;
      x_ += tau*w_;
      t_ += tau;
      a_=(a_+am)*0.5;
      w_=(w_+wm)*0.5;
   }

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[7];
   }
   Word met()const{ return "ADALF";}
};

// averaged asynchronous leap-frog (2nd order)
// not densified! New, change 2013/14.
// This is an inferior method that I deviced and tested at the
// very beginning of investigation of second order integrators
// for the time-dependend Schroedinger equation. I discarded the
// method and forgot it. Finally I re-invented it as the simplest
// version of the averaging method.

class KepAALF :public KepOsc{
// AALF: averaged asynchronous leapfrog
   R w_,a_;
   // Asynchronous leap frog methods hold two states psi and phi. Here these
   // are psi=R2(x_,v_) and phi=R2(w_,a_).
public:
   Word nameOf()const{ return "KepAALF";}
   KepAALF():KepOsc(),w_(v_),a_(f(x_)){}
   explicit KepAALF(R eps, R nPerRev=64, R E=0):
      KepOsc(eps,nPerRev,E),w_(v_),a_(f(x_)){}
   KepOsc* clone()const{ return new KepAALF(*this);}

   void shf_(R dx, R dv)
   { x_+=dx; v_+=dv; w_=v_; a_=f(x_);}

   void step_()
   {
      R tau=dt_*0.5;
      R xm=x_;
      R vm=v_;
      t_ += tau;
      x_ += tau*w_;
      v_ += tau*a_;
      w_ = v_;
      a_ = f(x_);
      v_ = vm + dt_*a_;
      x_ = xm + dt_*w_;
      t_ += tau;
   }

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[8];
   }
   Word met()const{ return "AALF";}
};

// This is the first order symplectic integrator for Hamiltonian systems in
// Sanz-Serna, Calvo p. 130(-1). It is computationally even simpler than
// the Euler method since it does not use a temporary memory variable.
// further it is reversible and symplectic. It is of first order just as
// Euler's method.
class KepHAM :public KepOsc{
// HAM: Hamilton
public:
   Word nameOf()const{ return "KepHAM";}
   KepHAM():KepOsc(){}
   explicit KepHAM(R eps, R nPerRev=64, R E=0):
      KepOsc(eps,nPerRev,E){}
   KepOsc* clone()const{ return new KepHAM(*this);}

   void shf_(R dx, R dv)
   { x_+=dx; v_+=dv;}

   void step_() // even simpler than Euler
   // arguably the simplest integrator for an ODE of second order.
   // Is symplectic!!!!
   {
      t_ += dt_;
      x_ += dt_*v_;
      v_ += dt_*f(x_); // this calls f with updated x_
   }

   Color col(Graph const& g)const{
      return KepOsc::colList(g)[9];
   }
   Word met()const{ return "HAM";}
};

// The file input select integration methods simply by integers; here we
// select the corresponding classes.
// Presently there are 9 integration methods.

Pp<KepOsc> singlet(Z method, R eps, R nPerRev, R E=0.)
{
   Pp<KepOsc> res;
   if (method==1) res=KepOsc(eps,nPerRev,E);
   else if (method==2) res=KepRK2(eps,nPerRev,E);
   else if (method==3) res=KepDMI(eps,nPerRev,E);
   else if (method==4) res=KepALF(eps,nPerRev,E);
   else if (method==5) res=KepDALF(eps,nPerRev,E);
   else if (method==6) res=KepLF(eps,nPerRev,E);
   else if (method==7) res=KepADALF(eps,nPerRev,E);
   else if (method==8)  res=KepAALF(eps,nPerRev,E);
   else if (method==9)  res=KepHAM(eps,nPerRev,E);
   else if (method==10)  res=KepExc(eps,nPerRev,E);
   return res;
}

Vp<KepOsc> generator(S<Z> const& methods, R eps, R nPerRev, R E=0.)
{
   Z d=methods.car();
   Vp<KepOsc> res(d);
   Z i=0;
   S<Z> m{methods}; // to avoid the dubious error on possible
   // conflict with constantness of method.
   Z meti;
   for (auto q=m.cbegin(); q!=m.cend(); ++q){
      meti=*q;
      i++;
      if (meti==1) res[i]=KepOsc(eps,nPerRev,E);
      else if (meti==2) res[i]=KepRK2(eps,nPerRev,E);
      else if (meti==3) res[i]=KepDMI(eps,nPerRev,E);
      else if (meti==4) res[i]=KepALF(eps,nPerRev,E);
      else if (meti==5) res[i]=KepDALF(eps,nPerRev,E);
      else if (meti==6) res[i]=KepLF(eps,nPerRev,E);
      else if (meti==7) res[i]=KepADALF(eps,nPerRev,E);
      else if (meti==8) res[i]=KepAALF(eps,nPerRev,E);
      else if (meti==9) res[i]=KepHAM(eps,nPerRev,E);
      else if (meti==10) res[i]=KepExc(eps,nPerRev,E);
   }
   return res;
}

//////////////////////////////////////////////////////////////////////////

void cpmplot(V<V<R2> > const& dat, Word const& fileName)
{
   Word loc("cpmplot(V<V<R2>>,Word)");
   Z mL=2;
   CPM_MA
   if (fileName.isVoid()){
      CPM_MZ
      return;
   }
   Z i,j,md=dat.dim();
   V<Z> vd(md);
   for (i=1;i<=md;++i){
      vd[i]=dat[i].dim();
   }
   S<Z> sd(vd);
   cpmassert(sd.car()==1,loc);
   Z nd=sd[1];
   OFileStream ofs(fileName);
   ofs()<<std::setprecision(CpmRoot::wrtPrc);
      // otherwise multiple precision mode would not be flexibly be used
   for (i=1;i<=nd;++i){
      ofs()<<endl;
      for (j=1;j<=md;++j){
         R2 y=dat[j][i];
         ofs()<<y[1]<<" "<<y[2]<<"    ";
      }
   }
   CPM_MZ
}

void app1(void) // orbit in phase space (v = momentum since m=1)
// simultaneously for a selected list of integration methods
{
   Message::subProgram="app1";
   Word loc("app1");
   Color cb=darkSky; // night sky blue. Each task appi
   //Color cb=morSky;
   RecordHandler rch("app1.ini");
      // This says that we expect data input to be available in file
      // app1.ini in a diretory which is determined by entries in
      // cpmsystemdependencies.h and cpmconfig.ini.
   Word sec="data"; // names 'rch' and 'sec' are to be employed
      // for the macro cpmrh to work
   R eps=0.25, fac=1.2, nPerRev=16, nRev=2, tWaitFinal=1;
   Z displayPeriod=1, pointSize=1;
   Word imageName,plotDataName;
   V<Z> methods;
// getting quantities from section 'data' of file app1.ini
   cpmrh(eps); // rh: read handler
   cpmrh(nPerRev);
   cpmrh(nRev);
   cpmrh(fac);
   cpmrh(tWaitFinal);
   cpmrhf(imageName);
   cpmrhf(plotDataName);
   cpmrh(displayPeriod);
   cpmrh(pointSize);
      // now the program has all its data
   cpmrh(methods);
   S<Z> met(methods); // eliminates doublicates
   Vp<KepOsc> vk=generator(met,eps,nPerRev);
   KepOsc koMem(eps,nPerRev);
   Z md=met.car();
   X3<R,Iv,Iv> os=vk(1).orbSize(); // cartesian product templates
      // are defined from X2 to X8
   Iv ivx=os.c2(); // Iv = interval is a valuable type
   Iv ivy=os.c3();
   ivx*=fac; // changing the size of an interval, while holding the center
      // fixed
   ivy*=fac; // same factor for the other direction
   Z nTot=cpmrnd(nPerRev*nRev);
   V<V<R2> > vvr(md,V<R2>(nTot)); // memory for all phase space positions
      // during run
   Graph gr(cb); // basic graphical 'window to the world'
   //Graph gr;
   gr.setXY(ivx,ivy);
   cpmdata<<"ivx = "<<ivx<<endl;
   cpmdata<<"ivy = "<<ivy<<endl;
   cpmdata<<"eps = "<<eps<<endl;
      // this associates physical sizes (in space) with the
      // 'pixeled' rectangle g
   R nInv=1./nTot; // for use in cpmprogress
   CyclicAlarm cy(displayPeriod); // device for triggering periodic
      // actions
   bool do_it=false;
   for (Z i=1;i<=nTot;++i){
      if (cy()) do_it=true; else do_it=false;
      for (Z j=vk.b();j<=vk.e();++j){
         if (do_it){
            Z ps=(j==1||j==9||j==10 ? pointSize : 1);
            vk(j).mark(gr,ps);
            cpmprogress("loop",i*nInv);
            gr.vis();
         }
         vk(j).step_();
         //cpmwait(0.1,2);
      }
      if (do_it){ // cy() becomes 'true' periodically
         gr.vis(); // transfers the whole memory frame
            // to screen. This expensive action should probably not be
            // orderd for each integration step. Object cy allows
            // to achieve this.
            // Progress indicator on pane 1 of he status bar.
      }
      if (i==nTot){
         V<Word> title{"A single phasespace orbit",
            "represented by various integrators"};
         V<Color> col{WHITE,WHITE};
         R down=0.6;
         gr.setText(title,col,down);
         gr.vis();
      }
   }
   Word img=gr.wrt(imageName); // writes the final frame to a ppm image file
      // if imageName is not the void word.
   CpmSystem::ppm2png(img);
   cpmplot(vvr,Message::extOut(plotDataName));
      // does nothing if the name is void
   cpmwait(tWaitFinal,2); // allows us to inspect the last frame
      // for the time (in seconds) put in here. The second argument
      // says that on pane 2 of the status bar there will be a
      // 'countdown' indication of the time left for inspection.
}

// Building block for app2(void). Not a stand-alone app which should be clear
// as there many arguments.
void app2a(S<Z> const& met, V<R> const& veps, V<R> const& vnPerRev,
   V<R> const& vE, Z const& nTot,  B const& backEvl, B const& noRefToExc,
   R const& acc,
   Z const& nFrames, Word const& imageName,
   Word const& movieName,
   R const& tWaitFinal, Z i, Z j, Z k, Z frameCount, B silent,
   V<R2>& vr2, Graph& gr)
{
   Vp<KepOsc> vk=generator(met,veps[j],vnPerRev[k],vE[i]);
  // if (!silent) gr.clr_(cb);
   KepOsc k0=vk(1);
   R t0=vk(1).t();
   // This defines the initial state.
   R2 p0=k0.phsPos();
   V<V<R2> > vvr(vk.dim(),V<R2>(nTot));
   Z nTotF=nTot*nFrames;
   R nInvF=1./nTotF;
   for (Z i=1;i<=nTot;++i){ // no graphical actions
      // within the loop. Terribly fast!!!
      KepOsc k1Act=vk(1); // conversion to base class is OK
         // since we call only functions which don't get re-defined
         // in derived classes
      R tAct=k1Act.t(); // the states are aware of the time to which they
         // belong: Schroedinger picture
      if (!backEvl){
         R2 posExt=k0.ephem(tAct,acc); // storing the exact phase space
                 // position
         for (Z j=vk.b();j<=vk.e();++j){
            R2 diff=noRefToExc ? vk(j).phsPos() : vk(j).phsPos()-posExt;
            if (silent) vr2<<diff;
            vvr[j][i]=diff;
         }
      }
      else{
         for (Z j=vk.b();j<=vk.e();++j){
            KepOsc vkj(vk(j));
            R2 pj=vkj.ephem(t0,acc);
            R2 diff=pj-p0;
            if (silent) vr2<<diff;
            vvr[j][i]=diff;
         }
      }
      for (Z j=vk.b();j<=vk.e();++j){
         vk(j).step_();
      }
      if (silent) cpmprogress("loop",(1_R*i*nFrames)*nInvF,1,2);
         // only 2 digits for progress indication to avoid
         // 'much adoo about nothing'
   }
   if (!silent){
      Z nm=met.car();
      V<Color> vc(nm);
      for (Z j=1;j<=nm;++j) vc[j]=vk(j).col(gr);
      gr.mark(vvr,vc); // powerful method for graphical
      // representation of polygons;
      gr.vis(true);
   }
}

void app2(void) // method error of the two integration methods
// represented by means of interaction picture dynamics. For
// applying this terminology, we interpret the dynamics as defined by
// the numerical methods as resulting from the exact Kepler dynamics
// by 'interaction'. This interaction picture motion is much slower for
// the more accurate direct midpoint method. The Runge Kutta motion
// finally explodes and reaches non-bound states for which the exact
// evolution back is no longer defined (by the method KepOsc::extEvl_).
// Therefore the option is implemented to deal only with the stable
// direct midpoint method.
{
   Message::subProgram="app2";
   Word loc("app2");
   Color cb=twiSky; // background
   RecordHandler rch("app2.ini");
   Word sec="data";
   V<R> El=V<R>(2),epsl=V<R>(2), nPerRevl=V<R>(2);
   R nRev=2., acc=1e-8, tWaitIntermediate=6,tWaitFinal=1;
   Z i,j,nFrames=100;
   R frameRate=40;
   Word imageName,movieName;
   V<Z> methods;
   B backEvl{true};
   B noRefToExc{false};
   cpmrhf(El);
   cpmrh(epsl);
   cpmrh(nPerRevl);
   cpmrh(nRev);
   cpmrh(nFrames);
   cpmrh(frameRate);
   cpmrh(acc);
   cpmrh(tWaitIntermediate);
   cpmrh(tWaitFinal);
   cpmrh(imageName);
   cpmrhf(movieName);
   cpmrh(methods);
   cpmrhf(backEvl);
   cpmrhf(noRefToExc);

   S<Z> met(methods); // eliminates doublicates
   Z nm=met.dim();
   cpmassert(nm>0,"app2()");
   V<R> vE,veps,vnPerRev;

   vE = (El[1]==El[2] || nFrames==0) ?
   V<R>{El[1]} : Iv(El[1],El[2]).cen(nFrames);

   veps = (epsl[1]==epsl[2] || nFrames==0) ?
      V<R>{epsl[1]} : Iv(epsl[1],epsl[2]).cen(nFrames);

   vnPerRev = (nPerRevl[1]==nPerRevl[2] || nFrames==0) ?
      V<R>{nPerRevl[1]} : Iv(nPerRevl[1],nPerRevl[2]).cen(nFrames);

   // First task
   V<R2> vr2;
   Graph gr; // white background
   B silent=B(1);
      // the only task is to set vr2 which sets the graphics window in the
      // second task
   Z frameCount=0;
   for (Z i=vE.b();i<=vE.e();++i){
      for (Z j=veps.b();j<=veps.e();++j){
         for (Z k=vnPerRev.b();k<=vnPerRev.e();++k){
            frameCount++;
            Z nTot=nRev*vnPerRev[k];
            app2a(met,veps,vnPerRev,vE,nTot,backEvl,noRefToExc,acc,nFrames,
            imageName,movieName,tWaitFinal,i,j,k,frameCount,silent,vr2,gr);
         }
      }
   }

   Z frameCountMem=frameCount;
   cpmdata<<"frameCountMem = "<<frameCountMem<<endl;
   R fcmi=1_R/frameCountMem;
   cpmdata<<"fcmi = "<<fcmi<<endl;

   gr.mark(vr2); // sets the scale according to vr2
   gr.vis();
   cpmwait(tWaitIntermediate,2);

   // second task
   frameCount=0;
   silent=B(0);
   gr.clr_(cb); // setting a dark background
   gr.setGridColor(DARKGRAY); // should not too conspicuous
   gr.setGridTextColor(WHITE); // should be clearly readable
      // controlling the appearance of movie frames

   for (Z i=vE.b();i<=vE.e();++i){
      for (Z j=veps.b();j<=veps.e();++j){
         for (Z k=vnPerRev.b();k<=vnPerRev.e();++k){
            frameCount++;
            Z nTot=nRev*vnPerRev[k];
            app2a(met,veps,vnPerRev,vE,nTot,backEvl,noRefToExc,acc,nFrames,
            imageName,movieName,tWaitFinal,i,j,k,frameCount,silent,vr2,gr);
            cpmprogress("loop",frameCount*fcmi,1,2);
         }
      }
   }

   bool makeMovie=movieName.dim()>0;
   if (makeMovie){ // creating video
      path dirMem = current_path();
      current_path(-Message::getOutDir2());
      path dirAct = current_path();
      cpmcerr<<"creating movie at "<<dirAct<<endl;
      ppm2mp4(movieName,frameRate);
      current_path(dirMem);
   }
}

void app3(void) // interaction picture dynamics for not a single
// point as an initial condition but for an ensemble of
// phase space points which form a closed curve. The deformation
// of this closed curve by interaction picture dynamics is monitored.
// Here the field of view in phase space follows the moving states.
// Otherwise they would slowly drift out. For the asynchronous
// leapfrog method (method 4) it was observed that the loop
// deformation can lead to crossing over. This puzzled me for a while
// since this cannot happen with mappings of the phase space into
// the phase space. However for the asynchronous leapfrog method
// the statespace (which is properly mapped) has additional degrees
// of freedom, which are not shown in the phase space curves. This
// resolves what seemed to be a logical contradiction. The details of
// the process are not clear to me. Finally the time evolution of the
// enclosed phase area is shown.
{
   Message::subProgram="app3";
   Word loc("app3");
   Color cb=darkSky; // background
   RecordHandler rch("app3.ini");
   Word sec="data";
   R eps=0.25, fac=1.2, nPerRev=16, acc=1e-8, E=0;
   R facX=1, facV=1, tWaitFinal1=1, tWaitFinal2=1, tWaitFrame=0.1;
   Z nTot=1000, nK=12, displayPeriod=1, method=1;
   B fileOnShow, backEvl;
   V<Z> filingLimits;

   cpmrh(eps);
   cpmrh(E);
   cpmrh(fac);
   cpmrh(nPerRev);
   cpmrh(nTot);
   cpmrh(backEvl);
   cpmrh(acc);
   cpmrh(tWaitFinal1);
   cpmrh(tWaitFinal2);
   cpmrh(tWaitFrame);
   cpmrh(nK);
   cpmrh(facX);
   cpmrh(facV);
   cpmrh(displayPeriod);
   cpmrh(method);
   cpmrhf(fileOnShow);
   cpmrhf(filingLimits);

   IvZ imgInterval; // void
   if (filingLimits.dim()>=2){
      imgInterval=IvZ(filingLimits[1],filingLimits[2]);
   }
   // It turned out to be necessary for being able to document the
   // phenomena that occur in runs to write screen images to file.
   // Since filing an image takes time, it is desirable to have a way
   // to file only a controlled subset of all screen frames.
   // Letting the program run without any filing one can make notes
   // of the frame numbers that contain the interesting events
   // (such as crossing over of loops). In a second run one then may
   // write these frames to file by suitable setting filingLimits.

   Pp<KepOsc> ko=singlet(method,eps,nPerRev,E);

   R x0=ko().x();
   R v0=ko().v();
   R t0=ko().t();
   R dt0=ko().dt();
   X3<R,Iv,Iv> os=ko().orbSize();

   Pp<KepOsc> kt=ko; // copy constructor works for Pp<>

   kt().step_();
   R xt=kt().x();
   R vt=kt().v();
   R dx=xt-x0;
   R dv=vt-v0;
   R r=R2(dx,dv).abs();
   R rx=r*facX;
   R rv=r*facV;
   Vp<KepOsc> vk(nK);
   vk.set(ko());
   R phi=0;
   R dphi=cpm2pi/(nK-1);
   R_Vector xx(nK);
   R_Vector vv(nK);

   Z i,j;
   for (i=vk.b();i<=vk.e();++i){
      vk(i).shf_(rx*cpmcos(phi),rv*cpmsin(phi));
      xx[i]=vk(i).x();
      vv[i]=vk(i).v();
      phi+=dphi;
   }
   R_Vector xx0=xx;
   R_Vector vv0=vv;

   R rf=r*fac;
   Iv ivx(-rf,rf);
   Iv ivy=ivx;
   Graph g(cb);
   g.setXY(ivx,ivy);
   R nInv=1./nTot;
   Z nPrev=nTot-1;
   Z imgCounter=0;
   CyclicAlarm cy(displayPeriod);
   R_Vector areas(nTot);
   R_Vector times(nTot,t0);
   for (i=1;i<=nTot;++i){
      if (i>1) times[i]=times[i-1]+dt0;
      for (j=vk.b();j<=vk.e();++j){ // sum over the contour
         vk(j).step_();
         R2 xvj= backEvl ? vk(j).ephem(t0,acc) : vk(j).phsPos();
         xx[j]=xvj.x1;
         vv[j]=xvj.x2;
      }
      R_Vector xxs=xx.subMean(); // subtracting the mean
      R_Vector vvs=vv.subMean();
      areas[i]=area(xxs,vvs); // CpmArrays::area(R_Vector,R_Vector)
         // gives the area enclosed by a polygon
      if (i<nPrev){
         if (cy()){
            g.draw(xxs,vvs,ko().col(g));
            g.setText("frame "&cpm(imgCounter++),0.05,0,ko().col(g));
            bool fileAct=false;
            if (fileOnShow) fileAct=imgInterval.hasElm(imgCounter);
            g.vis(fileAct);
            cpmprogress("loop",i*nInv);
            cpmwait(tWaitFrame,2);
            g.clr_(cb);
         }
      }
      else if (i==nPrev) g.clr_(cb); // that only a single
         // drawing is on the final display
      else{
         R redFac=0.5; // color of initial configuration should
            // be a bit less bright
         Color ci=ko().col(g)*redFac;
         g.draw(xxs,vvs,ko().col(g));
         g.draw(xx0.subMean(),vv0.subMean(),ci); // the final
         // image should also show the initial configuration
         cpmprogress("loop",i*nInv);
      }
   }
   V<Word> title=!backEvl ? V<Word>{"enclosed phase space area"} :
   V<Word>{"enclosed area with backevolution by exact dynamics"};
   V<Color> col{WHITE,WHITE};
   R down=0.66;
   g.setText(title,col,down);
   g.vis();
   g.wrt("image1");
   cpmwait(tWaitFinal1,2);
   if (tWaitFinal2>0){ // shows how the enclosed area evolves over time
      Graph gg(GREEN);
      gg.setXY(times,areas);
      gg.setWithOrigin(0);
      gg.mark(times,areas,Color(WHITE));
      gg.addText("Time evolution of enclosed area");
         // don't use setText here: this would purge the legend concerning
         // the x and y ranges
      gg.vis();
      gg.wrt("image2");
      cpmwait(tWaitFinal2,2);
   }
}

R cpmlog10s(R x)
{
   R y=cpmabs(x);
   if (y<1e-20) return -20;
   else return cpmlog10(y);
}

void app4(void)
// Analysis for order. Normal order turns out to be 2.
// However, if order is inferred from a run for which nPerRev is integer
// and nTot is an integer multiple of nPerRev one gets order 4 (and 2.87
// for method RK2). If one changes nTot only by a single step one falls
// to the normal value 2 of the order. This resonance-like behavior of
// order here observed in the case of a strictly periodic motion is not
// yet understood. Since a numerical method can be interpreted as a
// perturbation of the exact method, the KAM theory applies if the orbits
// of the exact method are periodic, just as in our case.
// Notice that so far the error was only recorded at the end of the
// trajectory. The new control 'endOnly' allows to form a error sum over
// the whole trajectory. This never shows orders significantly different
// from 2.

{
   Message::subProgram="app4";
   Word loc("app4");
   RecordHandler rch("app4.ini");
   Word sec="data"; // names 'rch' and 'sec' are to be employed
      // for the macro cpmrh to work
   R eps=0.25, tWaitFinal=1, nPerRev=16, acc=1e-8;
   Z i,j, nTot=16, nDoubling=2, method=2;
   B endOnly;
 // getting quantities from section 'data' of file app4.ini
   cpmrh(eps); // rh: read handler
   cpmrh(nPerRev);
   cpmrh(nTot);
   cpmrh(nDoubling);
   cpmrhf(endOnly);
   cpmrhf(acc);
   cpmrh(tWaitFinal);
   cpmrh(method);
   cpmrhf(cpmverbose);
   Z nData=nDoubling+1;
   R_Vector x(nData);
   R_Vector y(nData);
   Z sp=1; // The sampling period is taken as the number of steps in the
      // first run. Doubling the step number will  n o t  increase the
      // number of samples per revolution. Notice that for endOnly = 1
      // only the error at the end-point of the trajectory is taken into
      // account.
   Word wm; // word for method
   KepOsc koExt(eps);
   for (i=1;i<=nData;++i){
      Pp<KepOsc> ko=singlet(method,eps,nPerRev);
      if (i==1) wm=ko().met();
      CyclicAlarm cy(sp);
      R ei=0;
      for (j=1;j<=nTot;++j){
         ko().step_();
         if (endOnly ? j==nTot : cy()){
            R ts=ko().t();
            R xExt=koExt.ephem(ts,acc).x1;
            R xs=ko().x();
            ei+=cpmabs(xExt-xs);
         }
      }
      nTot*=2; // implementing halfing the time-step
      nPerRev*=2.0; // implementing halfing the time-step
      sp*=2; // this lets the total number of sampling points constant.
      x[i]=cpmlog10(ko().dt()); // log timestep
      y[i]=cpmlog10s(ei);
         // log absolute position error at tTot or errorsum at more
         // points
      cpmdata<<"x["<<i<<"] = "<<x[i]<<" y["<<i<<"] = "<<y[i]<<endl;
   }
   Z deg=1; // find a straight line
   R_Vector res=polyFit(x,y,deg,0.,1.,false);
   R order=cpmabs(res[4]);
   Graph g;
   g.setWithOrigin(0);
   g.addText("order="&cpm(order)&", method="&wm);
   g.mark(x,y);
   g.vis();
   g.wrt("image");
   cpmwait(tWaitFinal,2);
}

void app5(void)
// Collapsing error curves for different stepsizes by scaling them with
// an input-value for the order.
{
   Message::subProgram="app5";
   Word loc("app5");
   RecordHandler rch("app5.ini");
      // This says that we expect data input to be available in file
      // app5.ini in a diretory which is determined by entries in
      // cpmsystemdependencies.h and cpmconfig.ini.
   Word sec="data"; // names 'rch' and 'sec' are to be employed
      // for the macro cpmrh to work
   R eps=0.25, tWaitFinal=1, nPerRev=16, order=2, acc=1e-12;
   Z i,j, nTot=16, nDoubling=2, method=2;
// getting quantities from section 'data' of file app5.ini
   cpmrh(eps); // rh: read handler
   cpmrh(nPerRev);
   cpmrh(nTot);
   cpmrh(acc);
   cpmrh(nDoubling);
   cpmrh(tWaitFinal);
   cpmrhf(order);
   cpmrh(method);
   Z nData=nDoubling+1_Z;
   V< V<R2> > res(nData);
   KepOsc koExt(eps);
   R fac=1.;
   Word wm;
   for (i=1;i<=nData;++i){
      Pp<KepOsc> ko=singlet(method,eps,nPerRev);
      if (i==1) wm=ko().met();
      V<R2> resi(nTot);
      for (j=1;j<=nTot;++j){
         ko().step_();
         R tj=ko().t();
         R2 dpos=ko().phsPos()-koExt.ephem(tj,acc);
         R dphs=dpos.abs();
         resi[j]=R2(tj,dphs*fac);
      }
      res[i]=resi;
      nTot*=2_Z;
      nPerRev*=2_R;
      fac*=cpmpow(2_R,order);
   }
   Graph g;
   g.setWithOrigin(0);
   g.addText("method="&wm);
   g.mark(res,rainBow,1.1_R,B(1));
   g.vis();
   g.wrt("image");
   cpmwait(tWaitFinal,2);
}

void app6()
// Showing the time evolution of energy error for the integrators.
{
   Message::subProgram="app6";
   Word loc("app6");

   RecordHandler rch("app6.ini");
      // This says that we expect data input to be available in file
      // app1.ini in a diretory which is determined by entries in
      // cpmsystemdependencies.h and cpmconfig.ini.
   Word sec="data"; // names 'rch' and 'sec' are to be employed
      // for the macro cpmrh to work
   R eps=0.25, fac=1.2, nPerRev=16, tWaitFinal=1, E0=0.,acc=1.e-12;
   Z nTot=1000;
   V<Z> methods;
// getting quantities from section 'data' of file app1.ini
   cpmrh(eps); // rh: read handler
   cpmrh(nPerRev);
   cpmrh(nTot);
   cpmrh(tWaitFinal);
   cpmrh(methods);
   cpmrh(acc);
   cpmrh(E0);
   R t0=0.;
   S<Z> met(methods); // eliminates doublicates
   Vp<KepOsc> vk=generator(met,eps,nPerRev,E0);
   Z md=met.car();

   KepOsc ko(eps,nPerRev,E0,t0);
   R dt=ko.dt();
   R_Vector tL(nTot);
   for (Z i=1;i<=nTot;++i) tL[i]=(i-1)*dt; //tL[1]=0
   R_Vector xEx(nTot);
   R_Vector vEx(nTot);
   for (Z i=1;i<=nTot;++i){
      R2 xv=ko.ephem(tL[i],acc);
      xEx[i]=xv[1];
      vEx[i]=xv[2];
   }

   Graph g; // basic graphical 'window to the world'

   g.setAutoScale(true);
   g.setWithOrigin(1);

   R nInv=1./nTot; // for use in cpmprogress

   R_Matrix mat(md,nTot);
   V<Color> vc(md);
   for (Z i=1;i<=nTot;++i){
      for (Z j=1;j<=md;++j){
         if (i==1){
            vc[j]=vk(j).col(g);
         }
         mat[j][i]=cpmsqrt(cpmsqr(vk(j).x()-xEx[i])+cpmsqr(vk(j).v()-vEx[i]));
         vk(j).step_();
      }
   }
   g.markTr(tL,mat,vc);
   //g.wrt("image1");
   g.vis();// writes the final frame to a ppm image file
   g.wrt("image1");

   cpmwait(tWaitFinal,2); // allows us to inspect the last frame
      // for the time (in seconds) put in here. The second argument
      // says that on pane 2 of the status bar there will be a
      // 'countdown' indication of the time left for inspection.
}

//////////////////////////////// app7//////////////////////////////////////////

// Deviation from energy conservation of the various integrators for the
// Kepler oscillator.

void app7(void) // orbit in phase space (v = momentum since m=1)
// simultaneously for a selected list of integration methods

// Result:
{
   Message::subProgram="app7";
   Word loc("app7");

   RecordHandler rch("app7.ini");
      // This says that we expect data input to be available in file
      // app1.ini in a diretory which is determined by entries in
      // cpmsystemdependencies.h and cpmconfig.ini.
   Word sec="data"; // names 'rch' and 'sec' are to be employed
      // for the macro cpmrh to work
   R eps=0.25, fac=1.2, nPerRev=16, tWaitFinal=1;
   Z nTot=1000;
   V<Z> methods;
// getting quantities from section 'data' of file app1.ini
   cpmrh(eps); // rh: read handler
   cpmrh(nPerRev);
   cpmrh(nTot);
   cpmrh(tWaitFinal);
   cpmrh(methods);
   S<Z> met(methods); // eliminates doublicates
   Vp<KepOsc> vk=generator(met,eps,nPerRev);
   Z md=met.dim();

   Graph g; // basic graphical 'window to the world'

   g.setAutoScale(true);
   g.setText("rel. energy error");

   R nInv=1./nTot; // for use in cpmprogress

   R_Vector vec(nTot);
   R_Matrix mat(md,nTot);
   R e0;
   R e0Inv = e0==0_R ? 1_R : 1_R/e0;
   V<Color> vc(md);
   for (Z i=1;i<=nTot;++i){
      for (Z j=1;j<=md;++j){
         if (i==1){
            vc[j]=vk(j).col(g);
         }
         if (j==1) vec[i]=vk(1).t();
         if (j==1 && i==1) e0=vk(j).eTot();
         mat[j][i]=(vk(j).eTot()-e0)*e0Inv;
            // This is the relative deviation of total energy from its
            // initial value (called 'error' since constancy of energy is
            // the correct behaviour).
         vk(j).step_();
      }
   }
   g.markTr(vec,mat,vc);
   g.vis();
   g.wrt("image");// writes the final frame to a ppm and a png image file

   cpmwait(tWaitFinal,2); // allows us to inspect the last frame
      // for the time (in seconds) put in here. The second argument
      // says that on pane 2 of the status bar there will be a
      // 'countdown' indication of the time left for inspection.
}

//////////////////////////////// app8 //////////////////////////////////////////

// Deviation from reversability of the various integrators for the Kepler
// oscillator.

// Computes a trajectory together with its time reversal
// simultaneously for a selected list of integration methods.
// Output is the 'reversibility error' i.e. the difference between
// initial position and final position ('position' means location in phase
// space, i.e. R2(x_,v_)). Both values coincide if reversability is perfect.

// Result: The leafrog methods (i.e. LF, ALF, DALF) and the Stoermer-Verlet
// method (i.e. DMI) are reversible up to numerical noise.

void app8(void) //
// simultaneously for a selected list of integration methods
{
   Message::subProgram="app8";
   Word loc("app8");

   RecordHandler rch("app8.ini");
      // This says that we expect data input to be available in file
      // app8a.ini in a diretory which is determined by entries in
      // cpmsystemdependencies.h and cpmconfig.ini.
   Word sec="data"; // names 'rch' and 'sec' are to be employed
      // for the macro cpmrh to work
   R eps=0.25, fac=1.2, nPerRev=16,  E=0.;
  // Word imageName;
   R tWaitFinal;
   Z nTot=1000;
   V<Z> methods;
// getting quantities from section 'data' of file app1.ini
   cpmrh(E); // rh: read handler
   cpmrh(eps);
   cpmrh(nPerRev);
   //cpmrh(imageName);
   cpmrh(nTot);
   cpmrh(tWaitFinal);
   cpmrh(methods);
   S<Z> met(methods); // eliminates doublicates
   Vp<KepOsc> vk=generator(met,eps,nPerRev,E);
   Z md=met.dim();

   R_Vector dt(md);
   R_Vector xi(md);
   R_Vector vi(md);
   R_Vector xm(md);
   V<Word> metName(md);
   for (Z j=1;j<=md;++j){
    xi[j]=vk(j).x();
    vi[j]=vk(j).v();
   }

   for (Z i=1;i<=nTot;++i){
      for (Z j=1;j<=md;++j){
         vk(j).step_();
      }
   }

   for (Z j=1;j<=md;++j){
      dt[j]=vk(j).dt();
      xm[j]=vk(j).x();
      metName[j]=vk(j).met();
   }

   for (Z j=1;j<=md;++j) vk(j).rev_();

   for (Z i=1;i<=nTot;++i){
      for (Z j=1;j<=md;++j){
         vk(j).step_();
      }
   }

   R_Vector xf(md);
   R_Vector vf(md);
   R_Vector tf(md);

   for (Z j=1;j<=md;++j){
      xf[j]=vk(j).x();
      vf[j]=vk(j).v();
      tf[j]=vk(j).t();
   }

   R_Vector xL(md);
   for (Z j=1;j<=md;++j) xL[j]=1_R*j;

   R_Vector dxfi=(xf-xi).log10();
   R_Vector dxvfi(md);
   for (Z j=1;j<=md;++j){
      dxvfi[j]=cpmlog10(cpmsqrt(cpmsqr(xf[j]-xi[j])+cpmsqr(vf[j]-vi[j])));
   }
   cpmcerr<<"dt = "<<dt<<endl;
   cpmcerr<<"xi = "<<xi<<endl;
   cpmcerr<<"xm = "<<xm<<endl;
   cpmcerr<<"xf = "<<xf<<endl;
   cpmcerr<<"xf - xi = "<<xf-xi<<endl;
   cpmcerr<<"tf = "<<tf<<endl;
   cpmcerr<<"metName = "<<metName<<endl;
   Graph gr;
   gr.setAutoScale(true);
   Word title("Deviation from reversibility: log10(xf-xi)");
   gr.mark(xL,dxvfi,BLUE,0);
   gr.setNoGrid(true);
   gr.mark(xL,dxvfi,RED,2_Z,title,0.05_R);
   gr.vis();
//   Word img=gr.wrt(imageName); // writes the final frame to a ppm image file
   gr.wrt("image"); // writes the final frame to a ppm image file
      // if imageName is not the void word.
  // CpmSystem::ppm2png(img);
   cpmwait(tWaitFinal,2);
}


void app9(void) // A few oscillations of the Kepler oscillator.
// The red point shows the position of the oscillating particle
// along the x-axis and the total energy along the y-axis. The blue point marks
// the actual level of potential energy along the y-axis (the x-value is the
// same as that of the red point).
{
   Message::subProgram="app9";
   Word loc("app9");
   Color cb=darkSky; // background
   RecordHandler rch("app9.ini");
   Word sec="data";
   R eps=0.25, fac=1.2, nPerRev=16, nRev=2, acc=1e-12, E=0.;
   R frameRate=16.,  tWaitFinal1=1.;
   Word imageName;
   Word movieName;
   B fileOnShow;
   cpmrh(eps);
   cpmrh(E);
   cpmrh(nPerRev);
   cpmrh(nRev);
   cpmrh(acc);
   cpmrh(tWaitFinal1);
   cpmrh(frameRate);
   cpmrhf(imageName);
   cpmrhf(movieName);
   cpmrhf(fileOnShow);
   Z nTot=cpmrnd(nRev*nPerRev);
   R nTotInv=1_R/nTot;
   KepExc ko(eps,nPerRev,E);
   X3<R,Iv,Iv> os=ko.orbSize();
   R_Vector od=ko.orbDes();
   R a=od[1];
   R b=cpmsqrt(a); // since p = 1 ==> a = b*b
   R dt=ko.dt();
   Graph gr;
   R eT=ko.eTot();
   gr.setXY(os.c2(),Iv(-0.5_R,eT),1.1_R);
   gr.setNewGrid(1);
   R xL=os.c2().inf();
   R xU=os.c2().sup();
   R x0=ko.x();
   R_Func fpot{KepOsc::pot};
   bool makeMovie=movieName.dim()>0;

   for (Z i=1;i<=nTot;++i){
      gr.clearAndGrid();
      gr.draw(R2(xL,eT),R2(xU,eT),LIGHTGRAY);
      gr.mark(fpot,400_Z,LIGHTGRAY);
      gr.addText(Word("num. eccentricity")&"="&cpm(eps));
         // there is a version that uses classes LineArt and Camera to write
         // a greek letter epsilon here. But this destroyes the grid and the
         // grid-related legend. This version is conserved as the file
         // tut2_20230827b.cpp.
      R x=ko.x();
      R2 xPos=R2(x,ko.eTot());
      gr.mark(xPos,16_Z,RED);
      R2 xPot=R2(x,ko.ePot());
      gr.mark(xPot,16_Z,BLUE);
      gr.vis(fileOnShow||makeMovie);
      ko.step_();
      cpmprogress("loop",i*nTotInv);
   }

   if (makeMovie){ // creating video
      path dirMem = current_path();
      current_path(-Message::getOutDir2());
      path dirAct = current_path();
      cpmcerr<<"creating movie at "<<dirAct<<endl;
      ppm2mp4(movieName,frameRate);
      current_path(dirMem);
   }
}

void app10(void) // A few oscillations of the Kepler oscillator.
// The red point shows the position of the oscillating particle
// along the x-axis and the total energy along the y-axis. The blue point marks
// the actual level of potential energy along the y-axis (the x-value is the
// same as that of the red point).
{
   Message::subProgram="app10";
   Word loc("app10");
   Color cb=darkSky; // background
   RecordHandler rch("app10.ini");
   Word sec="data";
   R epsStart=0.1, epsStep=0.1, epsEnd=0.7, fac=1.2, nPerRev=16;
   R nRev=2, acc=1e-12, E=0., frameRate=16.,  tWaitFinal1=1.;
   Word imageName;
   Word movieName;
   B fileOnShow;
   cpmrh(epsStart);
   cpmrh(epsStep);
   cpmrh(epsEnd);
   cpmrh(E);
   cpmrh(nPerRev);
   cpmrh(nRev);
   cpmrh(acc);
   cpmrh(tWaitFinal1);
   cpmrh(frameRate);
   cpmrhf(imageName);
   cpmrhf(movieName);
   cpmrhf(fileOnShow);
   Z nTot=cpmrnd(nRev*nPerRev);
   R nTotInv=1_R/nTot;
   R_Vector epsList(epsStart,epsStep,epsEnd);
   Z nEps=epsList.dim();
   V<KepExc> koList(nEps);
   R_Vector xL(nEps);
   R_Vector xU(nEps);
   R_Vector eT(nEps);
   Graph gr;
   KepExc kof(epsEnd,nPerRev,E);
   R dtf=kof.dt();
   for (Z i=1;i<=nEps;++i){
      KepExc koi(epsList[i],nPerRev,E);
      koi.set_dt_(dtf);
      koList[i]=koi;
      X3<R,Iv,Iv> os=koi.orbSize();
      R_Vector od=koi.orbDes();
      R a=od[1];
      R b=cpmsqrt(a);
      eT[i]=koi.eTot();
      xL[i]=os.c2().inf();
      xU[i]=os.c2().sup();
      if (i==nEps){
         gr.setXY(os.c2(),Iv(-0.5_R,eT[i]),1.1_R);
         gr.adjXY();
         gr.clearAndGrid(10,10);
      }
   }
   R_Func fpot{KepOsc::pot};
   bool makeMovie=movieName.dim()>0;

   for (Z i=1;i<=nTot;++i){
      gr.clearAndGrid(10,10);
      gr.mark(fpot,400_Z,LIGHTGRAY);
      Word leg="numerical eccentricity: ";
      leg&=cpm(epsStart);
      leg&="(";
      leg&=cpm(epsStep);
      leg&=")";
      leg&=cpm(epsEnd);
      gr.addText(leg,0.55,0.1);
         // there is a version that uses classes LineArt and Camera to write
         // a greek letter epsilon here. But this destroyes the grid and the
         // grid-related legend. This version is conserved as the file
         // tut2_20230827b.cpp.
      for (Z j=1;j<=nEps;++j){
         gr.draw(R2(xL[j],eT[j]),R2(xU[j],eT[j]),LIGHTGRAY);
         R x=koList[j].x();
         R2 xPos=R2(x,koList[j].eTot());
         gr.fillDisk(xPos,6_R,RED);
         koList[j].step_();
      }
      gr.vis(fileOnShow||makeMovie);
      cpmprogress("loop",i*nTotInv);
   }

   if (makeMovie){ // creating video
      path dirMem = current_path();
      current_path(-Message::getOutDir2());
      path dirAct = current_path();
      cpmcerr<<"creating movie at "<<dirAct<<endl;
      ppm2mp4(movieName,frameRate);
      current_path(dirMem);
   }
}

void app11(void)
// Showing the relation between the elliptical Kepler orbit (along which the
// 'planet' --- symbolizes as a blue point --- runs) and the linear orbit of the
// Kepler oscillator (along which the 'oscillator particle' --- symbolizes as
// a red point --- runs). Notice that the motion of the blue point is determined
// (up to a factor which is determined by the numerical value of the period)
// by Kepler's second law of planetary motion and that the motion of the
// red point is coupled to that of the blue point as by gears.
{
   Message::subProgram="app11";
   Word loc("app11");
   Color cb=darkSky; // background
   RecordHandler rch("app11.ini");
   Word sec="data";
   R eps=0.25, frameRate=16.,  tWaitFinal1=1.;
   Z nPerRev=256;
   Word movieName;
   B fileOnShow;

   cpmrh(eps);
   cpmrh(nPerRev);
   cpmrh(tWaitFinal1);
   cpmrh(frameRate);
   cpmrhf(movieName);
   cpmrhf(fileOnShow);
   bool makeMovie=movieName.dim()>0;

   R a=1_R;
   V<R2> orbK=KepOsc::orb(a,eps,nPerRev);
   R ae=a*eps;
   R rMin=a*(1_R-eps);
   R2 posSun(ae,0_R);
   Graph gr;
   gr.setXY(Iv(-a,a),Iv(-a,a),1.1_R);
   gr.setAspRat(1_R);
   R nPerRevInv=1_R/nPerRev;
   for (Z i=1;i<=nPerRev;++i){
      R2 pli=orbK[i];
      R xPl=pli.x1-ae;
      R yPl=pli.x2;
      R2 pl(xPl,yPl); // now the sun is the origin
      R phi=C(xPl,yPl).arg();
      R r=C(xPl,yPl).abs();
      R2 plOsc(ae-r,0_R);
      VecGraph cir=VecGraph::circle(C(posSun[1]),r,200_Z);
      gr.clearAndGrid(4_Z,4_Z);
      gr.draw(pli,posSun,Color(BLACK));
      gr.setThcLin_(true); // the two orbits are drawn as thick lines
      gr.draw(R2(-a,0_R),R2(ae-rMin,0_R),Color(RED));
      gr.mark(orbK,V<Color>{BLUE});
      gr.setThcLin_(false);
      cir.mark(gr);
      gr.fillDisk(posSun,8_R,Color(YELLOW));
      gr.fillDisk(plOsc,6_R,Color(RED));
      gr.fillDisk(pli,6_R,Color(BLUE));
      gr.addText(Word("numerical eccentricity of Kepler orbit")&" = "&cpm(eps));
      gr.vis(fileOnShow||makeMovie);
      cpmprogress("loop",i*nPerRevInv);
   }
   cpmwait(tWaitFinal1,2);
   if (makeMovie){ // creating video
      path dirMem = current_path();
      current_path(-Message::getOutDir2());
      path dirAct = current_path();
      cpmcerr<<"creating movie at "<<dirAct<<endl;
      ppm2mp4(movieName,frameRate);
      current_path(dirMem);
   }
}

} // namespace CpmKepler

void info()
{
   cout<<"Studying the Kepler oscillator"<<endl
   <<" This function takes 0 to 1 integer arguments"<<endl
   <<" values 1 2 3 4 5 6 7 8 9 10 11 "<<endl
   <<" trigger one of the five basic modes, each other number"<<endl
   <<" triggers execution of these modes in succession."<<endl
   <<" No argument is equivalent to an argument of value 1.";
}

Z CpmApplication::main_(void)
// Selects the task from the command line, if there is a numerical
// argument following the function name. Further, the program tries to
// find and read tut2.ini. Here the filename is the name of the executable
// with all capital letters removed (so tut2MS.exe would als search for
// tut2.ini).
// Places to search for this are determined by cpmsystemdependencies.h
// and cpmconfig.ini, where the latter enjoys priority.
// If no tut2.ini can be found, the coded default selection takes
// place.
// The functionality of the tasks app1(), app2(), ..., app8() is controlled
// through files ap1.ini, app2.ini,..., app8.ini for which search places are
// determined by the same logic as those of tut2.ini.
{
   using namespace CpmKepler;
   title_="tut2";
   Message::program="tut2";
      // No longer determining the name of the ini-file
   Z mL=1;
   Z sel=1; // default selection
      // here test on matrix product
   cpmverbose=1;
   if (args_.valInd(2)){ // arrays allow testing the validity
         // of indexes without running into exceptions
      Word w2=args_[2];
      if (w2=="?"){
         info();
         return 0;
      }
      else sel=w2.toZ(); // conversion to Z (integer)
   }
   Word loc("CpmApplication::main_(), sel="&cpm(sel));
   CPM_MA

   if (sel==1) app1();
   else if (sel==2) app2();
   else if (sel==3) app3();
   else if (sel==4) app4();
   else if (sel==5) app5();
   else if (sel==6) app6(); //
   else if (sel==7) app7();
   else if (sel==8) app8();
   else if (sel==9) app9();
   else if (sel==10) app10();
   else if (sel==11) app11();
   else { // do everything sequentially
      app1();
      app2();
      app3();
      app4();
      app5();
      app6();
      app7();
      app8();
      app9();
      app10();
      app11();
   }
   cpmmessage("main_() returns the selection parameter sel = "&cpm(sel)&
   "\n so, a value different from 0 needs not to signal a problem");
   CPM_MZ
   return sel;
}

// end of tut2.cpp