//? eqm2.cpp
//? C+- by Ulrich Mutze. Status of work 2021-10-26.
//? Copyright (c) 2021 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.


#include <cpmc.h>
#include <cpmfa.h>
#include <cpmm.h>
#include <cpmgraph.h>
#include <cpmframes.h>
#include <cpmimagingtools.h>
#include <cpmlinalg.h>
#include <cpmcvector.h>
#include <cpmapplication.h>
#include <cpmalgorithms.h>
#include <cpmrecordhandler.h>
#include <cpmfile.h>
#include <cpmrfunc.h>
#include <cpmobservable.h>
#include <cpmoperator.h>
#include <compare>
#include <functional>
#include <set>
#include <algorithm>
#include <vector>
#include <valarray>
#include <cmath>

using namespace CpmRoot;
using namespace CpmRootX;
using namespace CpmFunctions;
using namespace CpmSystem;
using namespace CpmArrays;
using namespace CpmGraphics;
using namespace CpmImaging;
using namespace CpmApplication;
using namespace CpmAlgorithms;
using namespace CpmGeo;
using namespace CpmQuantumMechanics;
using namespace std;

/* Expansive Quantum Mechanics.
We consider a particularly simple and transparent quantum mechanical model
system in which an initially insolated subsystem A and a subsystem B
figurating as the evironment of A can be identified.
System A consists initially of a free particle in uniform motion.
System B consists initially of a large number of localized potential wells
each occupied by a particle in its ground state. All particles interact
with the potential wells and with all other particles via a universal pair
interaction. Physical space is idealized as a finite chain of equi-distant
points.
*/

namespace eqm2{ // expansive quantum mechanics
//#undef CPM_VEC
//#define CPM_VEC
   // the program test4 gave execution time
   // 69.83 using valarray and
   // 70.30 using vector
   // So there is no difference between the methods so far. Lets see whether
   // this is still valid if more than one particles are in the game.
   // The actual version of the program fisihes irregularly with a
   // std::bad_alloc exception for both settings of this macro at
   // approximately at system time 0.12. There seems to be a memory
   // leak. How can this occur with the main classes having constructors
   // and (default) destructors ? Befor this is not resolved the
   // present version of expansive quantum mechanics is not useful.
   // No leaking as can be seen by running the 'Systemueberwachung'
   // parallel to eqm2. One sees memory usage mostly oscillating between
   // 189 MB and 221 MB, very rarely exceeding 300 MB.
   // Is a time limit built in ??? No !!!
   // The problematic point is the introduction of the 4th particle.
   // In our case a single state would require 32*32 MB = 1GB.
   // Since we need two of those simultaneously we need 2 GB of RAM
   // Plausible that these can't be allocated.
   // One should try to work with mutating methods in place where I
   // used generating ones. (already implemented !!!)
   // 2023-02-20 not clear to what an degree this is still relevant.
   // 2023-03-16 Comparison of eqm2.cpp (multiindexes using function c()
   // for elegance) and eqm2Mem6.cpp (old treatment of multiindexes):
   // new elegant version 1152.31 s , old method 1159.65 s execution time
   // for the same input data:  cpmSel/eqm2/r20230316c new,
   // cpmSel/eqm2/r20230316b old. Both runs with optimization O2.
   // Quintessence: No reasopn for conserving the old version.

//#define CPM_VEC

#if defined(CPM_VEC)
   using VN = vector<N>;
   using VC = vector<C>;
   using VR = vector<R>;
#else
   using VN = valarray<N>;
   using VC = valarray<C>;
   using VR = valarray<R>;
#endif

R abs(VC const& x){
   R sum=0_R;
   N n{x.size()};
   for (N i=0; i<n; ++i) sum+=x[i].abs();
   return sqrt(sum);
}

// defining the arithmetic operations for std::vector which for
// std::valarray are defined in the standard library

#if defined(CPM_VEC)

vector<C> operator+(vector<C> const& v1, vector<C> const& v2)
{
   N n{v1.size()};
   vector<C> res(n);
   for (N i=0; i<n; ++i) res[i]=v1[i]+v2[i];
   return res;
}

vector<C> operator-(vector<C> const& v1, vector<C> const& v2)
{
   N n{v1.size()};
   vector<C> res(n);
   for (N i=0; i<n; ++i) res[i]=v1[i]-v2[i];
   return res;
}

vector<C> operator*(vector<C> const& v1, vector<C> const& v2)
{
   N n{v1.size()};
   vector<C> res(n);
   for (N i=0; i<n; ++i) res[i]=v1[i]*v2[i];
   return res;
}

vector<C> operator*(vector<C> const& v1, C const& z)
{
   N n{v1.size()};
   vector<C> res(n);
   for (N i=0; i<n; ++i) res[i]=v1[i]*z;
   return res;
}

vector<C>& operator+=(vector<C>& v1, vector<C> const& v2)
{
   N n{v1.size()};
   for (N i=0; i<n; ++i) v1[i]+=v2[i];
   return v1;
}

vector<C>& operator-=(vector<C>& v1, vector<C> const& v2)
{
   N n{v1.size()};
   for (N i=0; i<n; ++i) v1[i]-=v2[i];
   return v1;
}

vector<C>& operator*=(vector<C>& v, C z)
{
   N n{v.size()};
   for (N i=0; i<n; ++i) v[i]*=z;
   return v;
}

vector<C>& operator*=(vector<C>& v, R t)
{
   N n{v.size()};
   for (N i=0; i<n; ++i) v[i]*=t;
   return v;
}

#endif

void show(VN const& vn){
   cout<<"show(VN) started"<<endl;
   N sz=vn.size();
   cout<<" size = "<<sz<<endl;
   for (N i=0;i<sz;++i){
      cout<<"vn["<<i<<"] = "<<vn[i]<<endl;
   }
   cout<<"show(VN) done"<<endl;
}

//. scalar product
// see C& Wf::operator()(VN const& ind) for the main usage
N sp(VN const& vn1, VN const& vn2)
{
   N sum{0};
   N n{vn1.size()};
   for (N i=0;i<n;++i) sum+=vn1[i]*vn2[i];
   return sum;
}


// to multi-index
// i is the running index for which the multi-index representation has to
// be found. ind gives the ranges of sub-indices which make up the multi-
// index. Let be ind=(n0,n1,...,np-1). and i=(i_0,i_1,...,i_(np-1))
// Then the range of i_0 is [0,...,n0-1],
// the range of i_1 is [0,...,n1-1],
//   ...
// the range of i_(np-1) is [0,...,np-1]

VN toMultInd(N i, VN const& ind)
{
   N np=ind.size();
   auto q1=CpmRoot::div(i,ind[0]);
   if (np==1) return VN{q1.first};
   auto q2=CpmRoot::div(q1.second,ind[1]);
   if (np==2) return VN{q1.first,q2.first};
   auto q3=CpmRoot::div(q2.second,ind[2]);
   if (np==3) return VN{q1.first,q2.first,q3.first};
   auto q4=CpmRoot::div(q3.second,ind[3]);
   if (np==4) return VN{q1.first,q2.first,q3.first,q4.first};
   cpmassert(np<=4,"we won't consider more than 4 particles");
   return ind; // never reached
}

// dividing the interval [xL,xU] into m equal pieces gives an interval of
// length dIni(m,xL,xU)
R dIni(N m, R xL, R xU){
   return (xU-xL)/static_cast<R>(m);
}

// returns m^np in Fortran-notation
N szIni(N m, N np){
   return static_cast<N>(pow(m,np));
}

// function for initializing Bio::ind_
// This allows only to use biotopes which have the same size for
// all particles in it.
VN indIni(N m, N np)
{
   N i{1};
   N k;
#if defined(CPM_VEC)
   VN prel(np,1);
#else // then valarry
   VN prel(1,np);
#endif
   // preliminary result initialized with np numbers one
   cpmassert(np==prel.size(),"indIni(N m, N np)");
   while (i<np){
      prel[i]=prel[i-1]*m; i++;
   }
      // now prel holds the the np powers of m in the order of
      // increasing exponents
   VN ind(np);
   for (i=0,k=np-1; i<np; ++i,--k) ind[i]=prel[k];
      // now ind holds the np powers of m in the order of
      // decreasing exponents
   return ind;
}

VR xIni(N m, R xL, R d)
{
   VR x(m);
   R x0=xL+d*0.5;
   for (N i=0;i<m;++i) x[i]=d*i+x0;
   return x;
}

function<R(R)> fp{[](R x){ // defining a function in terms of an anonymous lambda
   // expression.
      R rho{0.5};
      R x2{x*x};
      R rho2{rho*rho};
      return x2 > rho2 ? R{0.} : x2-rho2;
   }
}; // notice that the definition of fp ends in ';' unlike the definition
   // of a functon.

N aN{1};
N sizeN=sizeof(aN);

R aR{1.};
N sizeR=sizeof(aR);

N sizeC=sizeR*2;


class Bio{ // biotope, only geometric data

   N m_{0}; // number of positions available to any of the np_
      // particles
   N np_{0};
      // number of particles. In the present program limited to 4
   R xL_{0.};
   R xU_{0.};
      // the m_ positions are the centers of the intervals obtained
      // by dividing the interval [xL,xU] into m_ congruent parts.
   R d_{0.}; // gap between neigboring particle positions

   N sz_{0};// m_^np_ in Fortran notation
      // This will be the dimension of Wf::x_

   VN ind_;
      // a list of np_ natural numbers which allow an economic
      // way of transforming multi-indexes into a single large index

   VR x_; // list of the positions available to any of the np_ partlcles.
      // thus a list of m_ real numbers

   Iv ivMax_;
      // we don't consider values xL_, xU_ for which the interval [xL_,xU_]
      // is not contained in in the interval ivMax_

public:
   Bio(N m, N np, R xL, R xU, Iv ivMax=Iv(-1000_R, 1000_R)):
      m_{m}, np_{np}, xL_{xL}, xU_{xU},
      d_{dIni(m_,xL_,xU_)}, sz_{szIni(m_,np_)}, ind_{indIni(m_,np_)},
      x_{xIni(m_,xL_,d_)}, ivMax_{ivMax} {}

   Bio(){} // required to satisfy the concept semiregular
   N m()const{ return m_;}
   N& m(){ return m_;}
   N np()const{ return np_;}
   R xL()const{ return xL_;}
   R xU()const{ return xU_;}
   R& xL(){ return xL_;}
   R& xU(){ return xU_;}
   R d()const{ return d_;}
   N sz()const{ return sz_;}
   VN ind()const{ return ind_;}
   VR x()const{ return x_;}
   Iv ivMax()const{ return ivMax_;}

   Bio shift(R sh)const;

   void toSingle_(){ np_=1; sz_=m_; ind_=indIni(m_,np_);}

   bool match(Bio const& b)const;

   bool contains(Bio const& b)const;

   // If the return value T2<N>(nL,nU) differs from (-1,-1) the biotope
   // *this contains the biotope b and has nL additional points
   // at the lower end of b and nU additional points at the upper end
   // end of b. Of course nL and/or nU may be equal to 0.
   T2<N> excess(Bio const& b)const;

   N memUse()const{ return (m_+3)*sizeR+(np_+3)*sizeN;}
      // plain memory used in bytes for an instance of Bio (not Wf !).

