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


#include <cpmgraph.h>
#include <cpmimagingtools.h>
#include <cpmsystemdependencies.h>
#include <cpmfile.h>

#include <vector>
#include <iomanip>

using CpmGraphics::Viewport;

using CpmGraphics::Frame;
using CpmGraphics::Graph;
using CpmGraphics::Fig;
//using CpmGraphics::gpScript;

using CpmLinAlg::R2;
using CpmLinAlg::Z2;
using CpmLinAlg::Z3;
using CpmGraphics::Rn2;
using std::vector;

using namespace CpmRoot;
using namespace CpmSystem;
using namespace CpmArrays;
using namespace CpmFunctions;
using namespace CpmImaging;

//////////////////////// function gpScript ///////////////////////////////

void CpmGraphics::gpScript(Word const& fn, Word const& title,
   Iv const& xRange, Iv const& yRange,
   Word const& xLabel, Word const& yLabel,
   Word const& output, Word const& datFile,
   V<Word> const& leg, V<Word> const& col,
   Word const& hPos, Word const& vPos)
{
   Z mL=3;
   Word loc("CpmGraphics::gpScript(Word,Word,Iv,Iv ... )");
   CPM_MA
   OFile res(fn);
   Z eowMem=Word::getEndOfWord();
   Z i;
   Z n=leg.dim();
   Word::setEndOfWord(0);

   res()<<"set style data lines"<<endl;
   for (i=1;i<=n;++i){
      res()<<"set linetype "<<i<<" lc rgb '"<<col.cyc(i)<<"'"<<endl;
   }
   if ((hPos=="left"||hPos=="center"||hPos=="right")
       &&(vPos=="bottom"||vPos=="center"||vPos=="top")){
          res()<<"set key "<<hPos<<" "<<vPos<<endl;
   }
   res()<<"set title '"<<title<<"'"<<endl;
   res()<<"set terminal pdf color enhanced"<<endl;
   res()<<"set size 1,1"<<endl;
   res()<<"set xrange [ "<<xRange[1]<<" : "<<xRange[2]<<" ]"<<endl;
   if (yRange[1]!=yRange[2]){
       res()<<"set yrange [ "<<yRange[1]<<" : "<<yRange[2]<<" ]"<<endl;
   }
   res()<<"set xlabel '"<<xLabel<<"'"<<endl;
   res()<<"set ylabel '"<<yLabel<<"'"<<endl;
   res()<<"set output '"<<output<<"'"<<endl;
   res()<<"plot \\"<<endl;
   for (i=1;i<=n-1;++i){
      res()<<"'"<<datFile<<"' using 1 : "<<i+1<<" title '"<<leg[i]<<"'"
      <<", \\"<<endl;
   }
   res()<<"'"<<datFile<<"' using 1 : "<<n+1<<" title '"<<leg[n]<<"'"<<endl;

   res()<<"set output"<<endl;
   Word::setEndOfWord(eowMem);
   CPM_MZ
}

///////////////////////////// class Graph ///////////////////////////////

namespace{
// used as initialization in iniGraph(...) , data may be changed by
// member functions.
   const R GRID_XPOS=0.05;
   const R GRID_YPOS=0.92;
   const R GRID_HEIGHT=15;
// some badly treated method parameters
   const Z N_EXCESS=100;
   // 2001-3-21 index out of range for N_EXCESS=40 obtained
   // now no longer needed in function polygonClip
   const R IMMENSE=1e99;
}

// basic new geometry functions

namespace{
   const R tinyRim=1e-6;
//   const R tinyRim=1e-3;
}

Z   Graph::to_i(R x, Z& excess)const
{
   excess=0;
   R xRel=(x-CL.x1)*dXInv;
   if (xRel<=0){
      if (xRel<-tinyRim) excess=-1;
      return 0;
   }
   else if (xRel>=1){
      if (xRel>tinyRim+1) excess=1;
      return nCol-1;
   }
   else return cpmtoz(xRel*nCol);
}

Z  Graph::to_j(R y, Z& excess)const
{
   excess=0;
   R yRel_=(CU.x2-y)*dYInv; // reversal of direction
      // as compared with to_x
   if (yRel_<=0){
      if (yRel_<-tinyRim) excess=+1; // sic, y too large!
      return 0;
   }
   else if (yRel_>=1){
      if (yRel_>tinyRim+1) excess=-1;
      return nRow-1;
   }
   else return cpmtoz(yRel_*nRow);
}

////////// end geometry

void Graph::iniGraph(R Xl, R Yl, R dX, R dY)
{
   setXY( Xl, Yl, dX, dY);
   gridColor=GRIDCOLOR;
   gridTextColor=WRITECOLOR;
   graphColor=RED;
   withOrigin=1;
   newGrid=1;
   autoScale=B(true);
   gridText=V<Word>(1);
   gridTextPosition.x1=GRID_XPOS;
   gridTextPosition.x2=GRID_YPOS;
   gridTextHeight_=GRID_HEIGHT;
}

Graph::Graph(R xl, R yl, R dx, R dy, const Color& color, R Xl, R Yl,
   R dX, R dY):Frame(xl,yl,dx,dy,color)
{
   iniGraph(Xl, Yl, dX, dY);
}

Graph::Graph(void):Frame()
{
   iniGraph(0,0,1,1);
}

Graph::Graph(const Color& c):Frame(c)
{
   iniGraph(0,0,1,1);
}

Graph::Graph(const Frame& f):
Frame(f)
{
   iniGraph(0,0,1,1);
}

Graph::Graph(Frame const& fr, R rR):
Frame(rR,rR*fr.getAspRat(),R(1)-rR*2,R(1)-rR*fr.getAspRat()*2)
{
   iniGraph(0,0,1,1);
}
bool Graph::prnOn( ostream& str) const
// return preliminary
//notice that content is not copied and retrieved !!!!
{
   using CpmRoot::write;
   Z mL=2;
   Word loc("Graph::prnOn(ostream)");
   CPM_MA
   cpmwt("Graph");
   Frame::prnOn(str);
   cpmwt("CL");
   CL.prnOn(str);
   cpmwt("CU");
   CU.prnOn(str);
   CPM_MZ
   return true;
}

bool Graph::scanFrom (istream& in)
// return preliminary
{
   using CpmRoot::read;
   Z mL=2;
   Word loc("Graph::scanFrom (istream)");
   CPM_MA
   Frame::scanFrom(in);
   CL.scanFrom(in);
   CU.scanFrom(in);
   updatedXdY();
   CPM_MZ
   return true;
}

void Graph::setGridColor(const Color& c)
{
   gridColor=c;
}

void Graph::setGridTextColor(const Color& c)
{
   gridTextColor=c;
}

void Graph::setGraphColor(const Color& c)
{
   graphColor=c;
}

void Graph::setWithOrigin(Z c)
{
   withOrigin=c;
}

void Graph::setNewGrid(Z c)
{
   newGrid=c;
}

void Graph::setGridText(Word a)
{
   gridText[1]=a;
}

void Graph::setGridTextPosition(const R2& pos)
{
   gridTextPosition=pos;
}

void Graph::setGridTextPosition(R xText, R yText)
// pos is understood in the Frame coordinate system (xF,yF)
{
   gridTextPosition=R2(xText,yText);
}

Z  Graph::iX(R X)const
{
   Z exc;
   return to_i(X,exc)+a0();
}

R  Graph::iX_(R x)const
{
   return to_i_(x)+a0();
}

Z  Graph::jY(R Y)const
// Y<=CL.x2 corresponds to the lower boundery of the Frame
// Y>=CU.x2 corresponds to the upper boundery of the Frame
{
   Z exc;
   return to_j(Y,exc)+b0();
}

R  Graph::jY_(R y)const
{
   return to_j_(y)+b0();
}

Z   Graph::iX0(R x)const
// x=XL corresponds to left boundery of the Frame
// x=XU corr. to right bound. of Frame
// for x<XL or x>XU we return -1
{
   if (x<CL.x1 || x> CU.x1) return -1;
   else return iX(x);
}

Z   Graph::jY0(R y)const
// y=YL corresponds to the lower boundery of the Frame
// y=YU corresponds to the upper boundery of the Frame
// for y<YL or y>YU we return -1
{
   if (y<CL.x2 || y> CU.x2) return -1;
   else return jY(y);
}

Z2 Graph::toPixel(const R2& p)const
// maps a point from object space into a pixel (discrete position on
// screen
{
   Z2 res(iX(p.x1),jY(p.x2));
   return res;
}

