//? classicalcrossway.cpp
//? C+- by Ulrich Mutze. Status of work 2022-12-27.
//? 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.


/*************************************************************************

This program tries to extend my Wolram Demonstration
      'Simulating Quantum Scattering'
/////////////////////////////////////////////////
   about
/////////////////////////////////////////////////
W file=quantum scattering.ini
Ws line1=Time evolution  of the crossway system.
Ws line2=...
Ws line3=...
Ws line4=...
Ws line5=...

These are the program specific data that go into the creation
of the start-up screen of the project. Other elements of
this start-up screen are common to all C+- programs which
use class CpmApplication::IniFileBasedApplication
(see file cpminifilebasapp.h) as an application framework.
The definition of function IniFileBasedApplication::announceTheWork()
in file cpminifilebasapp.cpp defines these additional elements.

Notice also the following section, which any ini-file is expected to
have and which allows to influence the program run and the
auto-generated documentation in many ways. Although
ini-files are not seen by the compiler, they read like
C++ code by having type indicators for their data entries.
Also comments are treated as in C++.
Valid type indicators are B, Z, R, W for bool, integer, real, and
character string and the corresponding lists Bs, Zs, Rs, Ws.
There is a further valid type indicator F, which introduces the name
of an ini-file to be included. So ini-files nest to an arbitrary level.
This feature was developed for an industrial project where all data
defining the kinematics of a complex robot were originally held in a
single ini-file.
With the number of data growing to several thousands and with
people from other groups becoming involved in adjusting the
data, it became mandatory to split the data to separatly
editable units and nevertheless ensure integrity of production runs
by having a single ini-file at work.
Several of my C+- programs make use of this facility for similar
reasons.

Text lines in which there is no data type indicator are section
headings. Sections behave very much like namespaces: Althoug
section 'about' (see above) has an entry 'file' a section with
different name may also have an entry named 'file'. Sections
don't nest, so that several F-included files may add to
the same section, which turned out to be very convenient in the
framework for which I developed the scheme.

////////////////////////////////////////////////////////////////////////
   selection
////////////////////////////////////////////////////////////////////////

R tWaitAnnounceScreen=0

Z sel=1
// a selector for the 'depth' of auto-generated documentation files.
// Control of auto-generated documentation files is based on the
// two entries sel, runName. These are the rules:
// 1. for sel<0 no documentation files will be created.
// 2. for sel==0 a light documentation file is created that
//   covers all the input so that the program run can be repeated.
// 3. for sel>0, in addition to the light documentation file
//    a more detailed one will be created that also can made to
//    list computational results.
// 4. All documentation files have names which contain runName
//    appended to the name of the calling program

Z cpmverbose=1
// controls the degree of detail with which program run
// information will be written to log file cpmcerr.txt

Z cpmdbg=2
// controls the reaction of the cpmassert macro:
// 1: warning on failing assertion
// 2: error stop on failing assertion

W runName=190901a

// identifier for the program run that will be appended to
// the names of documentation files.
// If the value of this quantity is "auto" this value will
// be replaced by an autogenerated name.

***********************************************************************/

#include <cpmbas1.h>
#include <cpmcmatrix.h>
#include <cpmgraph.h>
#include <cpmhistogr.h>
#include <cpmcrossway.h>

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 CpmQuantumMechanics;
using namespace std;
using namespace std::placeholders;
using namespace std::filesystem;

//namespace{ // anonymous namespace which hold all of classes
// ClassicalCrossWayApp and ClassicalCrossWay

///////////////////////////// class ClassicalCrossWay //////////////////////////


typedef CrossWay CW ;
typedef ClassicalCrossWay CCW ;
typedef ClassicalCrossWaySystem CCWS ;