   // most convenient way of writint to cpmcerr and cout (only if _CONSOLE is
   // set) and statusbar
   void write(std::ostringstream& ost)const{
      ost<<"m_="<<m_<<", np_="<<np_<<", xL_="<<xL_<<", xU_="<<xU_<<endl;
      cpmmessage(ost,1,true);
   }
};

   Bio Bio::shift(R sh)const{
      Bio res(*this);
      res.xL_=xL_+sh;
      res.xU_=xU_+sh;
      res.x_=xIni(res.m_,res.xL_,res.d_);
      return res;
   }

   bool Bio::match(Bio const& b)const{
      R tiny=1.e-12_R;
      bool bd=cpmabs(d_-b.d_)/d_<tiny;
      R fL=cpmabs(xL_-b.xL_)/d_;
      R fU=cpmabs(xU_-b.xU_)/d_;
      bool bLf=cpmabs(fL-cpmfloor(fL))<tiny;
      bool bLc=cpmabs(fL-cpmceil(fL))<tiny;
      bool bUf=cpmabs(fU-cpmfloor(fU))<tiny;
      bool bUc=cpmabs(fU-cpmceil(fU))<tiny;
      return bd&&(bLf||bLc)&&(bUf||bUc);
   }

   bool Bio::contains(Bio const& b)const{
      bool bm=match(b);
      bool bL= (xL_ <= b.xL_);
      bool bU= (xU_ >= b.xU_);
      return bm&&bL&&bU;
   }

   T2<N> Bio::excess(Bio const& b)const{
      if (!contains(b)) return T2<N>(-1_N,-1_N);
      return T2<N>(cpmtoz((b.xL_-xL_)/d_),cpmtoz((xU_-b.xU_)/d_));
   }

struct Dyn{ //dynamics: t,mass,potentials
   R t_{0.}; // introduced 2022-06-25
   R mass_{0.};
   mutable R tiny_{0.}; // aha! here it is. Following code never chances tiny !
      // comments say that this was different before.
   mutable B cyclic_; // initialized as false (important)

   R_Func pot_{};
   R_Func potPair_{};

   Dyn(){}

   Dyn(R mass, R tiny, R_Func const& pot, R_Func const& potPair):
      mass_{mass},tiny_{tiny},
      pot_(pot),potPair_(potPair){}

   R mass(){ return mass_;}
   R t()const{ return t_;}
   void set_t(R const& t){ t_=t;}

   N memUse()const{ return sizeR*2+1+sizeN*2;}
};

//////////////////// class Wf /////////////////////////////////////////////
// wave function, actually rep_ is the wave function

class Wf{
// independent variables
   Bio bio_;
   Dyn dyn_; // contains the time parameter, accessibla via Dyn::t() and
      // Dyn::set_t(R const&)
   VC rep_;
   mutable B readyL_;
   mutable B readyU_;
   mutable B alarmL_;
   mutable B alarmU_;
   mutable B bioFix_;
      // true if biotope has been set fixed (initialized as false)
   mutable R dt_{-1.};
   mutable R normH_{0_R};
   N smoothSteps{0_N};

public:
   Wf(){} // required to satisfy the concept semiregular
   Wf(Bio const& bio, Dyn const& dyn):
      bio_(bio), dyn_(dyn),rep_(bio_.sz()){}
#if defined(CPM_VEC)
   Wf(Bio const& bio, Dyn const& dyn, C const& t):
      bio_(bio), dyn_(dyn),rep_(bio_.sz(),t){}
#else
   Wf(Bio const& bio, Dyn const& dyn, C const& t):
      bio_(bio), dyn_(dyn),rep_(t,bio_.sz()){}
#endif

   Wf(Bio const& bio, Dyn const& dyn, VC const& rep):
      bio_(bio), dyn_(dyn),rep_(rep){
      cpmassert(rep_.size()==bio_.sz(),
      "dimension mismatch in constructor Wf(Bio,Dyn,vector<C>)")}

   Wf(Bio const& bio, Dyn const& dyn, F<R,C> const& wfunc):
      bio_(bio), dyn_(dyn),rep_(bio_.sz()){
      cpmassert(rep_.size()==bio_.sz(),
      "dimension mismatch in constructor Wf(Bio,Dyn,F<R,C>)");
      cpmassert(bio_.np()==1,"only one particle allowed");
         for (N i=0;i<bio_.sz();++i) rep_[i]=wfunc((bio_.x())[i]);
      }

   // Return a modified version of this. Modification consists of
   // replacing the wave function.
   Wf modify(VC const& rep)const{
      Wf res(*this);
      for (N i=0;i<bio_.sz();++i) res[i]=rep[i];
      return res;
   }

   // simplify
   void simplify_(VC const& rep){
      bio_.toSingle_();
      //cpmassert(rep.size()==m(),"dimension mismatch");
      cout<<"rep.size()="<<rep.size()<<", m()="<<m()<<endl;
      rep_=rep;
   }

   // applying the (1/4, 1/2, 1/4) smoothing kernel to a copy of *this
   // and return this modified copy. If n is larger than 0, the procedure gets
   // repeated n times.
   Wf smoothS()const;
      // S for single, smooth1 yet defined with different meaning

   Wf smooth(Z steps)const{
      if (steps==0) return *this;
      Wf res{*this};
      while (steps--) res=res.smoothS();
      return res;
   }

   void setCyclic(bool b)const{dyn_.cyclic_=b;}

 // call of Hamiltonian is necessary to update the alarm variables
   bool alrL()const;
   bool alrU()const;
   bool free()const{ return !alrL()&&!alrU();}
   bool freeL()const{ return !alrL();}
   bool freeU()const{ return !alrU();}

   N memUse()const{
      return dyn_.memUse()+bio_.memUse()+sizeC*sz()+sizeR+2;}
      // plain memory used in bytes

// How to access the primary data rep_ by means of a single index i.
   C const& operator[](N i)const{
      return rep_[i];
   }

   C& operator[](N i){
      return rep_[i];
   }

// nicer form for function cp (stands for 'component')
// unfortunately operator[] does not allow multiple arguments as functions do

   C const& c(N i)const{
     return rep_[i];
   }

   C const& c(N i, N j)const{
     return (*this)(VN({i,j}));
   }

   C const& c(N i, N j, N k)const{
     return (*this)(VN({i,j,k}));
   }

   C const& c(N i, N j, N k, N l)const{
     return (*this)(VN({i,j,k,l}));
   }

   C& c(N i){
     return rep_[i];
   }

   C& c(N i, N j){
     return (*this)(VN({i,j}));
   }

   C& c(N i, N j, N k){
     return (*this)(VN({i,j,k}));
   }

   C& c(N i, N j, N k, N l){
     return (*this)(VN({i,j,k,l}));
   }


// How to access the primary data rep_ by means of a multiindex ind.
   C const& operator()(VN const& ind)const{
      N s=sp(bio_.ind(),ind);
      return rep_[s];
   }

   C& operator()(VN const& ind){
      N s=sp(bio_.ind(),ind);
      return rep_[s];
   }

 // How multi-indexing and single indexing work together:
 // For any i such that 0<=i<sz_ define
 // VN vn = (*this)(i);
 // T val1 = (*this)[i];
 // T val2 = (*this)(vn);
 // Then we always have: val1 == val2.
 // This shows that vn is the multi-index which refers to the same part
 // of the data rep_ as the single index i
    VN operator()(N i)const{
      return toMultInd(i, bio_.ind());
   }

   bool match(Wf const& a)const{
      return bio_.m()==a.bio_.m() && bio_.np()==a.bio_.np() &&
      bio_.sz()==a.bio_.sz() && bio_.xL()==a.bio_.xL() &&
      bio_.xU()==a.bio_.xU();
   }

   VC const& rep()const{ return rep_;}
   Bio bio()const{ return bio_;}
   Dyn dyn()const{ return dyn_;}
   R t()const{ return dyn_.t();}
   void set_t(R const& t){dyn_.set_t(t);}
   N m()const{ return bio_.m();}
   N np()const{ return bio_.np();}
   R xL()const{ return bio_.xL();}
   R xU()const{ return bio_.xU();}
   Iv iv()const{ return Iv(bio_.xL(),bio_.xU());}
   R d()const{ return bio_.d();}
   N sz()const{ return bio_.sz();}
   VN ind()const{ return bio_.ind();}
   VR x()const{ return bio_.x();}
   R mass()const{ return dyn_.mass_;}
   R_Func pot()const{ return dyn_.pot_;}
   R_Func potPair()const{ return dyn_.potPair_;}
   bool alarmL()const{ return alrL();}
   bool alarmU()const{ return alrU();}
   R dt()const{ return dt_;}
   void setSmoothSteps_(N steps){ smoothSteps=steps;}
   void setAlarmL(bool b){ alarmL_=B{b};}
   void setAlarmU(bool b){ alarmU_=B{b};}
   void setPot_(R_Func const& f){ dyn_.pot_=f;}
   void setWaveFunc_(VC const& vc){ rep_=vc;}
   void cyclic(bool b){dyn_.cyclic_=b;}
   void cyclic_(bool b){dyn_.cyclic_=b;}
      // since the quantity that is changed is declared mutable
      // both names ('cyclic' and 'cyclic_') are adequate

   bool isCyclic()const{ return dyn_.cyclic_;}

// show
// Text representation of *this. Restricted to 4 particles.
   void show()const;

// Getting a one-particle density matrix over the biotope by 'tracing out'
// all particles except to one indexed p.
   C_Matrix traceOut2(N p=0)const;

// Getting a real-valued 'wave function' over the biotope by 'tracing out'
// all particles except the one indexed p. It is instructive to consider
// a free two-particle wave function in the case that the two particles
// each have a peak-like wave function with different peak positions for
// the two particles. Than traceOut(0) will show one of the two peaks and
// traceOut(1) will show the other one. None of the two traceOut functions
// will show two peaks! This irritaded me for some uncomfortable hours.
// Attempts to define a similar construction which gives a more realistic
// complex-valued wave function did not work so far.
   R_Vector traceToVec(N p)const;

// Returns the wave function as a C_Vector. Implemented only for the case
// that there is omnly one particle.
   C_Vector toVec()const;

//: mark
// Showing a real-valued 'wave function' over the biotope by 'tracing out'
// all particles except to one indexed p. Notice that valid indexes start
// with 0.
   void mark(Graph& gr, Iv ivx=Iv(0_R,0_R))const;

// shift in the sense of translation
   Wf shift(Z s)const;

//: norm
   R norm()const{
      R sum{0.};
      for (N i=0;i<sz();++i) sum+=rep_[i].abs2();
      return cpmsqrt(sum);
   }

   R abs2()const{
      R sum{0.};
      for (N i=0;i<sz();++i) sum+=rep_[i].abs2();
      return sum;
   }

   Wf normalize()const;

   //: normalize
   // Changes *this into a unit vector and returns the value of the
   // norm prior to normalization.
   R nor_();

// scalar product
   C operator|(Wf const& a)const{
      cpmassert(sz()==a.sz(),
         "scalar product needs equal number of components");
      C res;
      for (N i=0;i<=sz();++i) res+=(rep_[i]|a[i]);
      return res;
   }

   Wf operator+(Wf const& a)const{
      cpmassert(match(a),"trying to add non-matching Wfs");
      return Wf(bio_,dyn_,rep_+a.rep_);
   }