R2 Graph::toObjectSpace(const Z2& z )const
{
   Word loc("Graph::toObjectSpace(Z2)");
   cpmassert(pU.x1!=pL.x1,loc);
   cpmassert(pU.x2!=pL.x2,loc);
   R dx=(CU.x1-CL.x1)/(pU.x1-pL.x1);
   R dy=(CU.x2-CL.x2)/(pL.x2-pU.x2);   // positive with pL-pU
   R2 p;
   p.x1=CL.x1+dx*(z.x1+0.5-pL.x1);
   p.x2=CU.x2-dy*(z.x2+0.5-pU.x2);
    return p;
}

R2 Graph::get_xFyF(const R2& xGyG)const
{
   R xG=xGyG.x1, yG=xGyG.x2;
   R xF=(xG-CL.x1)*dXInv, yF=(yG-CL.x2)*dYInv;
   return R2(xF,yF);
}

R2 Graph::get_xGyG(const R2& xFyF)const
{
   R xF=xFyF.x1, yF=xFyF.x2;
   R xG=CL.x1+xF*dX, yG=CL.x2+yF*dY;
   return R2(xG,yG);
}

bool Graph::updatedXdY(void)
{
   if (!autoScaleY){
      dX=CU.x1-CL.x1;
      dXInv=cpminv(dX);
      if (dX<=0){
         cpmwarning("Graph::updatedXdY(): dX<=0: dX="&cpmwrite(dX));
         return false;
      }
   }
   dY=CU.x2-CL.x2;
   dYInv=cpminv(dY);
  // cout<<"Graph::updatedXdY: dX="<<dX<<" dY="<<dY<<endl;
   //cpmcheck(dYInv);
   if (dY<=0){
      cpmwarning("Graph::updatedXdY(): dY<=0: dY="&cpmwrite(dY));
      return false;
   }
   return true;
}

bool Graph::setXY(R X, R Y, R dX1, R dY1)
{
   R X_=X+dX1, Y_=Y+dY1;
   CpmRootX::order(X,X_);
   CpmRootX::order(Y,Y_);
   CL.x1=X; CU.x1=X_;
   CL.x2=Y; CU.x2=Y_;
   return updatedXdY();
}

bool Graph::setX(Iv iv, R fac)
{
   if (fac!=1_R) iv.stretch_(fac);
   CL.x1=iv.inf(); CU.x1=iv.sup();
   autoScale=B(false);
   return updatedXdY();
}

bool Graph::setY(Iv iv, R fac )
{
   if (fac!=1_R) iv.stretch_(fac);
   CL.x2=iv.inf(); CU.x2=iv.sup();
   autoScale=B(false);
   return updatedXdY();
}

bool Graph::setXY(Iv const& ivx, Iv const& ivy, R fac)
{
   Iv iv1=ivx;
   Iv iv2=ivy;
   iv1.stretch_(fac);
   iv2.stretch_(fac);
   CL.x1=iv1.inf(); CU.x1=iv1.sup();
   CL.x2=iv2.inf(); CU.x2=iv2.sup();
   autoScale=B(false);
   return updatedXdY();
}

// Setting the window according to the size of objects to be
// represented graphically. Overloading with differently
// typed arguments may be added later according to need

namespace{

void splitVal(R& a, R& b, R tol) // changed 2007-10-03
{
// precondition: a<=b and tol>0;  postcondition: a<b
   if (a>b) cpmerror("splitVal(R&,R&,R), a > b not allowed");
   if (a==b){
      R aa=cpmabs(a);
      R eps= aa==0 ? tol : aa*tol;
      a-=eps;
      b+=eps;
   }
}

} // namespace

bool Graph::setXY(const R_Vector& x, const R_Matrix& y)
{
   const Z mL=4;
   Word loc("Graph::setXY(R_Vector,R_Matrix)");
   CPM_MA

   const R fracTol=1e-6; // the object space should be a bit larger since
      // a horicontal line (x[i],y[k][i]), i=1,2,... can not define
      // exacly the boundary of the object space. Otherwise numerical
      // noise would influence whether this line gets drawn in
      // showGL(R_Vector,R_Matrix,...).
   Z n=x.dim(),  m=y.dim1(), n1=y.dim2(), i;
   if (n!=n1) cpmerror("Graph::setXY(R_Vector,R_Matrix): n!=n1");
   cpmassert(n>0,loc);

   Iv ivx=x.infSup();
   R x_min=ivx.inf(), x_max=ivx.sup();
   splitVal(x_min,x_max,fracTol);
   R d_x=x_max-x_min;

   R_Vector yMin(m);
   R_Vector yMax(m);
   for (i=1;i<=m;i++){
      Iv ivi=y[i].infSup();
      yMin[i]=ivi.inf();
      yMax[i]=ivi.sup();
   }
   R y_min=yMin.inf();
   R y_max=yMax.sup();
   if (!isVal(y_max)) cpmcerr<<yMax;
   splitVal(y_min,y_max,fracTol);
   if (withOrigin==1){ // This allows make the point y=0 appear in
               // the graph
      if (y_min>0.){
         y_min=0. ;
      }
      else if (y_max<0.){
         y_max=0. ;
      }
   }
   R d_y=y_max-y_min; // moved from *** to here, 2007-10-03
   R eps_x=fracTol*d_x;
   CL.x1=x_min-eps_x;
   CU.x1=x_max+eps_x;
   R eps_y=fracTol*d_y;
   CL.x2=y_min-eps_y;
   CU.x2=y_max+eps_y;

   bool res=updatedXdY();
   CPM_MZ
   return res;
}


bool Graph::setXY(const R_Vector& x, const R_Vector& y, Z method)
{
   const R secFac=1.01;
      // a 'security factor' which lets the x,y data being displayed
      // without ending up on the rim of the Graph
   if (method==0){
      R_Matrix yM(1);
      yM[1]=y;
      return setXY(x,yM);
   }
   else{ // so far, maybe other methods may have to be added
      Z n=x.dim(),i;
      cpmassert(n==y.dim(),"Graph::setXY(R_Vector,R_Vector,Z)");
      R x1=x[1];
      R y1=y[1];

      R_Vector xc=x;
      R_Vector yc=y;
      for (i=1;i<=n;++i){
         xc[i]-=x1;
         yc[i]-=y1;
         xc[i]=cpmabs(xc[i]);
         yc[i]=cpmabs(yc[i]);
      }
      R xMax=xc.sup();
      R yMax=yc.sup();
      xMax*=secFac;
      yMax*=secFac;
      R asp=getAspRat();
      R xComp=xMax/asp;
      R yComp=yMax;

      Iv ivx,ivy;
      if (yComp>=xComp){
         ivy=Iv(y1-yMax,y1+yMax);
         ivx=Iv(x1-yMax*asp,x1+yMax*asp);
      }
      else{
         ivx=Iv(x1-xMax,x1+xMax);
         ivy=Iv(y1-xComp,y1+xComp);
      }
      return Graph::setXY(ivx,ivy);
   }
}

bool Graph::setXY(const R_Matrix& y)
{
   R_Vector x=y[1];
   Z m=y.dim1()-1;
   if (m==1){ // y has only two lines and represents one curve
      R_Vector yy;
      yy=y[2];
      return setXY(x,yy);
   }
   else{   // y has more than two lines and represents a family
      // of curves
      R_Matrix yy(m);
      for (Z i=1;i<=m;i++) yy[i]=y[i+1];
      return setXY(x,yy);
   }
}

bool Graph::setXY(const V<R2>& pos, R enlarge)
{
   Z n=pos.dim();
   Vo<R> posx(n);
   Vo<R> posy(n);
   for (Z i=1;i<=n;i++){
      posx[i]=pos[i].x1;
      posy[i]=pos[i].x2;
   }
   T2<Z> limx=posx.indInfSup();
   T2<Z> limy=posy.indInfSup();
   Iv ivx(posx[limx.c1()],posx[limx.c2()]);
   Iv ivy(posy[limy.c1()],posy[limy.c2()]);
   ivx.stretch(enlarge); // experiment
   ivy.stretch(enlarge);
   return Graph::setXY(ivx,ivy);
}