void test1()
{
   Word loc("test1()");
   Z mL=1;
   CPM_MA
   RecordHandler rch("classicalcrossway.ini"); // all input data from file

   Word sec="biotope";
   R L1,L2;
   cpmrh(L1);
   cpmrh(L2);

   sec="initial conditions";
   R x10,sigma1;
   cpmrh(x10);
   cpmrh(sigma1);

   sec="system data";
   Z n1,n2,ni,nf,np;
   R m1,m2,widthRel,lambdaRel,dtRel,rho,eps;
   cpmrh(m1);
   cpmrh(m2);
   cpmrh(n1);
   cpmrh(n2);
   cpmrh(widthRel);
   cpmrh(lambdaRel);
   cpmrh(dtRel);
   cpmrh(rho);
   cpmrh(eps);
   cpmrh(ni);
   cpmrh(nf);
   cpmrh(np);

   sec="run data";
   Z nSteps,imgUpdate1,imgUpdate2;
   R gamma,tWait1,tWait2,tWait3,tWait4;
   Word movieName;
   cpmrh(nSteps);
   cpmrh(imgUpdate1);
   cpmrh(imgUpdate2);
   cpmrh(gamma);
   cpmrh(movieName);
   cpmrh(tWait1);
   cpmrh(tWait2);
   cpmrh(tWait3);
   cpmrh(tWait4);

   B makeMovie = movieName.dim()>0_Z ? B(1) : B(0);

   Color::setGamma(gamma);

   Z nSample = 5;
   Iv facRange(.9_R,1.1_R);
   Z jr=137_Z;
   for (Z is=1; is<=nSample; ++is){
      cpmdata<<"***************** is = "<<is<<endl;
      cout<<" is = "<<is<<endl;
      cpmvalue("is=",R(is),3);
      R facA = (is==1_Z ? 1_R : facRange.ran(jr++));
      R facB = (is==1_Z ? 1_R : facRange.ran(jr++));
      R facC = (is==1_Z ? 1_R : facRange.ran(jr++));
      cpmdata<<"facA ="<<facA<<endl;
      cpmdata<<"facB ="<<facB<<endl;
      cpmdata<<"facC ="<<facC<<endl;


      CrossWayData dat(L1,L2,m1,m2,n1,n2,x10,sigma1,
         widthRel,lambdaRel*facA,dtRel,rho*facB,eps*facC,ni,nf,np);
      cout<<"cw construction ahead"<<endl;
      CW cw=CrossWay::make(dat);
      cout<<"cw constructed"<<endl;
      R dt=cw.dt(dtRel);

//      Graph gr;
//      gr.setXY(cw.iv1(),cw.iv2());

      V<CCW> par=classical(dat); // np: number of trajectories ( not really 'particles'

      // as suggested by the 'p') to be created

      CyclicAlarm show1(imgUpdate1);
      for (Z i=1;i<=nSteps;++i){
         cw.step_(dt);
         if (show1()){
            cout<<" i = "<<i<<endl;
//            cw.mark(gr);
//            gr.vis(makeMovie);
         }
         for (Z j=1;j<=np;++j) par[j].step_(dt);
      }
      cpmwait(tWait1,2);

      R2 rCW=cw.pLpU();
      R2 rCCW=pLpU(par);
      cpmdata<<"rCW = "<<rCW<<endl;
      cpmdata<<"rCCW = "<<rCCW<<endl;
      R dis=rCW.dis(rCCW);
      cpmdata<<"dis = "<<dis<<endl;
      cout<<"dis = "<<dis<<endl;
   }
   CPM_MZ
   return;
}