   Wf& operator+=(Wf const& a){
      rep_+=a.rep_;
      return *this;
   }

   Wf operator-(Wf const& a)const{
      cpmassert(match(a),"trying to subtract non-matching Wfs");
      return Wf(bio_,dyn_,rep_-a.rep_);
   }

   Wf& operator-=(Wf const& a){
      rep_-=a.rep_;
      return *this;
   }

   Wf operator*(Wf const& a)const{
      cpmassert(match(a),"trying to multiply non-matching Wfs");
      return Wf(bio_,dyn_,rep_*a.rep_);
   }

   Wf operator*(C const& t)const{
      Wf res(*this);
      for (N i=0;i<sz();++i){ res[i]*=t;}
      return res;
   }

   Wf operator*(R t)const{
      Wf res(*this);
#if defined(CPM_VEC) // here vector is more flexible
      res.rep_*=t;
#else
      res.rep_*=C(t);
#endif
      return res;
   }

   Wf& operator*=(R t){
      rep_*=C(t);
      return *this;
   }

   Wf& boost_(R const& v){
      cpmassert(np()==1, "we ony boost the first particle");
      for (N i=0;i<m();++i){
         C zi=C::expi(v*(bio_.x())[i]*dyn_.mass_);
         rep_[i]*=zi;
      }
      return *this;
   }

   Wf boost(R const& v)const{
      cpmassert(np()==1, "we ony boost the first particle");
      Wf res(*this);
      for (N i=0;i<m();++i){
         C zi=C::expi(v*(bio_.x())[i]*dyn_.mass_);
         res[i]*=zi;
      }
      // since |zi|==1, boost does not change the norm of res
      return res;
   }

   //: randomize
   Wf ran(Z seed=0)const;

   // shift biotope
   // Returns a modified version of *this in which only the data member
   // bio_ is changed
   Wf shiftBio(R shf)const{
      Wf res(*this);
      res.bio_=bio_.shift(shf);
      return res;
   }

   //. extend at lower end
   // We add n points at the lower end of the biotope.
   // No particle added.
   // np_ is assumed to be not larger than 4.
   Wf extL(N n)const;

   //. extend at upper end
   // We add n points at the upper end of the biotope.
   // No particle added.
   // np_ is assumed to be not larger than 4.
   Wf extU(N n)const;

   //. extend biotope at upper end and add one particle which lives
   // on the added part of the biotope
   Wf extU(Wf const& a)const;

   //: tensor product
   // a is the wave function of a one-particle state on the same
   // biotope as *this
   // and the result of the function is original *this multiplied
   // tensorially by a
   Wf tenPrd(Wf const& a)const;

   // kinetic energy of free particles in the full form which allows
   // growth of the biotope
   Wf hf()const;

   // kinetic energy of free particles in a lean form that assumes
   // reflecting walls and forbids further changes of the biotope
   Wf hl()const;

   // kinetic energy of free particles
   Wf h()const{ return bioFix_ ? hl() : hf();}

   // potential energy in external field
   Wf v()const;

   // potental energy by pair interaction
   Wf vp()const;

   // Hamilton operator
   //Wf ham()const{ return (*this).h()+(*this).v()+(*this).vp(); }
   Wf ham()const{ return h()+v()+vp(); }

   Wf dot()const{ return ham()*C(R{0.},R{-1.});}

   // operator norm of the hamiltonian ham
   void hNorm_(Z iMax=25, R tWait=0);

   // operator norm of the kinetic hamiltonian (no potentials)
   R hKinNorm()const;

   R getTimeStep(R rel, Z iMax=16, R tWait=0.){
      Z mL=1;
      static Word loc("getTimeStep(R,Z,R)");
      CPM_MA
      hNorm_(iMax,tWait);
      dt_= R{2.}*rel/normH_;
      ostringstream ost;
      ost<<"dt_ = "<<dt_<<endl;
      cpmmessage(ost,2,true);
      CPM_MZ
      return dt_;
   }
};

   void Wf::show()const{
      if (np()==1){
         for (N i=0; i<m(); ++i){
            cout<<"rep_["<<i<<"]="<<rep_[i]<<endl;
         }
      }
      if (np()==2){
         for (N i=0; i<m(); ++i){
            for (N j=0; j<m(); ++j){
              cout<<"i="<<i<<", j="<<j<<", val="<<c(i,j)<<endl;
            }
         }
      }
      if (np()==3){
         for (N i=0; i<m(); ++i){
            for (N j=0; j<m(); ++j){
               for (N k=0; k<m(); ++k){
                  cout<<"i="<<i<<", j="<<j<<", k="<<k<<", val="
                  <<c(i,j,k)<<endl;
               }
            }
         }
      }
      if (np()==4){
         for (N i=0; i<m(); ++i){
            for (N j=0; j<m(); ++j){
               for (N k=0; k<m(); ++k){
                  for (N l=0; l<m(); ++l){
                     cout<<"i="<<i<<", j="<<j<<", k="<<k<<", l="<<l<<", val="
                     <<c(i,j,k,l)<<endl;
                  }
               }
            }
         }
      }
      else cpmassert(np()<=4,"number of particles larger than 4");
   }

   C_Vector Wf::toVec()const{
      cpmassert(np()==1,"function is defined only for np==1");
      C_Vector res(m());
      for (N i=0; i<m(); ++i){
         res[i+1]=rep_[i];
      }
      return res;
    }