bool Graph::setXY(V< V<R2> > const& pos, R enlargeX, R enlargeY)
// implementation without ordering and with fast iteration
// till 2007-11-24 I had a (if, else if, else)-construct
// which failed to work in rare cases. Now everything is waterproof.
{
 //  R rHuge=R1::max;
   R rHuge=1.e64;

   R xMin=rHuge;
   R yMin=rHuge;
   R xMax=-rHuge;
   R yMax=-rHuge;

   for (Z i=pos.b();i<=pos.e();i++){
     for (Z j=pos[i].b();j<=pos[i].e();++j){
         R2 xy=pos[i][j];
         R x=xy.x1;
         if (x<xMin) xMin=x;
         if (x>xMax) xMax=x;
         R y=xy.x2;
         if (y<yMin) yMin=y;
         if (y>yMax) yMax=y;
     }
   }
   Iv ivx(xMin,xMax);
   Iv ivy(yMin,yMax);
   ivx.stretch(enlargeX);
   ivy.stretch(enlargeY);
   return Graph::setXY(ivx,ivy);
}

bool Graph::setXY(const V<C>& pos, R enlarge, Z met)
{
   Z mL=3;
   Word loc("Graph::setXY(const V<C>& pos, R enlarge)");
   CPM_MA
   Z n=pos.dim();
   Vo<R> posx(n);
   Vo<R> posy(n);
   for (Z i=1;i<=n;i++){
      posx[i]=pos[i][1];
      posy[i]=pos[i][2];
   }
   T2<Z> limx=posx.indInfSup();
   T2<Z> limy=posy.indInfSup();
   Iv ivx(posx[limx.c1()],posx[limx.c2()]);
   Iv ivy(posy[limy.c1()],posy[limy.c2()]);
// avoiding zero sized window 2007-10-04
   const R tol=1e-6;
   R a=ivx[1], b=ivx[2];
   splitVal(a, b, tol);
   ivx=Iv(a,b);
   a=ivy[1]; b=ivy[2];
   splitVal(a, b, tol);
   ivy=Iv(a,b);
   bool res;
   if (met==1){
      Iv ivAll=ivx|ivy;
      ivAll.stretch(enlarge);
      res=Graph::setXY(ivAll,ivAll);
   }
   else{
      ivx.stretch(enlarge);
      ivy.stretch(enlarge);
      res=Graph::setXY(ivx,ivy);
   }
   CPM_MZ
   return res;
}

void Graph::adjustScaleX(void)
{
   R asp=screenAspRat(); // added 2005-08-23
   R scale=(CU.x2-CL.x2)/(cU.x2-cL.x2);
   R dX=scale*(cU.x1-cL.x1)*asp;
   CU.x1=CL.x1+dX;
   updatedXdY();
}

void Graph::adjustScaleY(void)
{
   R aspInv=screenAspRatInv();
   R scale=aspInv*(CU.x1-CL.x1)/(cU.x1-cL.x1);
   R dY=scale*(cU.x2-cL.x2);
   CU.x2=CL.x2+dY;
   updatedXdY();
}

bool Graph::scale(R magnification)
{
   R2 C0=0.5*(CU+CL);  // center
   R2 dw=CU-C0;
   dw*=magnification;
   CU=C0+dw;
   CL=C0-dw;
   return updatedXdY();
}

bool Graph::scaleY(R mag)
{
   Iv iy=getIvY();
   iy=iy.stretch(mag);
   return setY(iy);
}

void Graph::shift(R sx, R sy)
{
   R2 ds(sx*dX,sy*dY);
   CU+=ds;
   CL+=ds;
   updatedXdY();
}

void Graph::shiftAbs(R sx, R sy)
{
   R2 ds(sx,sy);
   CU+=ds;
   CL+=ds;
   updatedXdY();
}

void Graph::getObjectSpaceFrom(const Graph& gr)
{
   CL.x1=(gr.CL).x1;
   CL.x2=(gr.CL).x2;
   CU.x1=(gr.CU).x1;
   CU.x2=(gr.CU).x2;
   updatedXdY();
}

Z  Graph::adjust(R2& x)const
{
   R tinyx=cpmabs(CU.x1-CL.x1)*1e-10;  // 93-3-10 static R changed to R
   R tinyy=cpmabs(CU.x2-CL.x2)*1e-10;  // static leaves tinyx,tinyy the
   // same even if Xi,Yi where changed
   R x1=x.x1,x2=x.x2;
   Z  k=0;
   if (cpmabs(x1-CL.x1)<tinyx) { x.x1=CL.x1; k++;}
   if (cpmabs(x1-CU.x1)<tinyx) { x.x1=CU.x1; k++;}
   if (cpmabs(x2-CL.x2)<tinyy) { x.x2=CL.x2; k++;}
   if (cpmabs(x2-CU.x2)<tinyy) { x.x2=CU.x2; k++;}
   return k;
}

Z  Graph::boundary(const R2& p)const
{
   R2 u=p; adjust(u);    // otherwise not meaningfull
   R x_=u.x1,y_=u.x2;
   if (CL.x1<x_ && x_<CU.x1 && CL.x2<y_ && y_<CU.x2) return 0;
   if (y_==CL.x2 && CL.x1<=x_ && x_<CU.x1) return 1;
   if (x_==CU.x1 && CL.x2<=y_ && y_<CU.x2) return 2;
   if (y_==CU.x2 && CL.x1<x_ && x_<=CU.x1) return 3;
   if (x_==CL.x1 && CL.x2<y_ && y_<=CU.x2) return 4;
   return 5;
}

Z  Graph::same_boundary(R2 p, R2 q)
{
   adjust(p); adjust(q);      // no effect if the points are exactly
                        // on the boundary
   Z  j1=boundary(p);
   Z  j2=boundary(q);
   if (j1==0||j1==5||j2==0||j2==5) return 0;
   if (j1==j2) return 1; // fastest decision but not complete
   if (p.x1==CL.x1 && q.x1==CL.x1) return 1;
   if (p.x1==CU.x1 && q.x1==CU.x1) return 1;
   if (p.x2==CL.x2 && q.x2==CL.x2) return 1;
   if (p.x2==CU.x2 && q.x2==CU.x2) return 1; // this would be complete
   // also without the fast pre-decision
   return -1;
}

Z  Graph::next_corner(R2 b, R2& c)
{
   adjust(b); // no effect if b is precisely on boundary
   R x=b.x1,y=b.x2;
   if ((y==CL.x2) && CL.x1<=x && x<CU.x1)
      { c.x1=CU.x1; c.x2=CL.x2; return 2;}
   if ((x==CU.x1) && CL.x2<=y && y<CU.x2)
      { c.x1=CU.x1; c.x2=CU.x2; return 3;}
   if ((y==CU.x2) && CL.x1<x && x<=CU.x1)
      { c.x1=CL.x1; c.x2=CU.x2; return 4;}
   if ((x==CL.x1) && CL.x2<y && y<=CU.x2)
      { c.x1=CL.x1; c.x2=CL.x2; return 1;}
   c.x1=CU.x1; c.x2=CU.x2; return 0;
}

Z  Graph::line_bound_inters(R a, R b, R c, Rn2& y)
{
   Rn2 p(4), x(2);
   R2 q;
   if (a==0) {
      if (b==0) return 0;
      R y_=-c/b;
      if (CL.x2<=y_ && y_<=CU.x2){
         x[0].x1=CL.x1;
         x[0].x2=y_;
         x[1].x1=CU.x1;
         x[1].x2=y_;
         y=x;
         return 2;
      }
      else return 0;
   }
   if (b==0) {
      if (a==0) return 0;
      R x_=-c/a;
      if (CL.x1<=x_ && x_<=CU.x1) {
         x[0].x1=x_;
         x[0].x2=CL.x2;
         x[1].x1=x_;
         x[1].x2=CU.x2;
         y=x;
         return 2;
      }
      else return 0;
   }
   R b_=R(-1)/b, a_=R(-1)/a;
   p[0].x2=b_*(c+a*CL.x1); //y-value belonging to x=CL.x1
   p[0].x1=CL.x1;
   p[1].x2=b_*(c+a*CU.x1);  // y-value for CU.x1
   p[1].x1=CU.x1;
   p[2].x1=a_*(c+b*CL.x2);  //x-value for CL.x2
   p[2].x2=CL.x2;
   p[3].x1=a_*(c+b*CU.x2);  // x-value for CU.x2
   p[3].x2=CU.x2;

   Z  k=0;
   for (Z  i=0; i<4; i++) { q=p[i];
      adjust(q);
      Z  j=boundary(q);
      if (j>0 && j<5) {
         if (k<2){
            x[k]=q; k++;
         }
         else{
            return k;
         }
      }
   }
    y=x;
   return k;
}