void test2()
{
   Word loc("test2()");
   Z mL=1;
   CPM_MA

   RecordHandler rch("classicalcrossway.ini"); // all input data from file
   Word sec="system data";
   Z np,ns;
   cpmrh(np);
   cpmrh(ns);

   sec="run data";
   Z nSteps,imgUpdate1,imgUpdate2;
   R gamma,tWait1,tWait2,tWait3,tWait4;
   Word movieName;
   cpmrh(nSteps);
   cpmrh(imgUpdate1);
   cpmrh(imgUpdate2);
   cpmrh(gamma);
   cpmrh(movieName);
   cpmrh(tWait1);
   cpmrh(tWait2);
   cpmrh(tWait3);
   cpmrh(tWait4);
   Iv iv(.5_R,2_R);

   CrossWayData dat=CrossWayData::make("classicalcrossway.ini");
   V<CrossWayData> vrd=randomize(dat,iv,ns);
   for (Z i=1;i<=ns;++i){
      fpLpU(vrd[i],nSteps,np); // results on cpmdata
   }
   CPM_MA

   return;
}

///////////////////////// auxiliar class RZo ///////////////////////////////////

// Pairs of R and Z with order deremined by R alone.
class RZo{
   R x_{0_R};
   Z i_{0_Z};
   typedef RZo Type;
public:
   Word nameOf()const{return "RZo";}
   RZo(){}
   RZo(R x, Z i):x_(x),i_(i){}
   CPM_ORDER // essential
   CPM_IO // required for cooperation with S<>
   R x()const{return x_;}
   R i()const{return i_;}
};

Z RZo::com(RZo const& rz)const{ return Comp<R>()(x_,rz.x_);}
   // implements CPM_ORDER

// implements CPM_IO
bool RZo::prnOn(ostream& out)const
{
   if (!CpmRoot::writeTitle("RZo",out)) return false;
   if (!CpmRoot::write(x_,out)) return false;
   if (!CpmRoot::write(i_,out)) return false;
   return true;
}

bool RZo::scanFrom(istream& in)
{
   if (!CpmRoot::read(x_,in)) return false;
   if (!CpmRoot::read(i_,in)) return false;
   return true;
}