    C_Matrix Wf::traceOut2(N p)const{
      Z mL=1;
      Word loc("Wf::traceOut2("&cpm(p)&")const");
      CPM_MA
      cpmassert(p<np()," p out of range");
      C_Matrix y(m());
      if (np()==1){ //no influence of p
         for (N i=0; i<m(); ++i){
            for (N j=0; j<m(); ++j){
               y[i+1][j+1]=rep_[i]*~rep_[j];
            }
         }
         CPM_MZ
         return y;
      }
      if (np()==2){
         if (p==0){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     yij+=c(i,l)*~c(j,l);
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
         if (p==1){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     yij+=c(l,i)*~c(l,j);
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
      }
      if (np()==3){
         if (p==0){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        yij+=c(i,l,q)*~c(j,l,q);
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
         if (p==1){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        yij+=c(l,i,q)*~c(l,j,q);
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
         if (p==2){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        yij+=c(l,q,i)*~c(l,q,j);
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
      }
      if (np()==4){
         if (p==0){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        for (N r=0; r<m();++r){
                           yij+=c(i,l,q,r)*~c(j,l,q,r);
                        }
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
         if (p==1){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        for (N r=0; r<m();++r){
                           yij+=c(l,i,q,r)*~c(l,j,q,r);
                        }
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
         if (p==2){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        for (N r=0; r<m();++r){
                           yij+=c(l,q,i,r)*~c(l,q,j,r);
                        }
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
         if (p==3){
            for (N i=0; i<m(); ++i){
               for (N j=0; j<m(); ++j){
                  C yij;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        for (N r=0; r<m();++r){
                           yij+=c(l,q,r,i)*~c(l,q,r,j);
                        }
                     }
                  }
                  y[i+1][j+1]=yij;
               }
            }
            CPM_MZ
            return y;
         }
      }
      CPM_MZ
      return y;
   }

   R_Vector Wf::traceToVec(N p)const{
      Z mL=1;
      Word loc("Wf::traceOut1("&cpm(p)&")const");
      CPM_MA
      cpmassert(p<np()," p out of range");
      R_Vector y(m());
      if (np()==1){ //no influence of p
         for (N i=0; i<m(); ++i){
            for (N j=0; j<m(); ++j){
               y[i+1]=rep_[i].abs2();
            }
         }
         CPM_MZ
         return y;
      }
      if (np()==2){
         if (p==0){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  yi+=c(i,l).abs2();
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
         if (p==1){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  yi+=c(l,i).abs2();
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
      }
      if (np()==3){
         if (p==0){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  for (N q=0; q<m();++q){
                     yi+=c(i,l,q).abs2();
                  }
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
         if (p==1){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  for (N q=0; q<m();++q){
                     yi+=c(l,i,q).abs2();
                  }
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
         if (p==2){
            for (N i=0; i<m(); ++i){
              auto yi=0_R;
                  for (N l=0; l<m();++l){
                     for (N q=0; q<m();++q){
                        yi+=c(l,q,i).abs2();
                     }
                  }
                  y[i+1]=yi;
               }

            CPM_MZ
            return y;
         }
      }
      if (np()==4){
         if (p==0){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  for (N q=0; q<m();++q){
                     for (N r=0; r<m();++r){
                        yi+=c(i,l,q,r).abs2();
                     }
                  }
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
         if (p==1){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                 for (N q=0; q<m();++q){
                     for (N r=0; r<m();++r){
                        yi+=c(l,i,q,r).abs2();
                     }
                  }
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
         if (p==2){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  for (N q=0; q<m();++q){
                     for (N r=0; r<m();++r){
                        yi+=c(l,q,i,r).abs2();
                     }
                  }
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
         if (p==3){
            for (N i=0; i<m(); ++i){
               auto yi=0_R;
               for (N l=0; l<m();++l){
                  for (N q=0; q<m();++q){
                     for (N r=0; r<m();++r){
                        yi+=c(l,q,r,i).abs2();
                     }
                  }
               }
               y[i+1]=yi;
            }
            CPM_MZ
            return y;
         }
      }
      CPM_MZ
      return y;
   }

   void Wf::mark(Graph& gr, Iv ivx)const{
      Z mL=1;
      static Word loc("Wf::mark(Graph&,Iv)");
      CPM_MA
      R t1=cpmtime();
      ivx.abs()==0_R ? gr.setX(Iv(bio_.xL(), bio_.xU())) : gr.setX(ivx);
      R_Vector x(bio_.x());
      R frcNeg=0.2_R;
      R_Vector potL=dyn_.pot_.eval(x); // ad hoc down scaling of
         // potential for getting more space for the living part of the
         // graphics
      if (np()==1){ // Here we can show the true wave function. In all of the
         // following cases we show non-faithful modificatons of the wave
         // function which are enforced by the requirement to allow a graphical
         // representation.
         C_Vector y0=toVec();
         R_Vector y00=y0.real();
         R_Vector y01=y0.imag();
         R_Vector y02=y0.absC();
         R_Matrix y(4);
         y[1]=y00;
         y[2]=y01;
         y[3]=y02;
         y[4]=potL;
         V<Color> cl{GREEN,RED,BLACK,ORANGE};
         gr.mark(x,y,cl);
      }
      else if (np()==2){
         R_Vector y0=traceToVec(0);
         R_Vector y1=traceToVec(1);
         R_Matrix y(3);
         y[1]=y0;
         y[2]=y1;
         y[3]=potL;
         y=y.scaleNeg(frcNeg);
         V<Color> cl{BLACK,LIGHTMAGENTA,ORANGE};
         gr.mark(x,y,cl);
      }
      else if (np()==3){
         R_Vector y0=traceToVec(0);
         R_Vector y1=traceToVec(1);
         R_Vector y2=traceToVec(2);
         R_Matrix y(4);
         y[1]=y0;
         y[2]=y1;
         y[3]=y2;
         y[4]=potL;
         y=y.scaleNeg(frcNeg);
         V<Color> cl{BLACK,LIGHTMAGENTA,BLUE,ORANGE};
         gr.mark(x,y,cl);
      }
      else if (np()==4){
         R_Vector y0=traceToVec(0);
         R_Vector y1=traceToVec(1);
         R_Vector y2=traceToVec(2);
         R_Vector y3=traceToVec(3);

         R_Matrix y(5);
         y[1]=y0;
         y[2]=y1;
         y[3]=y2;
         y[4]=y3;
         y[5]=potL;
         y=y.scaleNeg(frcNeg);
         V<Color> cl{BLACK,LIGHTMAGENTA,BLUE,GREEN,ORANGE};
         gr.mark(x,y,cl);
      }
      R t2=cpmtime();
      cpmmessage(loc," done in "&cpm(t2-t1)&" s");
      CPM_MZ
   }

   N cycShift(N i, Z s, N m){
      Z is=(Z)i-s;
      while (is<0){ is+=m;}
      while (is>=m){is-=m;}
      return (N)is;
   }

   Wf Wf::shift(Z s)const{
      cpmassert(np()==1, "we only shift the first particle");
      Wf res(*this);
      for (N i=0;i<m();++i){
         N ic=cycShift(i,s,bio_.m());
         res[i]=rep_[ic];
      }
      return res;
   }

   Wf Wf::normalize()const{
      R r=0.;
      for (N i=0;i<sz();++i) r+=rep_[i].abs2();
      cpmassert(r>0.,"zero vectors can't be normalized");
      R fac=1_R/cpmsqrt(r);
      VC rep(rep_);
      for (N i=0;i<sz();++i){
         rep[i]*=fac;
      }
      return modify(rep);
   }

   R Wf::nor_(){
      R r{0.};
      for (N i=0;i<sz();++i) r+=rep_[i].abs2();
      R res=cpmsqrt(r);
      if (res!=0.){
         R fac=1_R/res;
         for (N i=0;i<sz();++i) rep_[i]*=fac;
      }
      return res;
   }

   Wf Wf::ran(Z seed)const
   {
      VC rep=rep_;
      for (N i=0;i<sz();++i) rep[i]=Root(rep_[i]).ran(seed);
      return Wf(bio_,dyn_,rep);
   }

   Wf Wf::extL(N n)const{
      Z mL=1;
      static Word loc("Wf::extL(N n)");
      R t1=cpmtime();
      CPM_MA
      Bio bioL(bio_.m()+n, bio_.np(), bio_.xL()-bio_.d()*n, bio_.xU());
      Iv ivMax=bio_.ivMax();
      R xLMax=ivMax.b();
      R xUMax=ivMax.e();
      bool bL=false;
      bool bU=false;
      if (bioL.xL()<= xLMax) {bioL.xL()=xLMax; bL=true;}
      if (bioL.xU()>= xUMax) {bioL.xU()=xLMax; bU=true;}
      N m=bioL.m();
      cpmassert(n>=3,"n should be larger than 2")
      Wf res(bioL,dyn_); // res carries the time information of *this
      if (np()==1){
         for (N i=0;i<m;++i) res[i] = i>=n ? rep_[i-n] : C{};
      }
      else if (np()==2){
         for (N i=0;i<m;++i){
            for (N j=0;j<m;++j){
               res.c(i,j) = (i>=n && j>=n) ? c(i-n,j-n) : C{};
            }
         }
      }
      else if (np()==3){
         for (N i=0;i<m;++i){
            for (N j=0;j<m;++j){
               for (N k=0;k<m;++k){
                  res.c(i,j,k) =
                  (i>=n && j>=n && k>=n) ? (c)(i-n,j-n,k-n) : C{};
               }
            }
         }
      }
      else if (np()==4){
         for (N i=0;i<m;++i){
            for (N j=0;j<m;++j){
               for (N k=0;k<m;++k){
                  for (N l=0;k<l;++l){
                     res.c(i,j,k,l) =
                        (i>=n && j>=n && k>=n && l>=n) ?
                           c(i-n,j-n,k-n,l-n) : C{};
                  }
               }
            }
         }
      }
      else cpmassert(np()<=4,"number of particles larger than 4");
      res.readyL_=true;
      res.alarmL_=false;
      res.smooth(smoothSteps);
      R t2=cpmtime();
      cpmmessage(loc," done in "&cpm(t2-t1)&" s");
      CPM_MZ
      return res;
   }

   Wf Wf::extU(N n)const{
      Z mL=1;
      static Word loc("Wf::extU(N n)");
      CPM_MA
      R t1=cpmtime();
      Bio bioU(bio_.m()+n, bio_.np(), bio_.xL(), bio_.xU()+bio_.d()*n);
      N m=bioU.m();
      N bm=bio_.m();
      VC rep(bioU.sz());
      Wf res(bioU,dyn_,rep); // res carries the time information of *this
      if (np()==1){
         for (N i=0;i<m;++i){
            res[i] = i<bm ? rep_[i] : C{};
         }
      }
      else if (np()==2){
         for (N i=0;i<m;++i){
            for (N j=0;j<m;++j){
               res.c(i,j) = (i<bm && j<bm) ? c(i,j) : C{};
            }
         }
      }
      else if (np()==3){
         for (N i=0;i<m;++i){
            for (N j=0;j<m;++j){
               for (N k=0;k<m;++k){
                  res.c(i,j,k) =
                  (i<bm && j<bm && k<bm) ? c(i,j,k) : C{};
               }
            }
         }
      }
      else if (np()==4){
         for (N i=0;i<m;++i){
            for (N j=0;j<m;++j){
               for (N k=0;k<m;++k){
                  for (N l=0;l<m;++l){
                     res.c(i,j,k,l) =
                     (i<bm && j<bm && k<bm && l<bm) ?
                        c(i,j,k,l) : C{};
                  }
               }
            }
         }
      }
      else cpmassert(np()<=4,"number of particles larger than 4");
      res.setAlarmU(false);
      res.readyU_=true;
      res.smooth(smoothSteps);
      R t2=cpmtime();
      cpmmessage(loc," done in "&cpm(t2-t1)&" s");
      CPM_MZ
      return res;
   }

   Wf Wf::extU(Wf const& a)const{
      Z mL=1;
      static Word loc("Wf::extU(Wf)");
      CPM_MA
      if (bioFix_){
         CPM_MZ
         return *this;
      }
      R t1=cpmtime();
      cpmassert(bio_.xU()==a.bio_.xL(),"no gap of biotopes allowed");
      cpmassert(bio_.d()==a.bio_.d(),
         "gap between neighbors should be constant");
      cpmassert(a.bio_.np()==1,"only one particle to be added");
      R noGBLim=15_R;

      N m=(*this).m()+a.m();
      N bm=bio_.m();
      N np=(*this).np()+1;
      R xL=(*this).xL();
      R xU=a.xU();
      Bio bioU(m,np,xL,xU);
      N sz=bioU.sz();
      N noBytes=sz*sizeC; // no: need of
      N szGB=1024_N*1024_N*1024_N; // so may bytes has a GB, my computer
         // offers 16.0 GB
      R noGB=noBytes*1_R/szGB;
    //  cpmassert(noGB<15_R,"too much RAM required");
      cout<<" need of GB = "<<noGB<<endl;
      if (noGB < noGBLim ){
         VC rep(sz);
         Wf res(bioU,dyn_,rep); // res carries the time
     // Wf res(Bio(m,np,xL,xU),dyn_,rep); // res carries the time
         // information of *this
         if (np==2){
            for (N i=0;i<m;++i){
               for (N j=0;j<m;++j){
                  C val{};
                  if (i<bm && j>=bm) val=rep_[i]*a[j-bm];
                  res.c(i,j) = val;
               }
            }
         }
         if (np==3){
            for (N i=0;i<m;++i){
               for (N j=0;j<m;++j){
                  for (N k=0;k<m;++k){
                     C val{};
                     if (i<bm && j<bm && k>=bm) val=c(i,j)*a[k-bm];
                     res.c(i,j,k) = val;
                  }
               }
            }
         }
         if (np==4){
            for (N i=0;i<m;++i){
               for (N j=0;j<m;++j){
                  for (N k=0;k<m;++k){
                     for (N l=0;l<m;++l){
                        C val{};
                        if (i<bm && j<bm && k<bm && l>=bm)
                           val=c(i,j,k)*a[l-bm];
                        res.c(i,j,k,l) = val;
                     }
                  }
               }
            }
         }
         else cpmassert(np<=4,"number of particles larger than 4");
         res.alarmL_=B(0);
         res.alarmU_=B(0);
         res.readyU_=true;
         res.smooth(smoothSteps);
         R t2=cpmtime();
         cpmmessage(loc," done in "&cpm(t2-t1)&" s");
         CPM_MZ
         return res;
      }
      else{
         dyn_.cyclic_=B(0);
         dyn_.tiny_=1.e12_R;
         bioFix_=B(1);
         if (bioFix_) cout<<"bioFix_ is true"<<endl;
         R t2=cpmtime();
         cpmmessage(loc," done in "&cpm(t2-t1)&" s");
         CPM_MZ
         return smooth(smoothSteps);
      }
   }

   Wf Wf::tenPrd(Wf const& a)const{
      Z mL=1;
      static Word loc("Wf Wf::tenPrd(Wf)");
      CPM_MA
      cpmassert(bio_.m()==a.bio_.m(),"biotope mismatch");
      cpmassert(bio_.xL()==a.bio_.xL(),"biotope mismatch");
      cpmassert(bio_.xU()==a.bio_.xU(),"biotope mismatch");
      cpmassert(a.bio_.np()==1,"only one particle to be added");

      Bio bioRes(m(),np()+1,xL(),xU());
      Wf res(bioRes,dyn_); // res carries the time information of *this
      if (res.np()==2){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               res.c(i,j) = rep_[i]*(a.rep_)[j];
            }
         }
      }
      if (res.np()==3){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  res.c(i,j,k) = c(i,j)*(a.rep_)[k];
               }
            }
         }
      }
      else if (res.np()==4){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  for (N l=0;l<m();++l){
                     res.c(i,j,k,l)=c(i,j,k)*(a.rep_)[l];
                  }
               }
            }
         }
      }
      else cpmassert(res.np()<=4,"number of particles larger than 4");
      CPM_MZ
      return res;
   }

bool Wf::alrL()const{
   Z mL=3;
   static Word loc("Wf::alrL()const");
   CPM_MA
   if (readyL_){
      CPM_MZ
      return alarmL_;}
   else{
      R tiny=dyn_.tiny_;
      N iL=0;
      alarmL_=false;
      R val;
      if (np()==1){
         val=rep_[iL].abs();
         if (val>tiny){ alarmL_=true;}
      }
      else if (np()==2){
         for (N i=0;i<m();++i){
            val=c(i,iL).abs();
            if (val>tiny){ alarmL_=true;}
         }
         for (N j=0;j<m();++j){
            val=c(iL,j).abs();
            if (val>tiny){ alarmL_=true;}
         }
      }
      else if (np()==3){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               val=c(i,j,iL).abs();
               if (val>tiny){ alarmL_=true; }
            }
         }
         for (N i=0;i<m();++i){
            for (N k=0;k<m();++k){
               val=c(i,iL,k).abs();
               if (val>tiny){ alarmL_=true;}
            }
         }
         for (N j=0;j<m();++j){
            for (N k=0;k<m();++k){
               val=c(iL,j,k).abs();
               if (val>tiny){ alarmL_=true; }
            }
         }
      }
      else if (np()==4){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  val=c(i,j,k,iL).abs();
                  if (val>tiny){ alarmL_=true;}
               }
            }
         }
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N l=0;l<m();++l){
                  val=c(i,j,iL,l).abs();
                  if (val>tiny){ alarmL_=true;}
               }
            }
         }
         for (N i=0;i<m();++i){
            for (N k=0;k<m();++k){
               for (N l=0;l<m();++l){
                  val=c(i,iL,k,l).abs();
                  if (val>tiny){ alarmL_=true;}
               }
            }
         }
         for (N j=0;j<m();++j){
            for (N k=0;k<m();++k){
               for (N l=0;l<m();++l){
                  val=c(iL,j,k,l).abs();
                  if (val>tiny){ alarmL_=true;}
               }
            }
         }
      }
      readyL_=true;
      CPM_MZ
      return alarmL_;
   }
}