Z  Graph::connect_on_boundary(R2 xi, R2 xf, Rn2& yc)
{
   R2 xa,p;
   Rn2 xc(6);
   adjust(xi); adjust(xf);
   Z  k1=boundary(xi);
   Z  k2=boundary(xf);
   if (k1==0 || k1==5)
      cpmwarning("point not on boundary in connect_on_b..");
   if (k2==0 || k2==5)
      cpmwarning("point not on boundary in connect_on_b..");
    xc[0]=xa=xi;
   Z  k=1;
iteration:
   Z  jj=same_boundary(xa,xf);
   if (jj==1) {
      xc[k]=xf;
      yc=xc;
      return k;
   }
   else { next_corner(xa,p);
         xc[k]=xa=p;
         k++;
         if (k>5) cpmerror("Graph","connect_on_boundary: k too large");
         goto iteration;
   }
}

R Graph::goto_boundary(R2 x, R2 e, R2& y)
{
   Rn d(4),d_pos(4);
   R2 xa;
   Z  i;
   xa=x;
   adjust(xa);
   if (e.x1==0. && e.x2==0.){
      y=x; return IMMENSE;
   }
   R r0=(e.x1==0.? IMMENSE : R(1)/e.x1);// one of those may be zero
   R r1=(e.x2==0.? IMMENSE : R(1)/e.x2);
   d[0]=(CL.x1-xa.x1)*r0;
   d[1]=(CU.x1-xa.x1)*r0;
   d[2]=(CL.x2-xa.x2)*r1;
   d[3]=(CU.x2-xa.x2)*r1;
   Z  k=0;
   // The points xa+d[i]*e belong to one of the defining infinite lines
   for (i=0;i<4;i++){
      if (d[i]>0){
         d_pos[k]=d[i] ;
         k++;
      }
   }
   //Now k is the number of positiv d's found.
   if (k==0) { y.x1=x.x1; y.x2=x.x2; return IMMENSE;}
   // this means that  in the direction of e no boundary point was found.
   // Maybe, in direction -e there may be one,
   // but this is of no interest.
   R d_=IMMENSE;
   for (i=0;i<k;i++){
      y=x+d_pos[i]*e;
      adjust(y);
      Z  z=boundary(y);
      if (z>0 && z<5 && (d_pos[i]<d_)){
         d_=d_pos[i];
      }
   }
   if (d_ == IMMENSE){
      y=x;
      return d_;
   }
   else {
      y=x+d_*e;
      adjust(y);
      return d_;
   }
}

bool Graph::indFunc(const R2& p)const
{
   if (p.x1>CU.x1) return false;
   if (p.x1<CL.x1) return false;
   if (p.x2>CU.x2) return false;
   if (p.x2<CL.x2) return false;
   return true;
}

void Graph::mark( const R2& p, const Color& color)const
{
   Z exc;
   Z i=to_i(p.x1,exc);
   if (exc) return;
   Z j=to_j(p.x2,exc);
   if (exc) return;
   img.putPel(i,j,color);
}

void Graph::mark(const R2& p, Z size, const Color& color)const
{
   Z exc;
   Z i=to_i(p.x1,exc);
   if (exc) return;
   Z j=to_j(p.x2,exc);
   if (exc) return;
   img.putPel(i,j,color);

   if (size==1) return;
   i++;
   img.putPel(i,j,color);

   if (size==2) return;
   j++;
   img.putPel(i,j,color);

   if (size==3) return;
   i--;
   img.putPel(i,j,color);

   if (size==4) return;
   i--;
   img.putPel(i,j,color);

   if (size==5) return;
   j--;
   img.putPel(i,j,color);

   if (size==6) return;
   j--;
   img.putPel(i,j,color);

   if (size==7) return;
   i++;
   img.putPel(i,j,color);

   if (size==8) return;
   i++;
   img.putPel(i,j,color);

   if (size==9) return;
   i++;
   img.putPel(i,j,color);

   if (size==10) return;
   j++;
   img.putPel(i,j,color);

   if (size==11) return;
   j++;
   img.putPel(i,j,color);

   if (size==12) return;
   j++;
   img.putPel(i,j,color);

   if (size==13) return;
   i--;
   img.putPel(i,j,color);

   if (size==14) return;
   i--;
   img.putPel(i,j,color);

   if (size==15) return;
   i--;
   img.putPel(i,j,color);

   if (size==16) return;
   i--;
   img.putPel(i,j,color);

    if (size==17) return;
   j--;
   img.putPel(i,j,color);

    if (size==18) return;
   j--;
   img.putPel(i,j,color);

    if (size==19) return;
   j--;
   img.putPel(i,j,color);

    if (size==20) return;
   j--;
   img.putPel(i,j,color);

    if (size==21) return;
   i++;
   img.putPel(i,j,color);

    if (size==22) return;
   i++;
   img.putPel(i,j,color);

    if (size==23) return;
   i++;
   img.putPel(i,j,color);

    if (size==24) return;
   i++;
   img.putPel(i,j,color);

    if (size==25) return;
   i++;
   img.putPel(i,j,color);

    if (size==26) return;
   j++;
   img.putPel(i,j,color);

    if (size==27) return;
   j++;
   img.putPel(i,j,color);

    if (size==28) return;
   j++;
   img.putPel(i,j,color);

    if (size==29) return;
   j++;
   img.putPel(i,j,color);

    return;
}

void Graph::draw( R2 const& p1,  R2 const& p2,
   Color const& c)
{
   R x1=iX_(p1.x1);
   R y1=jY_(p1.x2);
   R x2=iX_(p2.x1);
   R y2=jY_(p2.x2);
   line(x1,y1,x2,y2,c); // def in Frame
}

void Graph::fillCircle(const R2& center, R radius, const Color& c)
{
   R  x1=iX_(center.x1-radius);
   R  y1=jY_(center.x2-radius);
   R  x2=iX_(center.x1+radius);
   R  y2=jY_(center.x2+radius);
   Frame::fillEllipse(x1,y1,x2,y2,c);
}

void Graph::fillTriangle(R2 const& p1, R2 const& p2, R2 const& p3,
   Color const& color)
{
   R2 q1=natTrn(p1);
   R2 q2=natTrn(p2);
   R2 q3=natTrn(p3);
   Frame::fill(TriAng(q1,q2,q3),color);
}

void Graph::fillRectangle(R2 const& p1, R2 const& p2,
   Color const& color)
{
   R2 q1=natTrn(p1);
   R2 q2=natTrn(p2);
   Frame::fillRectangle(q1.x1,q1.x2,q2.x1,q2.x2,color);
}

void Graph::fillEllipse(R2 const& p1, R2 const& p2,
   Color const& color)
{
   R2 q1=natTrn(p1);
   R2 q2=natTrn(p2);
   Frame::fillEllipse(q1.x1,q1.x2,q2.x1,q2.x2,color);
}

void Graph::fillDisk(R2 const& pc, R r, Color const& color)
{
   R2 qc=natTrn(pc);
   Frame::fillDisk(qc.x1,qc.x2,r,color);
}

Z Graph::bundleClip(const R_Vector& a, const R_Vector& b,
      const R_Vector& c, Rn2& yClip)
{
   Rn2 x(2), xc(6);
   R2 x_old(CL.x1,CL.x2),x_new;
   Z k=0, n=a.dim();
   Rn2 y(2*n+N_EXCESS);

   Z i,i1=0,i2=1,i1_,j,p;
   for (i=1;i<=n;i++){
      j=line_bound_inters(a[i],b[i],c[i],x);
      if (j==2){
         x_new=x[i1];
         j=connect_on_boundary(x_old,x_new,xc);
         for (p=0;p<=j;p++){  // note that xc gives p+1 points
            y[k]=xc[p];
            k++;
         }
         y[k]=x[i2];
         k++;
         x_old=x[i2];
         i1_=i1; i1=i2; i2=i1_; // exchanging i1 and i2
      }
   }
   y=y.resize(k);
   yClip=y;
   return k;
}

void Graph::show(const R_Vector& a, const R_Vector& b, const R_Vector& c,
   const Color& color)
{
   Rn2 x(2);
   Z i,j, n=a.dim();

   for (i=1;i<=n;i++){
      j=line_bound_inters(a[i],b[i],c[i],x);
      if (j==2){
         draw(x[0],x[1],color);
      }
   }
   vis(writeOnShow);
}