void test3()
{
   Word loc("test3()");
   Z mL=1;
   CPM_MA

   RecordHandler rch("classicalcrossway.ini"); // all input data from file

// Notice that in the present project we don't use the service from
// inifilebasapplication.cpp and so have to handle the input from 'selection'
// ourselves.
   Word sec="selection";
   Word runName;
   cpmrh(runName);
   Message::setRunId(runName);

   sec="run data";
   Z nSteps,imgUpdate1,imgUpdate2,imgUpdate3,frameRate;
   R gamma,tWait1,tWait2,tWait3,tWait4;
   Word movieName;
   cpmrh(nSteps);
   cpmrh(imgUpdate1);
   cpmrh(imgUpdate2);
   cpmrh(imgUpdate3);
   cpmrh(frameRate);
   cpmrh(gamma);
   cpmrh(movieName);
   cpmrh(tWait1);
   cpmrh(tWait2);
   cpmrh(tWait3);
   cpmrh(tWait4);

   sec="topic3";
   R xLimL,xLimU,mag1,mag2,mag3;
   Z style1,style2,style3,nTrail1,nTrail2,nTrail3,is{0_Z};
   B aspSensitive,wos1,wos2,wos3; // aspect ratio sensitive
   B noSingleParticleRuns;

   cpmrh(nTrail1);
   cpmrh(nTrail2);
   cpmrh(nTrail3);
   cpmrh(style1);
   cpmrh(style2);
   cpmrh(style3);
   cpmrh(mag1);
   cpmrh(mag2);
   cpmrh(mag3);
   cpmrh(aspSensitive);
   cpmrh(wos1);
   cpmrh(wos2);
   cpmrh(wos3);
   cpmrh(noSingleParticleRuns);
   cpmrh(is);

   R asp= aspSensitive ? 1_R*Viewport::pelX()/Viewport::pelY() : 0_R;
   CrossWayData dat=CrossWayData::make("classicalcrossway.ini",asp);
   ClassicalCrossWaySystem ccws(dat);

   R dt=dat.dt();
   R rho=dat.rho;
   R L1h=dat.L1*.5_R;
   Z np=dat.np;

   if (is > 0_Z ){ // we know on what particle we will concentrate
      ClassicalCrossWaySystem ccws1=ccws.reduce(is);
      ccws1.evolve_(nSteps,style2,mag2,imgUpdate2,wos2,nTrail2);
      if (wos2) ppm2mp4(movieName,frameRate);
      return;
   }

   ClassicalCrossWaySystem ccws0=ccws;
   ccws.evolve_(nSteps,style1,mag1,imgUpdate1,wos1,nTrail1);
   cpmwait(tWait1,2);
   V<R> x(np);
   for (Z i=1; i<=np;++i) x[i]=ccws[i].x1();
   S<RZo> szL,szC,szU;
   for (Z i=1; i<=np;++i){
      R xi=x[i];
      if (xi<-rho){
         szL.add_(RZo(xi,i));
      }
      else if (xi>rho){
         szU.add_(RZo(xi,i));
      }
      else{
         szC.add_(RZo(xi,i));
      }
   }
   R npi=1_R/np;
   Z nL=szL.car(),nC=szC.car(),nU=szU.car();
   R pL=nL*npi,pC=nC*npi,pU=nU*npi;
   cpmdata<<"pL = "<<pL<<endl;
   cpmdata<<"pC = "<<pC<<endl;
   cpmdata<<"pU = "<<pU<<endl;

   Vo<RZo> vzL=szL.toV();
   Vo<RZo> vzC=szC.toV();
   Vo<RZo> vzU=szU.toV();

   vzL.sort_();
   vzC.sort_();
   vzU.sort_();

   Z ill=vzL.fir().i();
   Z il=vzL.last().i();
   Z ic=vzC.median().i();
   Z iu=vzU.fir().i();
   Z iuu=vzU.last().i();
   V<Z> vi{ill, il, ic, iu, iuu};
   cpmdata<<"vi = "<<vi<<endl;
   cpmdata<<"x[ill] = "<<x[ill]<<endl;
   cpmdata<<"x[il] = "<<x[il]<<endl;
   cpmdata<<"x[ic] = "<<x[ic]<<endl;
   cpmdata<<"x[iu] = "<<x[iu]<<endl;
   cpmdata<<"x[iuu] = "<<x[iuu]<<endl;
   if (wos1) ppm2mp4(movieName,frameRate);

   if (noSingleParticleRuns) return;

   for (Z i=1;i<=5;++i){
      Z is=vi[i];
      cpmvalue("is",is,3);
      ClassicalCrossWaySystem ccws1=ccws0.reduce(is);
      ccws1.evolve_(nSteps,style2,mag2,imgUpdate2,wos2,nTrail2);
   }
   cpmwait(tWait2,2);
   for (Z i=1;i<=5;++i){
      Z is=vi[i];
      cpmvalue("is",is,3);
      ClassicalCrossWaySystem ccws1=ccws0.reduce(is);
      ccws1.evolve_(nSteps,style3,mag3,imgUpdate3,wos3,nTrail3);
   }
   cpmwait(tWait3,2);
   if (wos1||wos2||wos3) ppm2mp4(movieName,frameRate);
   CPM_MA
}


Z CpmApplication::main_(void)
// In this form the C++ typic function main()(defined in file
// cpmapplication.cpp) expects to see the basic functionality of
// the program.
// Notice that main() in cpmapplication.cpp provides the program
// with its graphical window and a status bar.
// The fact that class QuantumScatteringApp is derived from class
// IniFileBasedApplication lets the following short code, although it
// looks like abstract nonsense, do very specific things.
{
   Word loc("Z CpmApplication::main_(void)");
   Z mL=1;
   CPM_MA
   RecordHandler rch("classicalcrossway.ini");
   Word sec="select topic";
   Z topic;
   cpmrh(topic);

   if (topic==1){
      test1();
   }
   else if (topic==2){
      test2();
   }
   else if (topic==3){
      test3();
   }
   else {
      cout<< "unvalid value of msel nothing to do"<<endl;
   }
   cpmwait(4);
   CPM_MZ
   return 0;
}