bool Wf::alrU()const{
   Z mL=3;
   static Word loc("Wf::alrU()const");
   CPM_MA
   if (readyU_){
      CPM_MZ
      return alarmU_;
   }
   else{
      R tiny=dyn_.tiny_;
      N iU=m()-1;
      alarmU_=false;
      R val;
      if (np()==1){
         val=rep_[iU].abs();
         if (val>tiny){ alarmU_=true;}
      }
      else if (np()==2){
         for (N i=0;i<m();++i){
            val=c(i,iU).abs();
            if (val>tiny){ alarmU_=true;}
         }
         for (N j=0;j<m();++j){
            val=c(iU,j).abs();
            if (val>tiny){ alarmU_=true;}
         }
      }
      else if (np()==3){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               val=c(i,j,iU).abs();
               if (val>tiny){ alarmU_=true;}
            }
         }
         for (N i=0;i<m();++i){
            for (N k=0;k<m();++k){
               val=c(i,iU,k).abs();
               if (val>tiny){ alarmU_=true;}
            }
         }
         for (N j=0;j<m();++j){
            for (N k=0;k<m();++k){
               val=c(iU,j,k).abs();
               if (val>tiny){ alarmU_=true;}
            }
         }
      }
      else if (np()==4){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  val=c(i,j,k,iU).abs();
                  if (val>tiny){ alarmU_=true; }
               }
            }
         }
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N l=0;l<m();++l){
                  val=c(i,j,iU,l).abs();
                  if (val>tiny){ alarmU_=true; }
               }
            }
         }
         for (N i=0;i<m();++i){
            for (N k=0;k<m();++k){
               for (N l=0;l<m();++l){
                  val=c(i,iU,k,l).abs();
                  if (val>tiny){ alarmU_=true;}
               }
            }
         }
         for (N j=0;j<m();++j){
            for (N k=0;k<m();++k){
               for (N l=0;l<m();++l){
                  val=c(iU,j,k,l).abs();
                  if (val>tiny){ alarmU_=true;}
               }
            }
         }
      }
      readyU_=true;
      CPM_MZ
      return alarmU_;
   }
}

   Wf Wf::hf()const{
      static Z beat1=0;
      alarmL_=false;
      alarmU_=false;
      readyL_=false;
      readyU_=false;
      R sumL{0.};
      R sumU{0.};
      R t1=cpmtime();
      Z cp=-1; // communication pane
    //  bool going=!dyn_.cyclic_; // h() is called for performing a time step
         // and not for working within norm determination
      cpmvalue("beat1: ",beat1,2);
      R fac=1_R/(2_R*d()*d()*mass());
      R tiny=dyn_.tiny_;
      N iU=m()-1;
      Wf res(*this);
      if (np()==1){
         R val=rep_[0].abs(); // avoid a second evaluation of abs() in
            // the block 'if (going){...}'
         if (val>tiny){
            sumL+=val;
            alarmL_=true;
         }
         val=rep_[iU].abs();
         if (val>tiny){
            alarmU_=true;
         }
         for (N i=0;i<m();++i){
            C amC = rep_[i]*R{2.}, amUi, amLi;
            if (dyn_.cyclic_){
               amUi= i+1 < m() ? rep_[i+1] : rep_[0];
               amLi= Z(i)-1 >=0 ? rep_[i-1] : rep_[iU];
            }
            else{
               amUi= i+1 < m() ? rep_[i+1] : C{};
               amLi= Z(i)-1 >=0 ? rep_[i-1] : C{};
            }
               // i-1 give unregular values for i of type N.
            res[i]=(amC-amUi-amLi)*fac;
         }
      }
      else if (np()==2){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
           //    VN vij({i,j});
               R val;
               if (i==0){
                  val=c(i,j).abs();
                  if (val>tiny){
                     alarmL_=true;
                  }
               }
               if (i==iU){
                  val=c(i,j).abs();
                  if (val>tiny){
                     alarmU_=true;
                  }
               }
               if (j==0){
                  val=c(i,j).abs();
                  if (val>tiny){
                     alarmL_=true;
                  }
               }
               if (j==iU){
                  val=c(i,j).abs();
                  if (val>tiny){
                     alarmU_=true;
                  }
               }
               C amC = c(i,j)*R{4.},amLi,amLj,amUi,amUj;
               if (dyn_.cyclic_){
                  amUi = i+1 < m() ? c(i+1,j) : c(0,j);
                  amLi = Z(i)-1 >= 0 ? c(i-1,j) : c(iU,j);
                  amUj = j+1 < m() ? c(i,j+1) : c(i,0);
                  amLj = Z(i)-1 >= 0 ? c(i,j-1) : c(i,iU);
               }
               else{
                  amUi = i+1 < m() ? c(i+1,j) : C{};
                  amLi = Z(i)-1 >= 0 ? c(i-1,j) : C{};
                  amUj = j+1 < m() ? c(i,j+1) : C{};
                  amLj = Z(i)-1 >= 0 ? c(i,j-1) : C{};
               }
               res.c(i,j)=(amC-amUi-amLi-amUj-amLj)*fac ;
            }
         }
      }
      else if (np()==3){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  VN vijk({i,j,k});
                  R val;
                  if (i==0){
                     val=c(i,j,k).abs();
                     if (val>tiny){
                        alarmL_=true;
                     }
                  }
                  if (i==iU){
                     val=c(i,j,k).abs();
                     if (val>tiny){
                        alarmU_=true;
                     }
                  }
                  if (j==0){
                     val=c(i,j,k).abs();
                     if (val>tiny){
                        alarmL_=true;
                     }
                  }
                  if (j==iU){
                     val=c(i,j,k).abs();
                     if (val>tiny){
                        alarmU_=true;
                     }
                  }
                  if (k==0){
                     val=c(i,j,k).abs();
                     if (val>tiny){
                        alarmL_=true;
                     }
                  }
                  if (k==iU){
                     val=c(i,j,k).abs();
                     if (val>tiny){
                        alarmU_=true;
                     }
                  }
                  C amC = c(i,j,k)*R{6.};
                  C amLi,amLj,amLk,amUi,amUj,amUk;
                  if (dyn_.cyclic_){
                     amUi = i+1 < m() ? c(i+1,j,k) : c(0,j,k);
                     amLi = Z(i)-1 >=0  ? c(i-1,j,k) : c(iU,j,k);
                     amUj = j+1 < m() ? c(i,j+1,k) : c(i,0,k);
                     amLj = Z(j)-1 >=0  ? c(i,j-1,k) : c(i,iU,k);
                     amUk = k+1 < m() ? c(i,j,k+1) : c(i,j,0);
                     amLk = Z(k)-1 >=0  ? c(i,j,k-1) : c(i,j,iU);
                  }
                  else{
                     amC = c(i,j,k)*R{6.};
                     amUi = i+1 < m() ? c(i+1,j,k) : C{};
                     amLi = Z(i)-1 >=0  ? c(i-1,j,k) : C{};
                     amUj = j+1 < m() ? c(i,j+1,k) : C{};
                     amLj = Z(j)-1 >=0  ? c(i,j-1,k) : C{};
                     amUk = k+1 < m() ? c(i,j,k+1) : C{};
                     amLk = Z(k)-1 >=0  ? c(i,j,k-1) : C{};
                  }
                  res.c(i,j,k)=(amC-amUi-amLi-amUj-amLj-amUk-amLk)*fac;
               }
            }
         }
      }
      else if (np()==4){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  for (N l=0;l<m();++l){
                     VN vijkl({i,j,k,l});
                     R val;
                     if (i==0){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmL_=true;
                        }
                     }
                     if (i==iU){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmU_=true;
                        }
                     }
                     if (j==0){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmL_=true;
                        }
                     }
                     if (j==iU){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmU_=true;
                        }
                     }
                     if (k==0){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmL_=true;
                        }
                     }
                     if (k==iU){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmU_=true;
                        }
                     }
                     if (l==0){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmL_=true;
                        }
                     }
                     if (l==iU){
                        val=c(i,j,k,l).abs();
                        if (val>tiny){
                           alarmU_=true;
                        }
                     }

                     C amC = c(i,j,k,l)*R{8.};
                     C amLi,amLj,amLk,amLl,amUi,amUj,amUk,amUl;
                     if (dyn_.cyclic_){
                        amUi = i+1 < m() ? c(i+1,j,k,l) : c(0,j,k);
                        amLi = Z(i)-1 >=0  ? c(i-1,j,k,l) : c(iU,j,k);
                        amUj = j+1 < m() ? c(i,j+1,k,l) : c(i,0,k);
                        amLj = Z(j)-1 >=0  ? c(i,j-1,k,l) : c(i,iU,k);
                        amUk = k+1 < m() ? c(i,j,k+1,l) : c(i,j,0);
                        amLk = Z(k)-1 >=0  ? c(i,j,k-1,l) : c(i,j,iU);
                        amUl = l+1 < m() ? c(i,j,k,l+1) : c(i,j,0);
                        amLl = Z(l)-1 >=0  ? c(i,j,k,l-1) : c(i,j,iU);
                     }
                     else{
                        amC = c(i,j,k,l)*R{8.};
                        amUi = i+1 < m() ? c(i+1,j,k,l) : C{};
                        amLi = Z(i)-1 >=0  ? c(i-1,j,k,l) : C{};
                        amUj = j+1 < m() ? c(i,j+1,k,l) : C{};
                        amLj = Z(j)-1 >=0  ? c(i,j-1,k,l) : C{};
                        amUk = k+1 < m() ? c(i,j,k+1,l) : C{};
                        amLk = Z(k)-1 >=0  ? c(i,j,k-1,l) : C{};
                        amUl = l+1 < m() ? c(i,j,k,l+1) : C{};
                        amLl = Z(l)-1 >=0  ? c(i,j,k,l-1) : C{};
                     }
                     res.c(i,j,k,l)=(amC-amUi-amLi-amUj-amLj-amUk-amLk
                        -amUl-amLl)*fac;
                  }
               }
            }
         }
      }
      else cpmassert(np()<=4,"Wf::h(): number of particles larger than 4");
      beat1++;
      R t2=cpmtime();
      R realTimeSpan=t2-t1;
      R beatRate=R{60.}/realTimeSpan;
      R speed = dt_< R{0.} ? R{0.} : dt_*R{1000.}/realTimeSpan;
      Z bioSize=m();
      Z nP=np();
      if (dt_>R{0.}){
         cpmvalues2("beats/min, speed: ",beatRate,speed,3);
      }
      else{
         cpmvalue("beats/min: ",beatRate,3);
      }
      cpmvalues2("bioSize, np: ",bioSize,nP,4);
      readyL_=true;
      readyU_=true;
      return res;
   }

   Wf Wf::smoothS()const{
      Z mL=1;
      static Word loc("Wf::smoothS()const");
      CPM_MA
      auto fac=0.25_R;
      Wf res(*this);
      if (np()==1){
         for (N i=0;i<m();++i){
            C amC = rep_[i]*2_R;
            C amUi= i+1 < m() ? rep_[i+1] : C{};
            C amLi= Z(i)-1 >=0 ? rep_[i-1] : C{};
           res[i]=(amC+amUi+amLi)*fac;
         }
      }
      else if (np()==2){
         fac=1_R/8_R;
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               C amC,amLi,amLj,amUi,amUj;
               amC = c(i,j)*4_R;
               amUi = i+1 < m() ? c(i+1,j) : C{};
               amLi = Z(i)-1 >= 0 ? c(i-1,j) : C{};
               amUj = j+1 < m() ? c(i,j+1) : C{};
               amLj = Z(i)-1 >= 0 ? c(i,j-1) : C{};
               res.c(i,j)=(amC+amUi+amLi+amUj+amLj)*fac ;
            }
         }
      }
      else if (np()==3){
         fac=1_R/12_R;
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  C amC,amLi,amLj,amLk,amUi,amUj,amUk;
                  amC = c(i,j,k)*6_R;
                  amUi = i+1 < m() ? c(i+1,j,k) : C{};
                  amLi = Z(i)-1 >=0  ? c(i-1,j,k) : C{};
                  amUj = j+1 < m() ? c(i,j+1,k) : C{};
                  amLj = Z(j)-1 >=0  ? c(i,j-1,k) : C{};
                  amUk = k+1 < m() ? c(i,j,k+1) : C{};
                  amLk = Z(k)-1 >=0  ? c(i,j,k-1) : C{};
                  res.c(i,j,k)=(amC+amUi+amLi+amUj+amLj-amUk-amLk)*fac;
               }
            }
         }
      }
      else if (np()==4){
         fac=1_R/16_R;
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  for (N l=0;l<m();++l){
                     C amC,amLi,amLj,amLk,amLl,amUi,amUj,amUk,amUl;
                     amC = c(i,j,k,l)*8_R;
                     amUi = i+1 < m() ? c(i+1,j,k,l) : C{};
                     amLi = Z(i)-1 >=0  ? c(i-1,j,k,l) : C{};
                     amUj = j+1 < m() ? c(i,j+1,k,l) : C{};
                     amLj = Z(j)-1 >=0  ? c(i,j-1,k,l) : C{};
                     amUk = k+1 < m() ? c(i,j,k+1,l) : C{};
                     amLk = Z(k)-1 >=0  ? c(i,j,k-1,l) : C{};
                     amUl = l+1 < m() ? c(i,j,k,l+1) : C{};
                     amLl = Z(l)-1 >=0  ? c(i,j,k,l-1) : C{};
                     res.c(i,j,k,l)=
                        (amC+amUi+amLi+amUj+amLj+amUk+amLk+amUl+amLl)*fac;
                  }
               }
            }
         }
      }
      else cpmassert(np()<=4,"Wf::smooth(): number of particles larger than 4");
      CPM_MZ
      return res;
   }

   Wf Wf::hl()const{
      static Z beat2=0;
      R t1=cpmtime();
      cpmvalue("beat2: ",beat2,2);
      R fac=1_R/(2_R*d()*d()*mass());
      N iU=m()-1;
      Wf res(*this);
      if (np()==1){
         for (N i=0;i<m();++i){
            C amC = rep_[i]*2_R, amUi, amLi;
            amUi= i+1 < m() ? rep_[i+1] : C{};
            amLi= Z(i)-1 >=0 ? rep_[i-1] : C{};
           res[i]=(amC-amUi-amLi)*fac;
         }
      }
      else if (np()==2){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               C amC = c(i,j)*4_R,amLi,amLj,amUi,amUj;
               amUi = i+1 < m() ? c(i+1,j) : C{};
               amLi = Z(i)-1 >= 0 ? c(i-1,j) : C{};
               amUj = j+1 < m() ? c(i,j+1) : C{};
               amLj = Z(i)-1 >= 0 ? c(i,j-1) : C{};
               res.c(i,j)=(amC-amUi-amLi-amUj-amLj)*fac ;
            }
         }
      }
      else if (np()==3){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  C amC,amLi,amLj,amLk,amUi,amUj,amUk;
                  amC = c(i,j,k)*R{6.};
                  amUi = i+1 < m() ? c(i+1,j,k) : C{};
                  amLi = Z(i)-1 >=0  ? c(i-1,j,k) : C{};
                  amUj = j+1 < m() ? c(i,j+1,k) : C{};
                  amLj = Z(j)-1 >=0  ? c(i,j-1,k) : C{};
                  amUk = k+1 < m() ? c(i,j,k+1) : C{};
                  amLk = Z(k)-1 >=0  ? c(i,j,k-1) : C{};
                  res.c(i,j,k)=(amC-amUi-amLi-amUj-amLj-amUk-amLk)*fac;
               }
            }
         }
      }
      else if (np()==4){
         for (N i=0;i<m();++i){
            for (N j=0;j<m();++j){
               for (N k=0;k<m();++k){
                  for (N l=0;l<m();++l){
                     C amC,amLi,amLj,amLk,amLl,amUi,amUj,amUk,amUl;
                     amC = c(i,j,k,l)*R{8.};
                     amUi = i+1 < m() ? c(i+1,j,k,l) : C{};
                     amLi = Z(i)-1 >=0  ? c(i-1,j,k,l) : C{};
                     amUj = j+1 < m() ? c(i,j+1,k,l) : C{};
                     amLj = Z(j)-1 >=0  ? c(i,j-1,k,l) : C{};
                     amUk = k+1 < m() ? c(i,j,k+1,l) : C{};
                     amLk = Z(k)-1 >=0  ? c(i,j,k-1,l) : C{};
                     amUl = l+1 < m() ? c(i,j,k,l+1) : C{};
                     amLl = Z(l)-1 >=0  ? c(i,j,k,l-1) : C{};
                     res.c(i,j,k,l)=(amC-amUi-amLi-amUj-amLj-amUk-amLk
                        -amUl-amLl)*fac;
                  }
               }
            }
         }
      }
      else cpmassert(np()<=4,"Wf::h(): number of particles larger than 4");
      beat2++;
      R t2=cpmtime();
      R realTimeSpan=t2-t1;
      R beatRate=R{60.}/realTimeSpan;
      R speed = dt_< 0_R ? 0_R : dt_*1000_R/realTimeSpan;
      Z bioSize=m();
      Z nP=np();
      if (dt_>0_R){
         cpmvalues2("beats/min, speed: ",beatRate,speed,3);
      }
      else{
         cpmvalue("beats/min: ",beatRate,3);
      }
      cpmvalues2("bioSize, np: ",bioSize,nP,4);
      return res;
   }


   Wf Wf::v()const{
      VC rep(rep_);
      for (N i=0;i<sz();++i){
         VN vi=toMultInd(i,ind());
         R pi=0.;
         for (N j=0;j<np();++j) pi+=dyn_.pot_(bio_.x()[vi[j]]);
         rep[i]*=pi;
      }
      return modify(rep);
   }

   // potental energy by pair interaction
   Wf Wf::vp()const{
      VC rep(rep_);
      for (N i=0;i<sz();++i){
         VN vi=toMultInd(i,ind());
         R pi=0.;
         for (N j=0; j<np(); ++j){
            for (N k=0; k<j; ++k){
               pi+=dyn_.potPair_(cpmabs(bio_.x()[vi[j]]-bio_.x()[vi[k]]));
            }
         }
         rep[i]*=pi;
      }
      return modify(rep);
   }

   void Wf::hNorm_(Z iMax, R tWait){
      if (bioFix_) return;
      Z mL=1;
      Word loc("void Wf::hNorm_(Z,R)");
      CPM_MA
      Graph gr;
      Wf psi(bio_,dyn_,C(1_R,1_R));
      psi=psi.ran();
      psi.nor_();
      if (tWait>0_R){
         psi.mark(gr);
         gr.vis();
         cpmwait(tWait,2);
      }
      R eta=0_R;
      R_Vector etaV(iMax);
    //  bool cycMem=dyn_.cyclic_;
      dyn_.cyclic_=true; // algorithm for operator norm requires an
         // hermitic operator. For running extensible dynamics cyclic
         // should be false (default value of B).
      for (Z i=1;i<=iMax;++i){
         psi.cyclic(true);
         psi=psi.ham(); // total hamiltonian (with potentials)
         eta=psi.nor_();
         if (tWait>0_R){
            psi.mark(gr);
            gr.vis();
            cpmwait(tWait,2);
         }
         ostringstream ost;
         ost<<" i= "<<i<<", eta ="<<eta<<endl;
         cpmmessage(ost,-1,true);
         etaV[i]=eta;
      }
      if (tWait>0_R){
         gr.mark(etaV);
         gr.vis();
         cpmwait(tWait,2);
      }
    //  dyn_.cyclic_=cycMem;
      R res=etaV[iMax];
      R res2=hKinNorm();
      ostringstream ost;
      ost<<"hNorm = "<<res<<endl;
      ost<<"hKinNorm = "<<res2<<endl;
      ost<<"number of particles = "<<np()<<endl;
      ost<<"size of *this in MB: "<<psi.memUse()/(1024.*1024)<<endl;
      cpmmessage(ost,-1,true);
      alarmL_=false;
      alarmU_=false;
      CPM_MZ
      normH_=res;
   }

   R Wf::hKinNorm()const{
      R fac=R{1.}/(R{2.}*d()*d()*mass());
      return fac*R{4.};
   }

   State WfToState(Wf const& wf){
      cpmassert(wf.np()==1,"WfToState(Wf)")
      State res((Z)wf.sz());
      for (N i=0;i<wf.sz();++i) res[(Z)i+1]=wf[i];
      return res;
   }

   Wf StateToWf(State const& s, Bio const& bio, Dyn const& dyn){
      Wf res(bio, dyn);
      for (Z i=1; i<=s.dim();++i) res[(N)(i-1)]=s[i];
      return res;
   }