void Graph::showConnected(const R_Vector& a,
   const R_Vector& b, const R_Vector& c, const Color& color)
{
   Z i,k;

   Rn2 y;
   k=bundleClip(a,b,c,y);
   for (i=0;i<k-1;i++){
      draw(y[i],y[i+1],color);
   }
   vis(writeOnShow);
}

void Graph::gridL(Z n_x, Z n_y)
{
   if (cpmverbose>4) cpmmessage("gridL() started");
   if (n_x<1) n_x=1;
   if (n_y<1) n_y=1;

   R h_x=1./n_x;
   R h_y=1./n_y;
   Z i,k;

   for (i=1; i<n_x; i++){ // vertical lines, lines of constant x
      k=ix(R(i)*h_x);
      line(k, pL.x2, k, pU.x2, gridColor);
   }

   for (i=1; i<n_y; i++){ // horizontal lines, lines of constant y
      k=jy(R(i)*h_y);
      line(pL.x1, k, pU.x1, k, gridColor);
   }

// vertical border lines
   line(pL.x1, pL.x2, pL.x1, pU.x2, gridColor);
   line(pU.x1, pL.x2, pU.x1, pU.x2, gridColor);

// horizontal border lines
   line(pL.x1, pL.x2, pU.x1, pL.x2, gridColor);
   line(pL.x1, pU.x2, pU.x1, pU.x2, gridColor);
   // here it is evident that we have a clean grid with no
   // 1-pixel gaps possible if lines actually fill their starting
   // and end pixel
}

void Graph::draw(const R_Vector& x, const R_Vector& y, const Color& c)
{
   Z mL=3;
   Word loc="Graph::draw(R_Vector,R_Vector,Color)";
   CPM_MA
   Z n=cpminf<Z>(x.dim(),y.dim());
   if (n==0){
      cpmwarning(loc&": n==0, nothing to draw");
   }
   else if (n==1)
      mark(R2(x[1],y[1]),c);
   else{
      for (Z i=1;i<n;++i) draw(R2(x[i],y[i]),R2(x[i+1],y[i+1]),c);
   }
   CPM_MZ
}

Z2 Graph::adjXY(Z divx, Z divy)
{
   Z mL=3;
   Word loc="Graph::adjXY(Z,Z)";
   CPM_MA
   R2 a;
   Z nX=lad_4(CL.x1,CU.x1,divx,a);
   CL.x1=a.x1;
   CU.x1=a.x2;
   Z nY=lad_4(CL.x2,CU.x2,divy,a);
   CL.x2=a.x1;
   CU.x2=a.x2;
   updatedXdY();
    // nX and nY are assured not to be 0
  // setAutoScale(false);
      // we don't like to destroy the carefully adjusted window in object space
   CPM_MZ
   return Z2(nX,nY);
}

void Graph::clearAndGrid(Z nX, Z nY)
{
   Z mL=3;
   Word loc("Graph::clearAndGrid(Z,Z)");
   CPM_MA
   Z nAcc=gridPrec;
      // parameter to ostringstream::setprecision, default value
      // was not always sufficient in my applications (added 2007-09-12)
      // also value 4 was not sufficient since it wrote 2008.5 as
      // 2008 in an application in which x was given in years.
      // Now this value is among the Graph's data and may be set at will.
   if (!noGridText){
      cpmcerr<<"Graph::clearAndGrid(Z nX, Z nY): noGridText = "<<noGridText<<endl;
      ostringstream ost1,ost2;
      ost1<<std::setprecision(nAcc);
      ost2<<std::setprecision(nAcc);
      if (nX==0 || noGrid){
         ost1<<toDouble(CL.x1)<<" < x < "<<toDouble(CU.x1);
      }
      else{
         R xStep=(CU.x1-CL.x1)/nX;
         ost1<<"x:"<<toDouble(CL.x1)<<"("<<toDouble(xStep)<<")"
         <<toDouble(CU.x1);
      }
      if (nY==0 || noGrid){
         ost2<<toDouble(CL.x2)<<" < y < "<<toDouble(CU.x2);
      }
      else{
         R yStep=(CU.x2-CL.x2)/nY;
         ost2<<"y:"<<toDouble(CL.x2)<<"("<<toDouble(yStep)<<")"
         <<toDouble(CU.x2);
      }
      R2 p1=gridTextPosition;
      R hT=gridTextHeight_;
      R dh=hT*2_R/h(); // shift of two text lines
      R2 delta_p(0_R,dh);
      R2 p2=p1-delta_p;
      R2 p3=p2-delta_p;

      Word w1=Word(ost1.str());
      Word w2=Word(ost2.str());
      Word w3=Word(" n = "&cpm(nPoints));

      TextCom t1(p1,w1,gridTextColor,hT);
      TextCom t2(p2,w2,gridTextColor,hT);
      TextCom t3(p3,w3,gridTextColor,hT);
      text_+=t1;
      text_+=t2;
      if (nPoints!=0_Z) text_+=t3;

   }
   if (newGrid && !noGrid){ // newGrid set to 1 in constructor
      // and normally never changed. Name misleading;
      // comes not only then in action if the previous
      // grid actually has to be changed
      /*if(!getMarkMode())*/ paint(colorActive); // experiment
      // in markmode one should be able to superimose graphics
      // on the other hand movies should not superimpose and this is
      // in my programs the dominant mode.
      gridL(nX,nY);
   }
   CPM_MZ
}

bool Graph::mark(const R_Vector& x, const R_Vector& y,
   const Color& c, Z bars, const Word& title, R xPosRel)
{
   Z mL=5;
   Word loc("Graph::show(R_Vector,R_Vector,Color,Z,Word,R)");
   CPM_MA
   if (title.dim()>0) addText(title, xPosRel);
   bool res=markGL(x,y,c,bars);
   CPM_MZ
   return res;
}

bool  Graph::mark(const R_Vector& y, const Color& c, Z bars,
    const Word& title, R xPosRel)
{
   Z mL=2;
   Word loc("Graph::mark(R_Vector,Color,Z,Word,R)");
   CPM_MA
   Z d=y.dim();
   R_Vector x(d);
   for (Z i=1;i<=d;i++) x[i]=i;
   bool res=mark(x,y,c,bars,title,xPosRel);
   CPM_MZ
   return res;
}

void Graph::addFrame(void)
{
   addFrm_(graphColor);
}

void Graph::addFrm_(Color const& cf)
{
   Z k=pL.x1;
   line(k,pL.x2,k,pU.x2,cf);
   k=pU.x1;
   line(k,pL.x2,k,pU.x2,cf);
   k=pL.x2;
   line(pL.x1,k,pU.x1,k,cf);
   k=pU.x2;
   line(pL.x1,k,pU.x1,k,cf);
}

bool Graph::mark(const R_Vector& x, const R_Matrix& y,
   const CpmArrays::V<Color>& cl, const Word& title, R xPosRel,
   bool xLog, bool yLog)
{
   Z mL=4;
   Word loc("Graph::mark(R_Vector,R_Matrix,V<Color>,Word,R,bool,bool)");
   CPM_MA
   if (title.dim()>0) addText(title, xPosRel);
      // Notice that setText cancels all previously set text
      // thus also the grid text. addText avoids this inconvenience.
   nPoints=x.dim(); // I had a case where Graph::nPoints helpfully shined up in the
   // legend of a diagram, although it was not registers at the point of
   // creation. Here nPoints was a memorized quantity.
   // cpmcerr<<" nPoint = "<<nPoints<<endl;
   bool res=markGL(x,y,cl,xLog,yLog);
   CPM_MZ
   return res;
}

bool Graph::markCf(R_Vector const& x, R_Matrix const& y, Iv const& ivX,
    CpmArrays::V<Color> const& cl, Word const& title, R xRelTitle,
    bool xLog, bool yLog)
{
   Z mL=4;
   Word loc("Graph::markCf(R_Vector const& x, R_Matrix const& y ...)");
   CPM_MA
   Z2 z2=x.indCon(ivX);
   R_Vector xcf=x.confine(z2);
   R_Matrix ycf=y.confine(z2);

   if (autoScaleY){
      Iv ivy=ycf.infSup();
      setY(ivy);
   }
   CPM_MZ
   return mark(xcf,ycf,cl,title,xRelTitle,xLog,yLog);
}