////////////////////////// class Sys //////////////////////////////////////
// Strategy for the 'essential operators' nach BS A Tour of C++ Third Edition
// Chapter 6, 'rule of zero'.

class Sys{
   R t_{0.};
   Wf x_;
   Wf v_;
   R dt_{0.};
public:
   Sys(){}
      // initialize all data members for which there is no default constructor.
      // Then the the default destructor can always be defined as above.
// Since we won't derive from Sys we get all 'essential operators' made by the
// compiler. If we would be willing to derive from Sys, we had to add the
// following declarations: 'rule of five'. My plan is to prefer templates
// over derivations.
//   virtual ~Sys() = default;
//   Sys(const Sys&)=default;
//   Sys(Sys&&)=default;
//   Sys& operator=(const Sys&)=default;
//   Sys& operator=(Sys&&)=default;
   Sys(Wf const& wf, R dtRel, R t=0.):t_{t}, x_{wf}, v_{x_.dot()},
      dt_{x_.getTimeStep(dtRel)}{}
   R t()const{ return t_;}
   void set_t(R t){t_=t;x_.set_t(t);v_.set_t(t);}
   Wf const& x()const{ return x_;}
   Wf& x(){ return x_;}
   Wf const& v()const{ return v_;}
   R dt()const{ return dt_;}
   R xL()const{ return x_.xL();}
   R xU()const{ return x_.xU();}
   bool free()const{ return x_.free();}
   bool freeL()const{ return x_.freeL();}
   bool freeU()const{ return x_.freeU();}
   bool isCyclic()const{ return x_.isCyclic();}
   void step_();
   void step0_();
   void mark(Graph& gr)const{
      gr.addText("t="&cpm(t_)&", size="&cpm(x_.sz())&", np="&cpm(x_.np()));
      x_.mark(gr);
   }
   void status(Graph& gr, R tWait=6.)const;
   void movToRim_(CyclicAlarm& ca,Graph& gr,B const& makeMovie,R tWait);
   N m()const{ return x_.m();}
   N np()const{ return x_.np();}
   N sz()const{ return x_.sz();}
   N memUse()const{ return x_.memUse()*2+sizeR*2;}
};

 void Sys::step_(){
   Z mL=1;
   static Word loc("Sys::step_()");
   CPM_MA
   R tau=dt_*0.25_R;
   R tau2=tau*2_R;
   t_+=tau;
   x_+=v_*tau;
   v_+=(x_.dot()-v_)*2_R;
   x_+=v_*tau2;
   t_+=tau2;
   v_+=(x_.dot()-v_)*2_R;
   x_+=v_*tau;
   t_+=tau;
   set_t(t_); // tansfers the new value of t_ also to the dyn_
      // objects contained in the Wf instances x_ and v_.
   CPM_MZ
 }

 void Sys::status(Graph& gr, R tWait)const{
   ostringstream ost;
   ost<<"******  status of Sys *************"<<endl;
   ost<<"t = "<<t()<<endl;
   ost<<"sz = "<<sz()<<endl;
   ost<<"np = "<<np()<<endl;
   ost<<"m = "<<m()<<endl;
   ost<<"alarmL = "<<x().alarmL()<<endl;
   ost<<"alarmU = "<<x().alarmU()<<endl;
   cpmmessage(ost,-1,true);
   mark(gr);
   gr.vis();
   cpmwait(tWait,2);
 }

void Sys::movToRim_(CyclicAlarm& ca,Graph& gr,B const& makeMovie,R tWait)
 {
   Z mL=1;
   Word loc( "Sys::movToRim_(...)" );
   CPM_MA
   if (!free()){
      status(gr,tWait);
      CPM_MZ
      return;
   }
   else{
      while (free()){
         step_(); // does not change sys.np()
         if ( ca()){
            mark(gr); // notes sys.t() and sys.sz() on the image
            gr.vis(makeMovie);
         }
      }
      status(gr,tWait);
      CPM_MZ
      return;
   }
 }

} // namespace eqm2;

using namespace eqm2;

// a full cycle of expanding quantum dynamics
// This seems to be the newer version. The new attempt of 2023-02-20
// starts from here. There is no test1() anymore.
// It became  replaced by test2().
void test2(){
   cout<<"test2() started"<<endl;
   Z mL=1;
   Word loc("eqm2::test2()");
   CPM_MA
// data to be read from "test2.ini"
// values set here may be overwritten by the values in "test2.ini"
   auto m=10_Z;
   auto xL=-1_R;
   auto xU=1_R;
   auto xLMax=-1000_R;
   auto xUMax=1000_R;
   auto mass=1_R;
   auto rhoRel=0.5_R; // all interactions have a reach rhoRel*(xU-xL)/2
   auto tiny=1.e-12_R;
   auto soft=.01_R;
   auto axisInspect=0_R;
   auto smoothSteps=0_N;
   auto lambdaPotRel=1_R;
   auto lambdaPotPairRel=1_R;
   auto vRel=0.25_R;
   auto tMax=0.01_R;
   auto dtRel=0.25_R;
   auto tWait1=1_R;
   auto tWait2=1_R;
   auto tWait3=1_R;
   auto tWait4=1_R;
   auto gamma=0.1_R;
   auto szMax=64_N;
   auto npMax=4_N;
   auto makeMovie=0_B;
   auto useReduction=0_B;
   auto imgPer=2_Z;
   RecordHandler rch("test2.ini");
   Word sec="system data";
   cpmrh(m);
   cpmrh(xL);
   cpmrh(xU);
   cpmrh(xLMax);
   cpmrh(xUMax);
   cpmrh(mass);
   cpmrh(rhoRel);
   cpmrh(tiny);
   cpmrh(soft);
   cpmrh(axisInspect);
   cpmrh(smoothSteps);
   cpmrh(lambdaPotRel);
   cpmrh(lambdaPotPairRel);
   cpmrh(vRel);
   sec="run control";// m() instead of sz() would not be ????
   cpmrh(tMax);
   cpmrh(dtRel);
   cpmrh(tWait1);
   cpmrh(tWait2);
   cpmrh(tWait3);
   cpmrh(tWait4);
   cpmrh(szMax);
   cpmrh(npMax);
   cpmrh(gamma);
   cpmrh(makeMovie);
   cpmrh(useReduction);
   cpmrh(imgPer);
// setting the initial biotope
   const Iv iv(xL,xU);
   const R shiftRange=iv.abs();
   const R gap=shiftRange/R(m);
//// providing a reduced interval to which the potential is confined (???)
// providing a reduced interval to which the initial wave function is confined


   const R potRange=iv.abs();
   const R xc=iv.cen();
   const R halfRange=potRange*0.5_R;
   const R rho=rhoRel*halfRange;
   const R rhoInv=1_R/rho;
   const R halfRangeInv=1_R/halfRange;
   QMBase::setUnits("natural",iv);
   QMBase::setMass(mass); // now all instances of of class State have the
      // mass as given by the input file test2.ini
// adapting the view port to the initial biotope

   Graph gr;
   gr.setX(iv);

// defining a potential well. The definition does no longer assume xL==-1
// and xU==+1.
   function<R(R)> fpot{[&xc,&halfRangeInv,&rhoRel](R x){ // defining a
   //function in terms of a lambda expression.
      R eta=(x-xc)*halfRangeInv;
         // eta has the value 0 at the midpoint of iv and the value 1 and
         // -1 at the rims. eta is a variable dependend on x. It also is the
         // value of x in a different coordinate system.
         // The value of fpot at x==0 is thus -rho2.
      R eta2{eta*eta};
      R rho2{rhoRel*rhoRel};
      return eta2 > rho2 ? 0_R : eta2-rho2;
   }};
   R_Func pot=fpot;
   R_Func pot0=fpot;
   R_Func potVac; // pot with constant value 0
   R eKinMax=QMBase::eMax(m);
   R eKinSel=eKinMax*vRel*vRel;
   R vMax=QMBase::vMax(m);
   R vb=vMax*vRel; // v boost
   R lambda=-lambdaPotRel*eKinSel/pot(0_R);
   gr.addText("function pot");
   gr.mark(pot);
   gr.visRep(makeMovie);
   cpmwait(tWait2,2);
   gr.clearText();
   gr.clr_();
   gr.paint(LIGHTGREEN);

   F<State,State> h(function<State(State)>{[&pot,&lambda,&iv]
   (State const& s){
      return s.sinParHam2(pot*lambda,iv.abs());
   }});
   Opr op((Z)m,h);
   Observable ob=op.toObservable();
   gr.setText("hamilton operator as matrix",0.1,0.2,WHITE);
   ob.mark(gr); // potentially changes gr
   gr.visRep(makeMovie);
   cpmwait(tWait3,2);
   gr.clearText();

   SpectralObject so=ob.spectralRepresentation();
   gr.addText("spectrum");
   so.markSpectrum(gr); // potentially changes gr
   gr.visRep(makeMovie);
   cpmwait(tWait3,2);
   gr.clearText();

   State groundstate=so[1];
   groundstate.mark(gr); // potentially changes gr
   gr.setText("raw groundstate",0.1,0.2);
   gr.visRep(makeMovie);
   cpmwait(tWait3,2);
   gr.clearText();

   gr.setX(iv); // retting shape of potentially changed gr

   R_Func winFunc=R_Func::tableMountain(iv,soft);
   gr.mark(winFunc);
   gr.setText("window function",0.1,0.2);
   gr.visRep(makeMovie);
   cpmwait(tWait3,2);
   gr.clearText();

   groundstate=groundstate.windowing(soft).normalize();
   groundstate.mark(gr);
   gr.setText("windowed groundstate",0.1,0.2);
   gr.visRep(makeMovie);
   cpmwait(tWait3,2);
   gr.clearText();

   R enrMin=so(1).real();
   R potMin=lambda*pot(0_R);
   R eTot=enrMin+eKinSel;
   ostringstream o0;
   o0<<"enrMin = "<<enrMin<<", eKinSel = "<<eKinSel<<", eTot = "<<eTot<<endl;
   cpmmessage(o0,1,true);

   function<R(R)> fpotPair{[&rhoInv](R x){ // x distance between
      // particles
      R y=cpmabs(x)*rhoInv;
      return y > 1 ? 0_R : (cpmcos(y*cpmpi)+1_R)*0.5_R;
   }};
   R_Func potPairPrelim=fpotPair;
   R_Func potPair=potPairPrelim*(lambdaPotPairRel*cpmabs(potMin));
   gr.clr_();
   gr.setX(iv);
   gr.mark(potPair);
   gr.setText("pair potential",0.1,0.2);
   gr.visRep(makeMovie);
   cpmwait(tWait4,2);
   gr.clearText();

// Now we start working with Wf instead of State.
// Our first state to be considered is a single particle state, so that
// potPair does not yet come in action.

   Bio bio(m,1_N,xL,xU,Iv(xLMax,xUMax));
   ostringstream ostBio;
   bio.write(ostBio);
   if (axisInspect > 0_R) gr.setY_(Iv(-axisInspect,axisInspect));
   Dyn dyn(mass,tiny,potVac,potPair);
   Wf wgs(bio,dyn); // wgs: wavefunction ground state, an instance of Wf
   wgs.setSmoothSteps_(smoothSteps);
   for (N i=0,j=1; j<=m; ++i,++j){
      wgs[i]=groundstate[j]; // now we no longer have vectors the indexing of
      // which starts with 1
   }
   wgs.mark(gr);
   gr.setText("initial wave function at rest",0.1,0.2);
   gr.visRep(makeMovie);
   cpmwait(tWait4,2);
   gr.clearText();

   Wf wf0=wgs.boost(vb);
   cpmmessage("initial wave function",1,true);
   wf0.mark(gr,iv); // initial wave function
   gr.setText("initial wave function moving",0.1,0.2);
   gr.visRep(makeMovie);
   cpmwait(tWait4,2);
   gr.clearText();

// now wf0 becomes part of a dynamical system with a defined time
// evolution.
   auto t_ini=0_R;
   Sys sys(wf0,dtRel,t_ini);
      // now the hamiltonian of Wf is responsible for the time evolution
   cout<<"sys.t()="<<sys.t()<<endl;
   sys.mark(gr);
   gr.visRep(makeMovie);
   cpmwait(tWait4,2);
   gr.clearText();
   ostringstream ostIni;
   ostIni<<" initial sys created "<<endl;
   cpmmessage(ostIni,1,true);
   CyclicAlarm ca(imgPer);
static Z procCount=0_Z;
proceed:
   procCount++;
   cout<<" procCount = "<<procCount<<" sys.t() = "<<sys.t()<<endl;
   cpmmessage("label 'proceed' passed, procCount = "&cpm(procCount),1,true);
   ostringstream ost0;
   ost0<<"*********** while loop ahead ***********"<<endl;
   ost0<<"sys.t()="<<sys.t()<<endl; // here 0 is OK for procCount == 1
   ost0<<"sys.dt()="<<sys.dt()<<endl;
   ost0<<"sys.xL()="<<sys.xL()<<endl;
   ost0<<"sys.xU()="<<sys.xU()<<endl;
   ost0<<"sys.m()="<<sys.m()<<endl;
   ost0<<"sys.np()="<<sys.np()<<endl;
   ost0<<"sys.sz()="<<sys.sz()<<endl;
   ost0<<"memUse="<<sys.memUse()<<" bytes"<<endl;
   ost0<<"sys.freeL()="<<sys.freeL()<<endl;
   ost0<<"sys.freeU()="<<sys.freeU()<<endl;
   cpmmessage(ost0,-1,true);
   while (sys.free() && sys.t()<=tMax && sys.m()<=szMax && sys.np()<=npMax)
   {
      sys.step_(); // only if this is done, sys.freeT() is meaningful
      if ( ca()){
         sys.mark(gr);
         gr.vis(makeMovie);
      }
   }

   ostringstream ost;
   ost<<"******  while loop done *************"<<endl;
   ost<<"****** alarmL = "<<sys.x().alarmL()<<endl;
   ost<<"****** alarmU = "<<sys.x().alarmU()<<endl;
   ost<<"After encountering rim, or passing time limit, or size limit"
   <<endl;
   cpmmessage(ost,-1,true);
   sys.mark(gr);
   gr.vis(makeMovie);
   cpmwait(tWait4,2);

   if (sys.t()>tMax){ // this is condition for finishing the program run
      cpmmessage("tMax exceeded: tMax="&cpm(tMax)&", t="&cpm(sys.t()));
      if (makeMovie) ppm2mp4(Word("video"));
      CPM_MZ
      return;
   }

   if (sys.sz()>szMax){ // we have to reduce size by reduction
      cpmmessage("szMax exceeded: szMax="&cpm(szMax)&
      ", sys.sz()="&cpm(sys.sz()));
      goto reduction; // reduces np to 1 and thus considerably reduces sz
   }
   if (sys.np()>npMax){ // we have to reduce size by reduction
      cpmmessage("np exceeded: npMax="&cpm(npMax)&
      ", sys.np()="&cpm(sys.np()));
      if (useReduction) goto reduction;
         // reduces np to 1 and thus considerably reduces sz
      else{
         if (makeMovie) ppm2mp4(Word("video"));
         CPM_MZ
         return;
      }
   }
// we reached the rim of the biotope and have to extend it. Since the size
// limit is not yet reached we may try this.
   Wf wf=sys.x();
   R t1=sys.t();
   cout<<"t1 = "<<t1<<" should not be 0"<<endl;
   if (!wf.freeL()){
      ostringstream ost1;
      ost1<<"if (!wf.freeL())"<<endl;
      cpmmessage(ost1,1,true);
      wf=wf.extL(m);
      N ws=wf.sz();
      if (ws>szMax){
         cpmmessage("szMax exceeded: szMax="&cpm(szMax)&
         ", wf.sz()="&cpm(ws));
         if (useReduction) goto reduction;
         // reduces np to 1 and thus considerably reduces sz
         else{
            if (makeMovie) ppm2mp4(Word("video"));
            CPM_MZ
            return;
         }
      }
      cpmmessage("new Sys ahead",1,true);
      sys=Sys(wf,dtRel,t1);
      cpmmessage("new Sys done",1,true);
      ostringstream ost2;
      ost2<<"if (!wf.freeL()) done"<<endl;
      cpmmessage(ost1,1,true);
   }
   //else if (!wf.freeU()){
   // this mistaken 'else' gives correct results in the probable case that only
   // one of the two functions freeU and freeL give the value 'false'. If both
   // do so the extension at the upper end is erroneously omitted. This error
   // was corrected 20232-03-19.
   if (!wf.freeU()){
      ostringstream ost1;
      ost1<<"if (!sys.freeU())"<<endl;
      cpmmessage(ost1,1,true);
      if (wf.np()==npMax){
         wf=wf.extU(m); // extending the biotope without adding a particle
         N ws=wf.sz();
         if (ws>szMax){
            cpmmessage("szMax exceeded: szMax="&cpm(szMax)&
            ", wf.sz()="&cpm(ws));
            if (useReduction) goto reduction;
         // reduces np to 1 and thus considerably reduces sz
            else{
               if (makeMovie) ppm2mp4(Word("video"));
               CPM_MZ
               return;
            }
         }
      }
      else{ // wf.np()<npMax
         R xU=wf.xU(); // upper end of biotope prior to extension
         R shf=xU-wgs.xL(); // wgs = wave function ground state
            // so shf and shiftRange may be rather different
         Wf wgsNew=wgs.shiftBio(shf);
         wf=wf.extU(wgsNew);
         cpmassert(wf.np()>1,
            "a new wave function was tensorially 'added'");
         N ws=wf.sz();
         if (ws>szMax){
            cpmmessage("szMax exceeded: szMax="&cpm(szMax)&
            ", wf.sz()="&cpm(ws));
            goto reduction;
         }
         R_Func potNew=pot0.shift(shiftRange)+wf.pot().shift(shiftRange);
         wf.setPot_(potNew);
      }
      cpmmessage("sys after extU etc.",1,true);
      cpmmessage("new Sys ahead");
      sys=Sys(wf,dtRel,t1);
      cpmmessage("new Sys done");
      gr.clearText();
      sys.mark(gr);
      gr.vis(makeMovie);
      cpmwait(tWait4,2);

      cpmassert(sys.freeL(),
         "why? action to ensure sys.freeL() was done before !");
      ostringstream ost2;
      ost2<<"if (!sys.freeU()) done"<<endl;
      cpmmessage(ost2,1,true);
   }
   if (sys.sz()<=szMax){
      ostringstream ost;
      ost<<"sys.sz()="<<sys.sz()<<", goto proceed"<<endl;
      cpmmessage(ost,4,true);
      goto proceed;
   }
// if we come to this point we have sys.sz>szMax and we have to reduce
// the wave function to that of a single particle state
   static Z redCount=0_Z;
reduction:
   ostringstream ostred;
   ostred<<"reduction started at sys.sz() = "<<sys.sz()<<endl;
   redCount++;
   cout<<" redCount = "<<redCount<<" sys.t() = "<<sys.t()<<endl;
   cpmmessage("label 'reduction' passed, redCount = "&cpm(redCount),1,true);
   cpmmessage(ostred,4,true);
   Wf wff=sys.x();
   R tff=sys.t();
   if (wff.np()>1){ // we need to create a one particle wave function
     C_Matrix densityOperator=wff.traceOut2(0); // we only look for the
      // projectile which is indexed 0
      Observable obs(densityOperator);
      SpectralObject spo=obs.spectralRepresentation();
      R_Vector sv=spo.getSpectrum();
      RandomSelector rs(sv);
      Z ri=rs();
      cout<<"ri="<<ri<<endl;
      State psi=spo[ri];
      gr.clearText();
      psi.mark(gr);
      gr.vis(true);
      cpmwait(tWait4,2);
      V<C> vc=psi.getRep();
      N n=vc.nc(); // nc: number of components
      cout<<"n="<<n<<endl;
      VC vc1(n);
      for (N i=0;i<n;++i) vc1[i]=vc.pi(i); // pi: plain index
      wff.simplify_(vc1);
      cpmmessage("new Sys ahead");
      sys=Sys(wff,dtRel,tff);
      cpmmessage("new Sys done");
   }
   cout<<"sys redefined after reduction, sys.sz() = "<<sys.sz()<<endl;
   gr.clearText();
   sys.mark(gr);
   gr.addText("sys after reduction");
   gr.vis(true);
   cpmwait(tWait4,2);
   if (sys.t()<tMax){
      sys.status(gr,tWait4);
      ostringstream ost3;
      ost3<<" t = "<<sys.t()<<", tMax = "<<tMax<<", still t < tMax"<<endl;
      ost3<<" goto proceed"<<endl;
      goto proceed;
   }
   if (makeMovie) ppm2mp4(Word("video"));
   CPM_MZ
   return;
}