bool Graph::markTr(R_Vector x, R_Matrix const& y, CpmArrays::V<Color> const& cl)
// Tr transparent
{
   Z mL=3;
   Word loc("Graph::markTr(R_Matrix,V<Color>)");
   CPM_MA
   if (autoScale){
      bool b=setXY(x,y);
      if (b==false){
         cpmwarning(loc&": no graph could be created");
         return false;
      }
      adjXY_(divX,divY); // here divX, divY are the actual values
         // which probably will be changed by functuion adXY_
      clearAndGrid(divX,divY);
   }
   else{
      Word xL=cpm(CL.x1);
      Word xU=cpm(CU.x1);
      Word yL=cpm(CL.x2,3);
      Word yU=cpm(CU.x2,3);
   }
   Vo<Color> clAct(cl);

   V<Color> exceptions(2);
   exceptions[1]=colorActive;
   exceptions[2]=gridColor;
   clAct=clAct.purge(exceptions);
      // clAct may be void after purge (2006-02-05)
   Z cAd=clAct.dim();
   if (cAd==0){
      cpmwarning(loc&
         ": color array was void after purging, replaced by RED"
      );
      clAct=V<Color>{RED};
   }
   Z m=y.dim1();
   for (Z i=1;i<=m;i++){
      draw(x,y[i],clAct(i,CYCLIC));
   }
   CPM_MZ
   return true;
}

bool Graph::markGL(R_Vector const& xi, R_Matrix const& yi,
   CpmArrays::V<Color> const& cl, bool xLog, bool yLog)
{
   Z mL=3;
   Word loc("Graph::markGL(R_Vector,R_Matrix,V<Color>)");
   CPM_MA
   R_Vector x=xLog ? xi.log10() : xi;
   R_Matrix y=yLog ? yi.log10() : yi;

   if (x.dim()!=y.dim2())  cpmerror (loc&": dimension mismatch");

   if (autoScale){
      bool b=setXY(x,y);
      if (b==false){
         cpmwarning(loc&": no graph could be created");
         return false;
      }
      adjXY_(divX,divY);
      clearAndGrid(divX,divY);
   }
   else{
      Word xL=cpm(CL.x1);
      Word xU=cpm(CU.x1);
      Word yL=cpm(CL.x2,3);
      Word yU=cpm(CU.x2,3);
     // clearAndGrid(10,10); // was 0,0
   }
   Vo<Color> clAct(cl);

   V<Color> exceptions(2);
   exceptions[1]=colorActive;
   exceptions[2]=gridColor;
   clAct=clAct.purge(exceptions);
      // clAct may be void after purge (2006-02-05)
   Z cAd=clAct.dim();
   if (cAd==0){
      cpmwarning(loc&
         ": color array was void after purging, replaced by RED"
      );
      clAct=V<Color>{RED};
   }
   Z m=y.dim1();
   for (Z i=1;i<=m;i++){
       draw(x,y[i],clAct(i,CYCLIC));
   }
   CPM_MZ
   return true;
}

bool Graph::markGL(R_Vector const& x, R_Vector const& y,
   Color const& c, Z bars)
{
   Z mL=5;
   Word loc("Graph::markGL(R_Vector,R_Vector,Color,Z)");
   CPM_MA
   if (x.dim()!=y.dim())
      cpmerror (loc&": dimension mismatch");

   if (autoScale){
      bool b=setXY(x,y);
      if (b==false){
         cpmwarning(loc&": no graph could be created");
         return false;
      }
      adjXY_(divX,divY);
      clearAndGrid(divX,divY);
   }
   else{
      Word xL=cpm(CL.x1);
      Word xU=cpm(CU.x1);
      Word yL=cpm(CL.x2,3);
      Word yU=cpm(CU.x2,3);
      adjXY_(divX,divY);
      clearAndGrid(divX,divY);
      // was clearAndGrid(0,0);
   }
   Color cFinal=( c==colorActive ? graphColor : c);
   if (bars==0){
       draw(x,y,cFinal);
   }
   else{
      Z n=x.dim();
      cpmassert(n==y.dim(),loc);
    //  bool thick=true;
    //  bool limited=true;
      for (Z i=1;i<=n;i++){
         R2 p1(x[i],0.);
         R2 p2(x[i],y[i]);
         draw(p1,p2,cFinal);
      }
   }
   CPM_MZ
   return true;
}

void Graph::mark1(const R_Matrix& a, const CpmArrays::V<Color>& cl,
    const Word& title, R xPosRel)
{
   Z mL=3;
   Word loc("Graph::mark1(const R_Matrix& , ...)");
   CPM_MA
   Z i,ma=a.dim1(),na=a.dim2();
   if (ma<1){
      cpmwarning("could not show a matrix of 0 lines");
   }
   else if (na<1){
      cpmwarning("could not show a matrix of 0 columns");
   }
   else{
      R_Vector x=a[1];
      Z my=ma-1;
      R_Matrix y(my);
      for (i=1;i<=my;++i) y[i]=a[i+1];
      mark(x,y,cl,title,xPosRel);
   }
   CPM_MZ
}

namespace{

   void lin(R2 p1, R2 p2, ostream& str)
   {
    //  Z np=1;
      str<<endl<<"   "<<"\\drawline[0]("<<p1[1]<<","<<p1[2]<<
         ")("<<p2[1]<<","<<p2[2]<<")";
   }

   void text(R2 pos, Word const& txt, R ox, R oy, ostream& str)
   {
      R px=pos[1]+ox;
      R py=pos[2]+oy;
      str<<endl<<"   \\put("<<
         px<<","<<py<<"){\\makebox(0,0)[bl]{"<<txt<<"}}";
   }

   void sym(R2 p, Z s, ostream& str)
   {
      if (s==-1)
         str<<endl<<
         "   \\put("<<p[1]<<","<<p[2]<<
         "){\\makebox(0,0){\\tiny $\\blacksquare$}}";
      else if (s==-2)
         str<<endl<<
         "   \\put("<<p[1]<<","<<p[2]<<
         "){\\makebox(0,0){\\tiny $\\blacktriangle$}}";
      else
         str<<endl<<
         "   \\put("<<p[1]<<","<<p[2]<<
         "){\\makebox(0,0){\\tiny $\\bullet$}}";
   }
   // -1 and -2 experiment, 1 and 2 are original
}

bool Graph::makeLaTexGraph(R_Matrix const& a, Fig const& fg,
   ostream& str, Iv px, Iv py)
{
   Z mL=3;
   Word loc("Graph::makeLaTexGraph(R_Matrix ...)");
   CPM_MA
   R hLeg=7;
   R width=px.abs();
   R height=py.abs();
   R raster=hLeg*2.5;

   Z nX=cpmrnd(width/raster);
   Z nY=cpmrnd((height-hLeg)/raster);
   Word grid_text=gridText[1];
   if (autoScale){
      cpmmessage(mL,"block autoScale started");
      setXY(a);
      R2 aa;
      Z divX=nX,divY=nY;
      nX=lad_4(CL.x1,CU.x1,divX,aa);
      CL.x1=aa.x1;
      CU.x1=aa.x2;
      nY=lad_4(CL.x2,CU.x2,divY,aa);
      CL.x2=aa.x1;
      CU.x2=aa.x2;
      updatedXdY();
      // nX and nY are assured not to be 0
      R xStep=(CU.x1-CL.x1)/nX;
      R yStep=(CU.x2-CL.x2)/nY;
      ostringstream ost;
      Z ad=a.dim();
      if (ad==2){
         ost<<"$x="<<fg.abs<<"\\;:\\;"<<CL.x1<<"("<<xStep<<")"<<CU.x1<<
         " \\quad \\quad  y="<<fg.ord[1]<<" (\\bullet)"
         <<"\\;:\\;"<<CL.x2<<"("<<yStep<<")"<<CU.x2<<"$";
         grid_text=Word(ost.str());
      }
      else if (ad==3){
         ost<<"$x="<<fg.abs<<"\\;:\\;"<<CL.x1<<"("<<xStep<<")"<<CU.x1<<
         " \\quad \\quad  y="<<fg.ord[1]<<
         " (\\blacksquare), \\:"<<fg.ord[2]<<" (\\blacktriangle) "
         <<"\\;:\\;"<<CL.x2<<"("<<yStep<<")"<<CU.x2<<"$";
         grid_text=Word(ost.str());
      }
      else ;

      cpmmessage(mL,"block autoScale done");
   }
   Iv rx=ranX();
   Iv ry=ranY();
   Iv pyU(py[1]+hLeg,py[2]); // U means used (for graphics)
   R_Func fx(rx,px);
   R_Func fy(ry,pyU);

   R ox=px[1];
   R oy=py[1];
   R lx=px[2]-ox;
   R ly=py[2]-oy;

   V<R> vpx=px.div(nX+1);
   V<R> vpy=pyU.div(nY+1);
   Z dx=vpx.dim();
   Z dy=vpy.dim();
   Z i,j;

   str.setf(std::ios_base::fixed,std::ios_base::floatfield);
   str<<endl<<"% Output of function Graph::makeLaTexGraph(...)";
   str<<endl<<"\\setlength{\\unitlength}{1mm}";
   str<<endl<<"\\begin{figure}";
   str<<endl<<"\\begin{center}";
   str<<endl<<"\\begin{picture}("<<lx<<","<<ly<<")("<<ox<<","<<oy<<")";
   str.precision(6);
   str<<endl<<"\\linethickness{0.01mm}";
   // making the grid
   for (i=1;i<=dx;++i) lin(R2(vpx[i],vpy[1]),R2(vpx[i],vpy[dy]),str);
   for (j=1;j<=dy;++j) lin(R2(vpx[1],vpy[j]),R2(vpx[dx],vpy[j]),str);
   // writing the gridtext
   lin(R2(vpx[1],pyU[1]),R2(vpx[dx],pyU[1]),str);
   lin(R2(vpx[1],py[1]),R2(vpx[dx],py[1]),str);
   lin(R2(vpx[1],py[1]),R2(vpx[1],pyU[1]),str);
   lin(R2(vpx[dx],py[1]),R2(vpx[dx],pyU[1]),str);
   text(R2(2,1),grid_text,ox,oy,str);

   Z m=a.dim1();
   Z n=a.dim2();
   for (i=2;i<=m;++i){
      Z sy=i-1;
      for (j=1;j<=n;++j){
         R2 p1(fx(a[1][j]),fy(a[i][j]));
         sym(p1,sy,str);

      }
   }

   str<<endl<<"\\end{picture}";
   str<<endl<<"\\caption{"<<fg.tit<<"}";
   str<<endl<<"\\end{center}";
   str<<endl<<"\\end{figure}";
   CPM_MZ
   return true;
}