void test3(){
   Z mL=1;
   Word loc("test3()");
   CPM_MA
   cout<<"test3() started"<<endl;
   Iv iv(-.5_R,.5_R);
   R_Func f=R_Func::tableMountain(iv,1e-12);
   Graph gr;
   gr.setXY(iv*1.1,Iv(-0.1,1.1));
   gr.mark(f);
   gr.vis();
   cpmwait(10.,2);
   cout<<"test3() done"<<endl;
   CPM_MZ
   return;
}

R f3(R const& x){return x*x*x;}

void test4(){
   Z mL=1;
   Word loc("test3()");
   CPM_MA
   cout<<"test4() started"<<endl;
   Viewport vp(600,300,"test4()");
   Frame fr;
   Graph gr(fr);
   gr.setX(Iv(-10.,10.));
   F<R,R> f(f3);
   gr.mark(f);
   gr.vis();
   Sys sys;
   cout<<"test4() done"<<endl;
   CPM_MZ
   return;
}

void test5(){
   Z mL=1;
   Word loc("test5()");
   CPM_MA
   cout<<"test5() started"<<endl;
   Sys sys;
   cout<<"test5() done"<<endl;
   CPM_MZ
}

void test6(){
   cout<<"test6() done"<<endl;
}

Z CpmApplication::main_(void)
// reading of variable sel from file eqm2.ini works only if
// program is called from the command line as (for sel=2)
// ~/e/cpm/eqm2 ./eqm2 2
// do not forget the './'
{
   Message::program="eqm2";
   cout<<"From smallExperiments/source/main_eqm2.cpp:"<<endl;
   cout<<"Z CpmApplication::main_(void) started"<<endl;

   Z mL=1;
   Z sel=-1; // no sel=1 !
   cpmverbose=1;
   R tWait=0.;
   bool all=false;

   Z ad=args_.dim();
   cout<<"args_ = "<<args_<<endl;

   if (ad>=2){
      cout<<"args_[2] = "<<args_[2]<<endl;
      sel=args_[2].toZ();
   }
   if (ad>=3){
      cout<<"args_[3] = "<<args_[3]<<endl;
      cpmverbose=args_[3].toZ();
   }
   if (ad>=4){
      cout<<"args_[4] = "<<args_[4]<<endl;
      tWait=args_[4].toR();
   }

   Word loc("CpmApplication::main_(), sel="&cpm(sel));
   CPM_MA
   cout<<"CpmApplication::main_(), sel= "<<sel<<endl;
   if (sel==2 || all) test2();
   if (sel==3 || all) test3();
   if (sel==4 || all) test4();
   if (sel==5 || all) test5();
   if (sel==6 || all) test6();
   cout<<"Z CpmApplication::main_(void) done, sel = "<<sel<<endl;
   CPM_MZ
   return sel;
}