bool Graph::show(const F<R,R>& f, Z n, const Color& cf,
   const V<R>& xList, const Color& cx)
{
   using CpmGeo::Iv; // interval
   Z mL=3;
   static Word loc("show(F<R,R>,Z,Color,V<R>,Color)");
   CPM_MA
   if (n<=1) n=2;
   Z i;
   Iv ivx=getIvX();
   R_Vector x(ivx,n);
   R_Matrix y(2,n);
   for (i=1; i<=n; i++){
      y[1][i]=x[i];
      y[2][i]=f(x[i]);
   }
   V<Color> vcf(1,cf);
   show(y,vcf);
// now vertical lines for selected x-values
// for analysis of curve behavior
   Z nx=xList.dim();
   if (nx==0){
      CPM_MZ
      return false;
   }
   else{
      Iv ivy=getIvY();
      R y1=ivy[1];
      R y2=ivy[2];
      for (i=1;i<=nx;i++){
         R xi=xList[i];
         if (ivx.contains(xi)) draw(R2(xi,y1),R2(xi,y2),cx);
      }
   }
   CPM_MZ
   return true; // trivial return
}

bool Graph::mark(const F<R,R>& f, Z n, const Color& cf)
{
   using CpmGeo::Iv; // interval
   Z mL=4;
   static Word loc("mark(F<R,R>,Z,Color)");
   CPM_MA
   if (n<=1) n=2;
   Z i;
   Iv ivx=getIvX();
   R_Vector x(ivx,n);
   R_Matrix y(2,n);
   for (i=1; i<=n; i++){
      y[1][i]=x[i];
      y[2][i]=f(x[i]);
   }
   V<Color> vcf(1,cf);
   mark1(y,vcf); // not mark !
   CPM_MZ
   return true; // trivial return
}

void  Graph::mark( V<R> const& xList, V<Color> const& cList)
{
   Z nx=xList.dim();
   if (nx==0) return;
   Iv ivx=getIvX();
   Iv ivy=getIvY();
   R y1=ivy[1];
   R y2=ivy[2];
   for (Z i=xList.b();i<=xList.e();++i){
      R xi=xList.cui(i);
      if (ivx.contains(xi))
         draw(R2(xi,y1),R2(xi,y2),cList(i,CYCLIC));
   }
}

void Graph::show(Z n, R pl, R pu,
                R f1(R z1), R f2(R z2), const Color& c)
{
   Z i;
   R h;
   if (n<=1) n=2;
   R_Vector x(n);
   R_Vector y(n);
   h=(pu-pl)/(n-1);
    R p_=pl;
   for (i=1; i<=n; i++){
      x[i]=f1(p_); y[i]=f2(p_); p_+=h;
   }
   show(x,y,c);
}

// Common code for implementing the next two functions.
// An attempt to use a function taking a function object as argument
// failed since Graph would be a parameter of such a function object
// and had to be constant in this role. (2001-4-20)

#define CPM_SHR_COD1\
   R2 cLowerLeft(XL,YL);\
   R2 cUpperRight(XU,YU);\
   Z2 pU_(iX(XU),jY(YU));\
   Z2 pL_(iX(XL),jY(YL));\
   R dy_=YU-YL;\
   if (!(dy_>0.)) cpmerror("Graph::forAll(...): dy_<=0.");\
   R dx_=XU-XL;\
   if (!(dx_>0.)) cpmerror("Graph::forAll(...): dx_<=0.");\
   Z nRow_=pL_.x2-pU_.x2+1;\
   if (!(nRow_>0)) cpmerror("Graph::forAll(...): !(nRow_>0)");\
   Z nCol_=pU_.x1-pL_.x1+1;\
   if (!(nCol_>0)) cpmerror("Graph::forAll(...): !(nCol_>0)");\
   Z i,j;\
   R yStep=dy_/nRow_;\
   R xStep=dx_/nCol_;\
   R xStart=cLowerLeft.x1;\
   R yStart=cUpperRight.x2;\
   xStart+=xStep*0.5;\
   yStart-=yStep*0.5;\
   R xAct=xStart;\
   R yAct=yStart;\
   Z counter=0; R dyInv=R(1)/nRow_;\
   for (j=pU_.x2;j<=pL_.x2;j++){\
      counter++; vis();cpmprogress("line by line",dyInv*counter);\
      for (i=pL_.x1;i<=pU_.x1;i++){


#define CPM_SHR_COD2\
          putPixel(i,j,Color(vr,vg,vb));\
          xAct+=xStep;\
      }\
      xAct=xStart;\
      yAct-=yStep;\
   }\
   if (!noDisplay){\
      vis(writeOnShow);\
   }

void Graph::show(const F<RR,V<R> >& f, R XL, R XU, R YL, R YU,
   bool noDisplay )
{
   CPM_SHR_COD1
          RR xy(xAct,yAct);
          V<R> val=f(xy);
          R vr=val[1];
          R vg=val[2];
          R vb=val[3];
   CPM_SHR_COD2
}

void Graph::show(const F<R2, R3 >& f, R XL, R XU, R YL, R YU,
   bool noDisplay )
{
   Z mL=3;
   Word loc("Graph::show(F<R2,R3>,R,R,R,R,bool)");
   CPM_MA
   CPM_SHR_COD1
          R2 xy(xAct,yAct);
          R3 val=f(xy);
          R vr=val[1];
          R vg=val[2];
          R vb=val[3];
   CPM_SHR_COD2
   CPM_MZ
}

void Graph::show(const R2_Func& fR, const R2_Func& fG, const R2_Func& fB,
   R XL, R XU, R YL, R YU, bool noDisplay )
{
   CPM_SHR_COD1
          RR xy(xAct,yAct);
          R vr=fR(xy);
          R vg=fG(xy);
          R vb=fB(xy);
   CPM_SHR_COD2
}

void Graph::show(const R2_Func& fR, const R2_Func& fG, const R2_Func& fB,
   const Response& rR, const Response& rG, const Response& rB,
   Iv roiX, Iv roiY, bool noDisplay )
{
   Z mL=3;
   Word loc("Graph::show(R2_Func,..Response,..Iv,Iv)");
   cpmmessage(mL,loc&" started");
   R XL=roiX.inf();
   R XU=roiX.sup();
   R YL=roiY.inf();
   R YU=roiY.sup();
   CPM_SHR_COD1
          RR xy(xAct,yAct);
          R vr=rR(fR(xy));
          R vg=rG(fG(xy));
          R vb=rB(fB(xy));
   CPM_SHR_COD2
   cpmmessage(mL,loc&" done");
}

#undef CPM_SHR_COD2

void Graph::show(R2_Func const& f, Response const& rsp, Iv roix, Iv roiy)
{
   Z mL=3;
   Word loc("Graph::show(R2_Func,Response,Iv,Iv)");
   CPM_MA
   Iv roiX=(roix.isVoid() ? getIvX() : roix);
   Iv roiY=(roiy.isVoid() ? getIvY() : roiy);
   R XL=roiX.inf();
   R XU=roiX.sup();
   R YL=roiY.inf();
   R YU=roiY.sup();
   CPM_SHR_COD1
      RR xy(xAct,yAct);
      R val=f(xy);
      Color c=rsp.color(val);
      putPixel(i,j,c);
          xAct+=xStep;
      }
      xAct=xStart;
      yAct-=yStep;
   }
   vis(writeOnShow);
   CPM_MZ
}

#undef CPM_SHR_COD1


void Graph::show(F<C,C> const& f, R gamma, bool zeroIsZero)
{
   using CpmImaging::Color;
   R rMax=-1;
   //R rMin=CpmRootX::hugeNumber;
   R rMin=1.e64;

   V<R> cx=centersX();
   V<R> cy=centersY();
   Z nx=dim_x();
   Z ny=dim_y();
   VV<C> val(nx,ny);
   Z i,j;
   C fz;
   // notice that our indexes start with 1
   for (i=1;i<=nx;i++){
      for (j=1;j<=ny;j++){
         C z(cx[i],cy[j]);
         fz=f(z);
         val[i][j]=fz;
         R r=fz.abs();
         if (r>rMax) rMax=r;
         if (r<rMin) rMin=r;
      }
   }
   if (zeroIsZero) rMin=0;
   Iv range(rMin,rMax,true);
   Color::setGamma(gamma);
   for (i=1;i<=nx;i++){
      for (j=1;j<=ny;j++){
         fz=val[i][j];
         Color c(fz,range); // depends on gamma
         put_xy(i,j,c); // simpel indexing due to
            // a newly defined function put_xy
      }
   }
   vis(writeOnShow);
}

void Graph::show(F<R2,C> const& f, R gamma, bool zeroIsZero)
{
   using CpmImaging::Color;
   R rMax=-1;
   //R rMin=CpmRootX::hugeNumber;
   R rMin=1.e64;

   V<R> cx=centersX();
   V<R> cy=centersY();
   Z nx=dim_x();
   Z ny=dim_y();
   VV<C> val(nx,ny);
   Z i,j;
   C fz;
   // notice that our indexes start with 1
   for (i=1;i<=nx;i++){
      for (j=1;j<=ny;j++){
         R2 z(cx[i],cy[j]);
         fz=f(z);
         val[i][j]=fz;
         R r=fz.abs();
         if (r>rMax) rMax=r;
         if (r<rMin) rMin=r;
      }
   }
   if (zeroIsZero) rMin=0;
   Iv range(rMin,rMax,true);
   Color::setGamma(gamma);
   for (i=1;i<=nx;i++){
      for (j=1;j<=ny;j++){
         fz=val[i][j];
         Color c(fz,range); // depends on gamma
         put_xy(i,j,c); // simpel indexing due to
            // a newly defined function put_xy
      }
   }
   vis(writeOnShow);
}

void Graph::show(const R2_Func& fR, const R2_Func& fG, const R2_Func& fB)
{
   show(fR,fG,fB,CL.x1,CU.x1,CL.x2,CU.x2);
}

void Graph::show(const F<RR,V<R> >& f)
{
   show(f,CL.x1,CU.x1,CL.x2,CU.x2);
}

bool Graph::mark(V< V<R2> > const& dat, V<Color> const& cl, R secFac,
B adaptToData)
{
   Z mL=3;
   Word loc("Graph::mark(V<V<R2> >, V<Color>,R)");
   CPM_MA
   static Z2 steps;
   if (adaptToData){
      adaptToData=false;
      bool b=setXY(dat,1_R,secFac);
      if (b==false){
         cpmwarning(loc&": no graph could be created");
         CPM_MZ
         return false;
      }
      adjXY_(divX,divY);
      clearAndGrid(divX,divY);
   }
   else{
      adjXY_(divX,divY);
      clearAndGrid(divX,divY);
   }
   Vo<Color> clAct(cl);
   V<Color> exceptions(2);
   exceptions[1]=colorActive;
   exceptions[2]=gridColor;
   clAct=clAct.purge(exceptions);
   Z cAd=clAct.dim();
   if (cAd<cl.dim()){
      cpmwarning("colors changed due to gridColor or colorActive",true);
      cpmwait(4.,2);
   }
   if (cAd==0){
      cpmwarning(loc&
         ": color array was void after purging, replaced by RED"
      );
      clAct=V<Color>{RED};
   }
   vector<R2>::const_iterator j,jF;
   for (Z i=dat.b();i<=dat.e();i++){ // loop over the colors
      Color ci=clAct(i,CYCLIC);
      //Color ci=cl[i];
      vector<R2> dati=dat[i].std();
      j=dati.begin();
      jF=dati.end();
      while (j!=jF){ // elegant iteration
          R2 p1=*j++; // now j is the next one
          if (j!=jF) draw(p1,*j,ci); // j valid if used
      }
   }
//   Z dim1=dat.dim(); ~ 50% slower
//   Z dim2=dat[1].dim();
//   for (Z i=1;i<=dim1;i++){ // loop over the colors
//      Color ci=clAct(i,CYCLIC);
//      for (Z j=2;j<=dim2;++j){
//         draw(dat[i][j-1],dat[i][j],ci);
//      }
//   }
   CPM_MZ
   return true;
}

bool Graph::mark(R_Matrix const& dat, V<Color> const& cl,
                 R secFac, Z modeXY)
// This is efectively the previous function where the array of type
// V< V<R2> > to be inserted as argument into the previous function
// is formed from matrix dat in a manner that depends on the argument
// modeXY.
// Matrix dat is assumed to have the following layout:
//
// t1 x11 y11 x12 y12 x13 y13 ... x1n y1n
//    ...
// tm xm1 ym1 xm2 ym2 xm3 ym3 ... xmn ymn
//
// Data of this form arise when a system of n-particles is simulated in
// the course of times (t1,t2,...,tm). Here (x,y) is a 2D position of a
// particle, which may result from a 3D position by projection.
// A convenient visualization should collect the (x,y)-pairs of any
// selected particle, ordered  according to their associated time ti
// as a curve. This curve thus results from the polygon
//    Pi := ((x1i,y1i),(x2i,y2i),...,(xmi,ymi)),
// where the points belong to the times (t1,t2,...,tm).
// Representation of these polygons is done for modeXY==0.
// For modeXY==1 we consider and represent the polygons
//    Xi := ((t1,x1i),(t2,x2i),...,(tm,xmi))
// and for modeXY==2 we consider and represent the polygons
//    Yi := ((t1,y1i),(t2,y2i),...,(tm,ymi)).
{
   Z mL = 1;
   Word loc("Graph::mark(R_Matrix, V<Color>,R,Z)");
   CPM_MA
   Z m=dat.dim1(), n=dat.dim2(),i,j;
   cpmassert(m>0,loc);
   cpmassert(n>0,loc); // n needs to be a positive odd integer
   n-=1; // now n is to be even
   cpmassert(n%2==0,loc); // n needs to be even
   n/=2; // now n is the n

   V< V<R2> > dat2(n, V<R2>(m));
   if (modeXY==1){
      for (i=1;i<=m;++i){
         for (j=1;j<=n;++j){
            Z j2=j*2;
            dat2[j][i]=R2(dat[i][1],dat[i][j2]);
         }
      }
   }
   else if (modeXY==2){
      for (i=1;i<=m;++i){
         for (j=1;j<=n;++j){
            Z j2=j*2+1;
            dat2[j][i]=R2(dat[i][1],dat[i][j2]);
         }
      }
   }
   else{ // regularly modeXY==0
      for (i=1;i<=m;++i){
         for (j=1;j<=n;++j){
            Z j2=j*2;
            dat2[j][i]=R2(dat[i][j2],dat[i][j2+1]);
         }
      }
   }
   bool res=mark(dat2,cl,secFac);
   CPM_MZ
   return res;
}