//========================================

#include "bib_fred.h"


//---------------------------------------------------------
/*---ecrit le contenu du vecteur dans un fichier --
format : taille, contenu

*/


void Ecrit_fichier(cx_vec & V, ofstream & fich)
{

    fich<<V.n_rows<<endl;

    for(int i=0; i<V.n_rows; i++)
        fich<<V(i)<<endl;

}
//-----------------------------------------------
void Ecrit_fichier(cx_vec &V,const char * nom_fich)
{
    cout<<" Ecrit cx_vec dans fichier: "<<nom_fich<<endl;
    ofstream fich(nom_fich);
    Ecrit_fichier(V,fich);

}


//------------------------------------------
void Lit_fichier(cx_vec &V,ifstream & fich)
{

    int size;
    fich>>size;

    V.set_size(size);

    for(int i=0; i<size; i++)
        fich>>V(i);

}
//--------------------------------------------
void Lit_fichier(cx_vec &V,const char * nom_fich)
{
    ifstream fich(nom_fich);
    Lit_fichier(V,fich);
}


//============================
/*
Entrée:
 matrice complexe ou reelle M
Sortie:
 val_p, vect _p  ordonnées par module croissant
 */
/*
void Diago(cx_vec & val_p,cx_mat vect_p, cx_mat & M)
{
  eig_gen(val_p, vect_p, M); // diagonalise


  // cout<<"A3="<<A3<<endl<<" eigenvalues="<<val_p<<endl;
  // for(int i=0;i<10;i++)
  //   cout<<norm(val_p(i))<<", "<<flush;
  // cout<<endl;


  //... ordonne les valeurs propres (et vect p)  par module decroissant

  uvec indices = sort_index(abs(val_p)); // liste des indices
  val_p=val_p(indices); // vecteur ordonné
  vect_p=vect_p.cols(indices); // vecteur ordonné
}


void Diago(cx_vec & val_p,cx_mat vect_p, mat & M)
{
  eig_gen(val_p, vect_p, M); // diagonalise


  // cout<<"A3="<<A3<<endl<<" eigenvalues="<<val_p<<endl;
  // for(int i=0;i<10;i++)
  //   cout<<norm(val_p(i))<<", "<<flush;
  // cout<<endl;


  //... ordonne les valeurs propres (et vect p)  par module decroissant

  uvec indices = sort_index(abs(val_p)); // liste des indices
  val_p=val_p(indices); // vecteur ordonné
  vect_p=vect_p.cols(indices); // vecteur ordonné
}
*/
//=======================================


//=========================== GRAPHISME ============================

#include <TCanvas.h>
#include <TGraph.h>
#include <TPolyMarker.h>
#include <TEllipse.h>
#include <TH1.h>
#include <TROOT.h>

#include <TStyle.h>
#include <TView.h>

#include <TPolyMarker3D.h>

//========= dessin ========================
TH1D*  Dessin(vec &V,const char * titre,int superp,const char *opt,int opt_log)
{
    return Dessin(V,"","",1,V.size(),titre,superp,opt,-1111,-1111,1,opt_log);
}

//========= dessin ========================
TH1D*  Dessin(vec &V,const char * titre,const char *opt)
{
    return Dessin(V,"","",1,V.size(),titre,0,opt);
}
//========= dessin ========================
TH1D*  Dessin(vec &V,const char * xtitre,const char * ytitre,double xmin,double xmax,const char * titre, int superp,const char *opt,double ymin, double ymax,int col,int opt_log)
{



    //... recherche si il y a deja une fenetre "titre"...
    TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject(titre);
    if (c==0) // il n'y en a pas déjà, crée nouvelle fenetre
    {
//		cout<<"cree nouvelle fenetre dans vecteur.cc, titre="<<titre<<endl;
        c = new TCanvas(titre,titre,400,400);
        c->SetFillColor(0);
        superp=0;
    }
    else
    {
//		cout<<"fenetre deja existante:"<<titre<<endl;
    }



    return Dessin(V,xtitre, ytitre,xmin, xmax, titre, c,superp,opt,ymin,ymax,col,opt_log);
}
//========= dessin ========================
TH1D *  Dessin(vec & V,TCanvas * c)
{
    return Dessin(V,"","",1,V.size(),"Vecteur", c,0);
}
//========= dessin ========================
//  c : fenetre ou dessiner
// superp: =0 si premier figure, 1,2.. si superposition (style de ligne)
// ymin, ymax pour fixer. (-1111 si libre)
TH1D*  Dessin(vec& V,const char * xtitre,const char * ytitre,double xmin,double xmax,const char * titre,TCanvas * c,int superp,const char* opt,double ymin, double ymax,int col,int opt_log)
{
    int nx(V.size());



    ostringstream titre_histo;
    titre_histo<<titre<<superp<<ends;

    TH1D *h1=(TH1D*) gROOT->FindObject(titre_histo.str().c_str());

    if(h1) // "titre" existe dejà
    {
//		cout<<"on detruit h1 dans vecteur.cc"<<endl;
        delete h1;
    }


    h1=new TH1D(titre_histo.str().c_str(),titre,nx,xmin,xmax);



    h1->SetStats(kFALSE); // pas de panneau stat
    h1->SetMinimum(ymin);
    h1->SetMaximum(ymax);
    h1->SetXTitle(xtitre);
    h1->SetYTitle(ytitre);


    for (int  i=0; i<nx; i++) // remplissage
    {
        double z=V[i];
        h1->SetBinContent(i+1,z);

    }


//---- definit le marker P ----
    h1->SetMarkerColor(1);
    h1->SetMarkerStyle(6); // 3:etoile  6: petit cercle
    h1->SetLineStyle(1); // 1=solid, 2=dash, 3=dash-dot, 4=dot-dot
    h1->SetLineColor(col);

    c->cd();



    if (superp==0)
        h1->Draw(opt); // L: ligne  P: marker
    else
    {
        int style=1;

        if (superp>=1)
            style=(superp-1%4)+1;
        h1->SetLineStyle(1);//style); // 1=solid, 2=dash, 3=dash-dot, 4=dot-dot
        col=style;
        h1->SetLineColor(col);
        h1->Draw("same"); // L: ligne entre  P: marker
    }

    if (opt_log) c->SetLogy(1); //echelle LOG


    c->Update();
    return h1;
}
//========= dessin 2D (x[i],y[i] = objet)========================
TCanvas *  Dessin(vec &x,vec &y,const char * xtitre,const char * ytitre,const char * titre)
{

    if (x.size() !=y.size())
    {
        cerr<<"Il faut x.I=y.I dans Vecteur.cc"<<endl;
        return 0;
    }

    TCanvas *c = new TCanvas(titre,titre,400,400);
    c->SetFillColor(0);
    int nx(x.size());


    float *xx=new float[nx];
    float *yy=new float[nx];

    for (int i=0; i<nx; i++)
    {
        xx[i] =x[i];
        yy[i] =y[i];
    }


    TGraph  *gr= new  TGraph(nx,xx,yy);

    /*
    (gr->GetXaxis())->SetLabel(xtitre);
    (gr->GetYaxis())->SetLabel(ytitre);
    */
    c->cd();

    gr->Draw("AP"); // A: axis L:ligne , P=marker *: etoile

    c->Update();

    return c;
}

//========= dessin 3D (x[i] = objet, y[i], z[i])========================
TCanvas *  Dessin(vec &x,vec &y,vec &z, const char * xtitre,const char * ytitre,const char * ztitre,const char * titre)
{

    if ((y.size() !=x.size())||(z.size() !=x.size()))
    {
        cerr<<"Il faut taille y.I=I dans Vecteur.cc"<<endl;
        return 0;
    }

    TCanvas *c = new TCanvas(titre,titre,400,400);
    c->SetFillColor(0);
    int nx(x.size());


// creating view
    TView *view = TView::CreateView(1,0,0);

    double xmin=x.min(),ymin=y.min(),zmin=z.min(),xmax=x.max(),ymax=y.max(),zmax=z.max();
    view->SetRange( xmin,ymin,zmin,xmax,ymax,zmax);

    TPolyMarker3D *pm3d = new TPolyMarker3D(x.size());



    for( int i = 0; i <x.size(); i++ )
    {
        double xp=x[i], yp= y[i], zp= z[i] ;
        //double xp=0.5*(1.+cos(2.*Pi_*i/(1.*I))),  yp=0.5*(1.+sin(2.*Pi_*i/(1.*I))) , zp=0;
        pm3d->SetPoint( i,xp,yp,zp);
    }

    // set marker size, color & style
    pm3d->SetMarkerSize( 1 );
    pm3d->SetMarkerColor(2);
    pm3d->SetMarkerStyle( 1 );

    //draw
    c->cd();
    pm3d->Draw();





    c->Update();

    return c;
}
//========= dessin sur cercle========================
void  Dessin_phases(vec &V,const char * titre)
{
    TCanvas *c = new TCanvas(titre,titre,400,400);
    c->SetFillColor(0);
    int nx(V.size());


    float *x=new float[nx+2];
    float *y=new float[nx+2];

    for (int i=0; i<nx; i++)
    {
        x[i] = cos(V[i]);
        y[i] = sin(V[i]);
    }
    x[nx]=1;
    y[nx]=1; // pour cadrer
    x[nx+1]=-1;
    y[nx+1]=-1;

    TGraph  *gr= new  TGraph(nx+2,x,y);

    gr->Draw("A*");
    TEllipse *e=new  TEllipse(0,0,1);
    e->Draw();

    c->Update();

}


//========= dessin  sur colonne a l abcisse x========================
void  Dessin_col(vec &V,double xx, TCanvas *c,int col, int style)
{
    c->cd();
    float *x=new float[V.size()];
    float *y=new float[V.size()];

    double xm,ym,xM,yM;
    c->GetRange(xm,ym,xM,yM); //limites de la fenetre

    for (int i=0; i<V.size(); i++)
    {
        x[i] = xx;
        y[i] = V[i];

        if (y[i]>yM) y[i]=yM; //si dépassement vertical
        if (y[i]<ym) y[i]=ym;
    }



    TPolyMarker *pp= new TPolyMarker(V.size(),x,y);

    pp->SetMarkerColor(col);
    pp->SetMarkerStyle(style ); //6: point  8: petit cercle
    pp->Draw();


    /*
    TGraph  *gr= new  TGraph(I,x,y);
    gr->Draw("sameAC*");
    */
    //   c->Update();
}



//========= dessin ========================
// dessin des points dans le plan complexe
TCanvas *  Dessin(cx_vec &V,const char * titre)
{
//... recherche si il y a deja une fenetre "titre"...
    TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject(titre);
    if(c==0) // il n'y en a pas déjà, crée nouvelle fenetre
    {
        c = new TCanvas(titre,titre,400,400);
        c->SetFillColor(0);

    }
    else
    {
//      cout<<"fenetre deja existante:"<<titre<<endl;
        c->Clear();
    }



    c->cd();

	double  xmin = min((vec)real(V)), xmax=max((vec)real(V)), ymin=min((vec)imag(V)), ymax=max((vec)imag(V));

// double xmin(-2),xmax(2),ymin(-2),ymax(2);

    c->DrawFrame(xmin,ymin,xmax,ymax,titre);

    for (int i=0; i<V.size(); i++)
    {
        double x = (double)real(V[i]);
        double y= (double)imag(V[i]);
        point(x,y,kRed,8);// color,style



    }



    //-------- dessin du cercle unité -----------------
    circle(0,0,1);

    //-----------------------------------
    c->Update();

    return c;
}
/*
//==================================
// Dessine les histogrammes des phases et des modules des elements du vecteur
void Dessin_histo_phases_modules(cx_vec &V,double module_min, double module_max)
{
    cout<<"Nbre de valeurs ="<<V.size()<<endl;

    double phase_min=-Pi_, phase_max=Pi_;
    int nbin=V.size()/3+1;
    vec h_phases=((Phase(V)).Remplit_histo(phase_min,phase_max,nbin,"histo of phases"); // histogramme
    h_phases.Dessin("#phi","P",phase_min,phase_max,"phases");



   // int nbin=(*this).I/20+1;
    Vecteur mod=(V.Module());
   module_max=mod.Maximum();
   module_min=mod.Minimum();
    Vecteur h_modules=mod.Remplit_histo(module_min,module_max,nbin,"histo of modules"); // histogramme
    h_modules.Dessin("|z|","P",module_min,module_max,"modules");



}

*/





//========= Matrices ===================






//========= dessin ========================
TCanvas*  Dessin(mat &M,const char * titre,const  char *option)
{

    return Dessin(M,"i",1,M.n_cols,"j",1,M.n_rows,titre,option);
}
//========= dessin ========================
TCanvas*  Dessin(mat &M,const char * xtitre,double xmin,double xmax,const char * ytitre,double ymin,double ymax,const char * titre,const char * option,double zmin, double zmax)
{


//... recherche si il y a deja une fenetre "titre"...
    int superp=0;

    TCanvas *c;
    c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject(titre);

    if(c==0) // il n'y en a pas déjà, crée nouvelle fenetre
    {
        cout<<"cree fenetre dans matrice.cc"<<endl;
        c = new TCanvas(titre,titre,500,500);
        c->SetFillColor(0);

        superp=0;
    }
    else
    {
        cout<<"fenetre deja existante:"<<titre<<endl;
    }



    Dessin(M,xtitre,xmin,xmax,ytitre,ymin,ymax,titre,option,c,superp,zmin,zmax);

    return c;
}
//========= dessin ========================
//  c : fenetre ou dessiner
// superp: =0 si premiere figure, 1 si superposition
// zmin, zmax pour fixer. (-1111 si libre)
TCanvas*  Dessin(mat &M,const char * xtitre,double xmin,double xmax,const char * ytitre,double ymin,double ymax,const char * titre,const  char * option,TCanvas * c,int superp,double zmin, double zmax)
{
    //-------


    char options_2[30];
    strcpy(options_2,option);
    if(superp!=0)
    {
        strcat(options_2,"same");
        c->Clear();
    }

//-----------------------

    int nx(M.n_cols),ny(M.n_rows);



    TH2D *h=new TH2D(titre,titre,nx,xmin,xmax,ny,ymin,ymax);

    h->SetXTitle(xtitre);
    h->SetYTitle(ytitre);

    for( int i=0; i<M.n_rows; i++) // remplissage
        for(int j=0; j<M.n_cols; j++)
        {
            double z=M(i,j);


            h->SetCellContent(i+1,j+1,z);
        }
    //-----palette de gris -------pour:TH2F -SetContour(7); Draw("cont" ou "surf2");------
    /*
      Int_t colors[]={1,0,19,18,17,15,13,1,1}; // (ordre:des indices dans la palette 4,5,6,7,0,1,3)
      gStyle->SetPalette(8,colors); // classe TStyle
    h->SetContour(7);
    */

//gStyle->SetPalette(9,0); h->SetContour(10); // classe TStyle, palette de gris.
//gStyle->SetPalette(30,0); h->SetContour(30); // classe TStyle, palette de gris.


    gStyle->SetPalette(1,0); // jolie palette bleu->rouge
    gStyle->SetPadGridX(kTRUE); // grilles en X,Y
    gStyle->SetPadGridY(kTRUE);
    gPad->SetCrosshair(); // dessine une croix sous la souris



// h->SetEntries(1); // pour voir les lignes avec "cont".

//zmin=1e-4; cout<<"Min 1e-4 dans matrice.cc"<<endl;

// zmax=1; //@@
    h->SetMaximum(zmax);
    h->SetMinimum(zmin);

    c->Clear(); // efface les dessins sous-jacents
    h->SetStats(kFALSE); // pas de panneau stat

//..... en 3D: car cont ne marche pas!
    /*
    gPad->SetPhi(0.0);    // angle de position (action sur le pad actuel)
    gPad->SetTheta(90.0);
    h->Draw("asurf2");  // a: sans graduations
    */
//    c->SetLogz(1); // echelle log


//....... peu de divisions sur les axes ....
// h->GetXaxis()->SetNdivisions(5);  h->GetYaxis()->SetNdivisions(5);  h->GetZaxis()->SetNdivisions(5);

    c->cd();
    h->Draw(options_2);


//-------- dessin du cercle unité -----------------
// circle(0,0,1);


    c->Update();
    return c;
}




//--------------------------
/*
Dessine les lignes de niveau de || (z-P)^-1 ||_2  dans le plan complexe z.

Algo: ameliorer l'algo en  remplacant P par une matrice Triangulaire: P= O T O^-1
 avec O orthogonale,
ref: Threfethen's book.


 */
/*
void Dessin_Pseudo_Spectre(mat & M,double xmin,double xmax,double ymin,double ymax)
{


	int Nx(200), Ny(200); // >1, taille du tableau de z=x+iy

	mat R(Nx,Ny);// tableau de la norme de la resolvente

	int N=M.n_rows;
	if(M.n_cols!=N)
	{
		cout<<"Pb taille matrice dans armadillo.cc"<<endl;
		exit(1);
	}

	for(int ix=0;ix<Nx;ix++)
		for(int iy=0;iy<Ny;iy++)
		{
			double x=ix/(Nx-1.)*(xmax-xmin)+xmin;
			double y=iy/(Ny-1.)*(ymax-ymin)+ymin;
			Complex z=x+I_*y;

			mat Id=mat(N,N,fill::eye); // matrice identité
			vec S=svd(z*Id-M); // valeurs singulieres , ordre decroissant
			double smin=S(N-1); // la plus petite

			if(smin>0)
				R(ix,iy)=1./smin; //norme de la resolvente en z
			else
				// sinon c'est une val propre
				R(ix,iy)=0;
		} // for ix,iy


//....

	Dessin(R,"Re(z)",xmin,xmax,"Im(z)",ymin,ymax,"pseudo-spectre");




}

*/



//=====================================

// OPERATIONS SUR LES VECTOR





//=====================
// A typical alternative that forces a reallocation is to use swap:
/*
void Clear(vector & V)
{
vector<T>()->swap(V);
}
*/
//==========================
// Affiche vector
ostream & operator <<(ostream & sortie,  const vector<double> & V)
{

    sortie<<"["<<flush;

    for(int i=0; i<V.size(); i++)
    {
        if(i>0)
            sortie<<","<<flush;
        sortie<<V[i]<<flush;
    }
    cout<<"]"<<endl;

    return sortie;
}

//==========================
// Affiche vector
ostream & operator <<(ostream & sortie,  const vector<int> & V)
{

    sortie<<"["<<flush;

    for(int i=0; i<V.size(); i++)
    {
        if(i>0)
            sortie<<","<<flush;
        sortie<<V[i]<<flush;
    }
    cout<<"]"<<endl;

    return sortie;
}
//==========================
// Affiche vector
ostream & operator <<(ostream & sortie,  const vector<vector<int> > & L)
{
    for(int l=0; l<L.size(); l++)
    {
        sortie<<L[l];
    }

    return sortie;
}

//==========================
// Affiche vector
ostream & operator <<(ostream & sortie,  const vector<string> & V)
{

    sortie<<"["<<flush;

    for(int i=0; i<V.size(); i++)
    {
        if(i>0)
            sortie<<","<<flush;
        sortie<<V[i]<<flush;
    }
    cout<<"]"<<endl;

    return sortie;
}

//==========================
// Affiche vector
ostream & operator <<(ostream & sortie,  const vector<vector<string> > & L)
{
    for(int l=0; l<L.size(); l++)
    {
        sortie<<L[l];
    }

    return sortie;
}

//=============================
bool Subset(vector<int> V1,vector<int> V2)
{
    /*
    sortie:
       1 si V1 est un sous ensemble de V2
       0 sinon.

    Exemple:
       vector<int> V1; V1.push_back(2); V1.push_back(1);  V1.push_back(3);
       vector<int> V2; V2.push_back(1); V2.push_back(1);  V2.push_back(3);
       vector<int> V3; V3.push_back(3); V3.push_back(1);

       cout<<V1<<endl;    cout<<V2<<endl;    cout<<V3<<endl;

       cout<<"V1 in V2 :"<<Subset(V1,V2)<<endl;
       cout<<"V3 in V1 :"<<Subset(V3,V1)<<endl;
       cout<<"V3 in V2 :"<<Subset(V3,V2)<<endl;


     */

    sort(V1.begin(),V1.end());
    sort(V2.begin(),V2.end());
    //  cout<<"V1 sorted:"<<V1<<endl;
    //cout<<"V2 sorted:"<<V2<<endl;

    bool res= includes(V2.begin(),V2.end(), V1.begin(), V1.end());
//cout<<"res="<<res<<endl;
    return  res;
}

//===========================
vector<int> Sym_diff(vector<int> A,vector<int> B)
{
    /*
    Renvoit elements presents dans un seul vecteur (avec multiplicité)

    Exemple

    vector<int> V1; V1.push_back(2); V1.push_back(1);  V1.push_back(3);
    vector<int> V2; V2.push_back(1); V2.push_back(1);  V2.push_back(4);
    cout<<V1<<endl;
    cout<<V2<<endl;
    cout<<Sym_diff(V1,V2)<<endl;

    renvoit
    [2,1,3]
    [1,1,4]
    [1,2,3,4]
     */

    sort(A.begin(), A.end());
    sort(B.begin(), B.end());
    // allocate the smallest size of A,B as maximum size
    vector<int> c(A.size() < B.size() ? B.size() : A.size());
    vector<int>::iterator i;
    i = set_symmetric_difference(A.begin(), A.end(), B.begin(), B.end(), c.begin());
    return vector<int>(c.begin(), i);
}


//===========================
vector<int> diff(vector<int> A,vector<int> B)
{
    /*
    Renvoit elements presents A et pas dans B

    Exemple

    vector<int> V1; V1.push_back(2); V1.push_back(1);  V1.push_back(1);  V1.push_back(3);
    vector<int> V2; V2.push_back(1);   V2.push_back(4);
    cout<<V1<<endl;
    cout<<V2<<endl;
    cout<<diff(V1,V2)<<endl;

    renvoit
    [2,1,1,3]
    [1,4]
    [1,2,3]
     */

    sort(A.begin(), A.end());
    sort(B.begin(), B.end());
    // allocate the smallest size of A,B as maximum size
    vector<int> c(A.size() < B.size() ? B.size() : A.size());
    vector<int>::iterator i;
    i = set_difference(A.begin(), A.end(), B.begin(), B.end(), c.begin());
    return vector<int>(c.begin(), i);
}

//===========================
vector<int> Intersection(vector<int> A,vector<int> B)
{
    /*
    Renvoit elements presents dans les deux vecteurs (avec multiplicité)

    Exemple

    vector<int> V1; V1.push_back(2); V1.push_back(1);  V1.push_back(3);
    vector<int> V2; V2.push_back(1); V2.push_back(1);  V2.push_back(4);
    cout<<V1<<endl;
    cout<<V2<<endl;
    cout<<Intersection(V1,V2)<<endl;

    renvoit
    [2,1,3]
    [1,1,4]
    [1]
     */

    sort(A.begin(), A.end());
    sort(B.begin(), B.end());
    // allocate the smallest size of A,B as maximum size
    vector<int> c(A.size() < B.size() ? B.size() : A.size());
    vector<int>::iterator i;
    i = set_intersection(A.begin(), A.end(), B.begin(), B.end(), c.begin());
    return vector<int>(c.begin(), i);
}

//===========================
vector<int> Union(vector<int> A,vector<int> B)
{
    /*
    Renvoit union des elements des  deux vecteurs

    Exemple

    vector<int> V1; V1.push_back(2); V1.push_back(1);  V1.push_back(3);
    vector<int> V2; V2.push_back(1); V2.push_back(1);  V2.push_back(4);
    cout<<V1<<endl;
    cout<<V2<<endl;
    cout<<Union(V1,V2)<<endl;

    renvoit
    [2,1,3]
    [1,1,4]
    [1,1,2,3,4]
     */

    sort(A.begin(), A.end());
    sort(B.begin(), B.end());
    // allocate the smallest size of A,B as maximum size
    vector<int> c(A.size() < B.size() ? B.size() : A.size());
    vector<int>::iterator i;
    i = set_union(A.begin(), A.end(), B.begin(), B.end(), c.begin());
    return vector<int>(c.begin(), i);

}

//===========================
// renvoit true si j est pas contenu dans V, false si contenu
int Not_find_j_in_V(int j,vector<int> V)
{
    return (find(V.begin(), V.end(), j) == V.end());
}


//===================================


#include <iostream>
#include <iomanip>
#include <sstream>

#include <TLatex.h>
#include <TSystem.h>
#include <TLatex.h>
#include <TArc.h>
#include <TF1.h>
#include <TArrow.h>


double x_v(0),y_v(0); // coordonnees en memoire
int color_v(1); // couleur
int width_v(1); // largeur de ligne
int style_v(1); // style de ligne

//------------------------
void move2(double x,double y)
{
    x_v=x;
    y_v=y;
}
//------------------------
void draw2(double x,double y)
{
    TLine *l= new TLine(x_v,y_v,x,y);
    x_v=x;
    y_v=y;
    l->SetLineColor(color_v);
    l->SetLineWidth(width_v);
    l->SetLineStyle(style_v);
    l->Draw();

// objets_v.Add(l); // rajoute aux objets
//objets_v[objets_v.GetLast()]->Draw(); // dessine
}
//------------------------
void draw_arrow2(double x,double y)
{
    TArrow *l= new TArrow(x_v,y_v,x,y);
    x_v=x;
    y_v=y;
    l->SetLineColor(color_v);
    l->SetLineWidth(width_v);
    l->SetLineStyle(style_v);
    l->SetArrowSize(0.02);
    l->Draw("|>");

// objets_v.Add(l); // rajoute aux objets
//objets_v[objets_v.GetLast()]->Draw(); // dessine
}
//-----------------
void draw2(double x,double y,int col,int width,int style)
{
    color_v=col;
    width_v=width;
    style_v=style;
    draw2(x,y);
}

//---------------------------------
void drawline(double x1,double y1,double x2,double y2,int col,int width,int style)
{

    TLine *l= new TLine(x1,y1,x2,y2);
    x_v=x2;
    y_v=y2;
    l->SetLineColor(col);
    l->SetLineWidth(width);
    l->SetLineStyle(style);
    l->Draw();
}

void drawtext(double x,double y, const char * s, double size)
{
    TLatex *text=new TLatex();
    text->SetTextSize(size);
    text->DrawLatex(x,y,s);
}

//-----------------------
void color(int c)
{
    color_v=c;
}

//=================================
TArc * circle(double x,double y, double r,int col, double phi1, double phi2, int FillStyle)
{
    TArc *e=new TArc(x,y,r,phi1,phi2);

    e->SetLineColor(col);
    e->SetFillStyle(FillStyle);
    e->SetFillColor(col);
    e->Draw("only"); // pour ne pas tracer les rayons.
	return e;
}
//=========================================
// dessin arc de cercle (angles en radian)par segments de droites. sens=+1: direct
void Dessin_arc(double xc,double yc,double R,double a1,double a2,int sens)
{
    if(sens==-1)
    {
        double at=a2;
        a2=a1;
        a1=at;
    }


    while(a2<a1)
        a2+=2*M_PI;

    while(a2>a1+2*M_PI)
        a2-=2*M_PI;


    double x0,y0,x1,y1;
    TLine *l;
    x0=xc+R*cos(a1);
    y0=yc+R*sin(a1);

    double  a=a1;

    while (a<a2)
    {
        a=a+0.05;
        if(a>a2) a=a2;
        x1=xc+R*cos(a);
        y1=yc+R*sin(a);

        l= new TLine(x0,y0,x1,y1);
        l->SetLineColor(BLACK);
        l->Draw();
        x0=x1;
        y0=y1;
    }




}
//========================================================
// Dessin point (petit cercle: style=6)
TMarker * point(double x,double y,  int col,int style,int size)
{
    TMarker *pp=new TMarker(x,y,style); // style=6:petit cercle
    pp->SetMarkerColor(col);
    pp->SetMarkerSize(size);
    pp->Draw();
    x_v=x;
    y_v=y;
    return pp;
}

//======Assigne x_v,y_v ===============
void Set_xv_yv(double xv,double yv)
{
    x_v=xv;
    y_v=yv;
}
//================================================
/* Sauve fichier titre.gif dans repertoire rep (defaut ="."), de la fenetre c en mode batch
 */
void Sauve_gif(TPad *c,const char * titre,const char *rep,int xsize,int ysize)
{
    /*
    ostringstream chaine1,chaine2,chaine3,chaine4,chaine5,chaine6,chaine7;

    chaine1<<rep<<"/"<<titre<<".eps"<<ends;   // chaine1="rep/titre.eps"
    //chaine2<<"pstopnm  -portrait -xsize  "<<xsize<<"  -ysize "<<ysize<<" "<<chaine1.str()<<";"<<ends;
    chaine2<<"pstopnm  -portrait "<<chaine1.str().c_str()<<";"<<ends;
    chaine3<<rep<<"/"<<titre<<".eps001.ppm"<<ends;   // chaine3="rep/titre.eps001.ppm"
    chaine4<<rep<<"/"<<titre<<".gif"<<ends;   // chaine4="rep/titre.gif"
    chaine5<<" ppmtogif "<<chaine3.str()<<" > "<<chaine4.str().c_str()<<";"<<ends;
    chaine6<<"rm -f "<<chaine3.str().c_str()<<";"<<ends;
    chaine7<<" rm -f "<<chaine1.str().c_str()<<";"<<ends;


    c->SaveAs(chaine1.str().c_str());
    gSystem->Exec(chaine2.str().c_str());
    gSystem->Exec(chaine5.str().c_str());
        gSystem->Exec(chaine6.str().c_str());
        gSystem->Exec(chaine7.str().c_str());
    */
    ostringstream chaine4;
    chaine4<<rep<<"/"<<titre<<".gif"<<ends;   // chaine4="rep/titre.gif"

    c->SaveAs(chaine4.str().c_str());

}

//============================================
/* Crée fichier d'animation nom_fichier.gif dans repertoire rep="." (par default).
   (si i<imax, une image gif est crée. Si i=imax, elles sont unifiée pour former un seul fichier d'animation.
  fen: fenetre du dessin
  i: indice courant
  imin,imax: intervalle de i
  nom_fichier: "nom"  nom du fichier final animation
  texte: texte à afficher  sur la fenetre
  rep: repertoire final du fichier
  periode (en 1/100 sec)
*/
void Animation(TPad *fen,int i,int imin,int imax,const char * nom_fichier,int periode,const char * texte, const char * rep)
{
    int opt_jpg=0; //creation du film mp4 à partir des images .jpg


    //.... affiche le texte ....
    fen->cd();
    //  TText text(0.3,0.95,texte); text.SetNDC(kTRUE);
    TLatex text(0.3,0.95,texte);
    text.SetNDC(kTRUE);
    text.Draw();
    //fen->Update();


//..... sauve image png..................



    ostringstream chaine_gif;
    // cela impose 5 chiffres et rajoute des 0 devant si il faut
    chaine_gif<<rep<<"/"<<nom_fichier<<"_"<<setfill ('0') <<setw (5)<<i<<".png"<<ends;
    fen->SaveAs(chaine_gif.str().c_str());

    //.... sauve image jpeg aussi ...

    ostringstream chaine_jpg;
    if(opt_jpg)
    {
        chaine_jpg<<rep<<"/"<<nom_fichier<<"_"<<setfill ('0') <<setw (5)<<i<<".jpg"<<ends;
        fen->SaveAs(chaine_jpg.str().c_str());
    }
//...........................


    if(i==imax)
    {
        //----------- creation du film gif animé à partir des images .gif ------------
        /*
        cout<<"On fait le film.gif.."<<endl;
        ostringstream com;//commande
        com<<"gifsicle  --colors 256 -d="<<periode<<" "; // cadence temporelle
        for(int it=1;it<=(imax-imin+1); it++)
        com<<rep<<"/"<<nom_fichier<<"_"<<setfill ('0') <<setw (5)<<it<<".gif ";
           com<<"  > "<<rep<<"/"<<nom_fichier<<".gif;  gifsicle -b -O2 "<<rep<<"/"<<nom_fichier<<".gif"<<ends;
           gSystem->Exec(com.str().c_str());
        */
        //----------- creation du film gif animé à partir des images png (en convertissant en gif et utilise gifsicle)  ------------
        cout<<"On fait le film.gif.."<<endl;
        ostringstream com;//commande
        com<<"convert "<<nom_fichier<<"_*.png  GIF:- | gifsicle --delay="<<periode<<" --loop --optimize=2 --colors=256 --multifile - > "<<nom_fichier<<".gif"<<ends;

        gSystem->Exec(com.str().c_str());


        ostringstream com2;//commande pour effacers fichiers intermediares
        com2<<"rm -f "<<rep<<"/"<<nom_fichier<<"_*.png "<<ends;
        gSystem->Exec(com2.str().c_str());

        cout<<"Les fichiers animation: "<<rep<<"/"<<nom_fichier<<".gif  ont ete cree. Voir avec xanim."<<endl;

        //----------- creation du film mp4 à partir des images .jpg ------------
        if(opt_jpg)
        {
            int fps=100./periode ;// nbre images /sec.
            ostringstream com3;
            com3<<"    mencoder -mf fps="<<fps<<" \"mf://*.jpg\" -o "<<rep<<"/"<<nom_fichier<<".avi  -ovc lavc -lavcopts vcodec=mjpeg"<<ends;
            //com3<<"mencoder "<<rep<<"/"<<nom_fichier<<".gif  -ovc lavc -o "<<rep<<"/"<<nom_fichier<<".avi"<<ends;
            gSystem->Exec(com3.str().c_str());


            ostringstream com4;//commande pour effacers fichiers intermediares
            com4<<"rm -f "<<rep<<"/"<<nom_fichier<<"_*.jpg "<<ends;
            gSystem->Exec(com4.str().c_str());
            cout<<"Les fichiers animation: "<<rep<<"/"<<nom_fichier<<".avi  ont été créé. Voir avec mplayer."<<endl;
        }

//-------------------------------------------------

    }

}

// ----------
void Recree_Animation(int imin,int imax,const char * nom_fichier,int periode)
{

    string rep=".";
    {
        cout<<"On fait le film.."<<endl;
        ostringstream com;//commande
        com<<"gifsicle  --colors 256 -d="<<periode<<" "; // cadence temporelle

        for(int it=imin; it<=imax; it++)
            com<<rep<<"/"<<nom_fichier<<"_"<<setfill ('0') <<setw (5)<<it<<".gif ";

        com<<"  > "<<rep<<"/"<<nom_fichier<<".gif;  gifsicle -b -O2 ./"<<nom_fichier<<".gif"<<ends;
        gSystem->Exec(com.str().c_str());

        cout<<com.str().c_str()<<endl;


        cout<<"Les fichiers animation: "<<rep<<"/"<<nom_fichier<<".gif a été créé. Voir avec animate ."<<endl;
    }

}


//================================================================
// Pour dessiner une fonction F(x)
double (*P_fonct)(double);

//----------------- y=f(x) ---------------------
double fonction_fx(double * xx,double *par)
{
    double x=xx[0];
    return P_fonct(x);
}
//-------------------------------------------------

void Dessin_Fonction_x(double (*fonct)(double),double xmin,double xmax,const char *titre, const char * opt_dessin)
{
    P_fonct=fonct;

    TCanvas *c=new TCanvas(titre,titre,5);


    TF1 *f1 = new TF1("f1",fonction_fx,xmin,xmax,0); // "type C
    f1->Draw(opt_dessin);
    c->Update();

}



//-------constructeur --------
TH1D_f::TH1D_f(const char* name,const char* title, int nbinsx,double xmin,double xmax):TH1D(name,title, nbinsx,xmin,xmax)
{
    cout<<"Constructeur de TH1D_f"<<endl;
    this->SetStats(0);
    //    numero_s=n_s_i;
}
TH1D_f::TH1D_f(TH1D h)
{
    *this=h;
}
//==========================================
//------evenement souris---------------------
void TH1D_f::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{

    _Event(event,px,py);

}



//-------constructeur --------
TH2F_f::TH2F_f(const char* name,const char* title, int nbinsx,double xmin,double xmax,int nbinsy,double ymin,double ymax):TH2F(name,title, nbinsx,xmin,xmax,nbinsy,ymin,ymax)
{
    cout<<"Constructeur de TH2F_f"<<endl;
    this->SetStats(0);


}



void Event_triade(TH2F_f * h,Int_t event, Int_t px, Int_t py);
void Event_Note(TH2F_f * h,Int_t event, Int_t px, Int_t py);

//==========================================
//------evenement souris---------------------
void TH2F_f::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{

    _Event(event,px,py);


}





//==============================================

//==========================================================
/*... signal( fichier wav) -> vecteur

Entrée: nom du fichier wav
         signal: vecteur pour la sortie
		 opt_select=1: selectionne intervalle de temps (s1,s2) en sec.
                   =0: prend tout le fichier


Sortie SampleRate
 modifie aussi N,Dt,T.
*/
int Lecture_fichier_wav(const char * nom_fichier, vec & signal, int opt_select,double s1,double s2)
{
    cout<<endl<<"=====Crée tableau à partir du fichier: "<<nom_fichier<<endl;

    ifstream pfich(nom_fichier);
    if(!pfich.good())
        cerr<<"Ouverture du fichier:"<<nom_fichier<<" a echouée."<<endl;



//......Lit l'entete du fichier wav.....................................;

    char ChunkID[4];
    pfich.seekg(0,ios::beg);
    pfich.read(ChunkID,4);
    char Format[4];
    pfich.seekg(8,ios::beg);
    pfich.read(Format,(sizeof(char)*4));
    short int NumChannels;
    pfich.seekg(22);
    pfich.read((char*)&NumChannels,sizeof(short int));
    unsigned int SampleRate;
    pfich.seekg(24);
    pfich.read((char*)&SampleRate,sizeof(unsigned int));
    short int BitsPerSample;
    pfich.seekg(34);
    pfich.read((char*)&BitsPerSample,sizeof(short int));
    int NumberOfData;
    pfich.seekg(40);
    pfich.read((char*)&NumberOfData,sizeof(int));
    short int BlockAlign;
    pfich.seekg(32);
    pfich.read((char*)&BlockAlign,sizeof(short int));


    cout<<"ChunkID="<<ChunkID<<endl;     //"RIFF"
    cout<<"format="<<Format<<endl;     //"WAVE"
    cout<<"Nombre de canaux="<<NumChannels<<endl;  // 1 ou 2
    cout<<"Frequence d'echantillonage="<<SampleRate<<" Hz"<<endl;  // 8000 ou 44100 ou...
    cout<<"Bits par echantillon="<<BitsPerSample<<endl;  // 8 ou 16 ou ...
    cout<<"Nombre de données ="<<NumberOfData<<endl;  //
    cout<<"Nombre d'octets par echantillon="<<BlockAlign<<endl; // se deduit de cidessus

//......Deduit les informations..............................

    int NumSamples=NumberOfData/NumChannels/(BitsPerSample/8); // nombre d'echantillons

    int N=NumSamples; //nbre d'echantillons
    double Dt=1.0/SampleRate; // pas en temps (durée entre deux echantillons)
    double T=Dt*NumSamples; // durée totale du signal

    cout<<"Nombre d' echantillons ="<<NumSamples<<endl;  //
    cout<<"Pas de temps dt ="<<Dt<<" s."<<endl;
    cout<<"Durée T ="<<T<<" s."<<endl;


//----------------------------
    if(NumChannels==2) // si stereo
        cout<<"Signal stereo: On ne prendra que le cannal gauche"<<endl;


//........Lecture des données.................................
    int ia=1,ib=N,Nnew(N); //position indice sur fichier

    if(opt_select)
    {
        ia=(int)(s1/T*N);
        ib=(int)(s2/T*N);
        Nnew=ib-ia+1;
    }

    signal.resize(Nnew);

    for(int is=1; is<=Nnew; is++)
    {
        double res(0);
        int pos=(is-1+ia-1);
        if((pos>=0)&&(pos<N)) //lit donnée sur fichier
        {
            pos=pos*BlockAlign+44;  // position en octets depuis le debut
            //place le pointeur a la position pos (en octets)du fichier
            pfich.seekg(pos);

            //lit valeur à cette position
            if(BitsPerSample==8) // unsigned char 0->255
            {
                unsigned char cx;
                pfich.read((char *)&cx,1);
                res=(cx+128);

            }
            else if(BitsPerSample==16) // unsigned char 0->255
            {
                short int  cx;
                pfich.read((char*)&cx,2);
                res=(cx);
            }

        }//if pos
        signal(is-1)=res;
        //signal[is]=signal[is]/1000;
    }// for is
//....... ferme fichier ....................;
    pfich.close();
//	fclose(pfich);

//............Modifie s1,s2,T,N..........;
    if(opt_select)
    {
        s2=s2-s1;
        s1=0;
        T=s2;
        N=Nnew;
    }


    cout<<"Lecture fichier finie."<<endl;
    cout<<"--------------------"<<endl<<endl;

    return SampleRate;
}



//==========================================================
/*... tableau -> signal( fichier wav, 16- bits, Mono)
 */

void Ecriture_fichier_wav(vec & signal,double Dt,const char* nom_fichier)
{
    cout<<endl<<"=====Crée un fichier wav à partir du tableau ====fichier:"<<nom_fichier<<endl;




    //...... ouvre le fichier en ecriture..........................

    ofstream pfich(nom_fichier);
    if (!pfich.is_open())
    {

        cerr<<" erreur ecriture fichier:"<<nom_fichier<<endl;
        return;
    }

    //------------ prépare les données (lire Utils/ondelettes/format_wav.html)------------------
    const char *ChunkID="RIFF";
    unsigned int ChunkSize; // affecté après
    const char *Format="WAVE";
    const char *Subchunk1ID="fmt ";
    unsigned int Subchunk1Size=16;
    short int AudioFormat=1;
    short int NumChannels=1; // :Mono
    unsigned int SampleRate=( unsigned int)(1/Dt);
    unsigned int ByteRate;     // affecté après
    short int BlockAlign; // affecté après
    short int  BitsPerSample=16; // 16- bits
    const char *Subchunk2ID="data";
    unsigned int Subchunk2Size;  // affecté après


    //........................
    int NumSamples=signal.size();
    double T=signal.size() * Dt;


    Subchunk2Size = (unsigned int) (NumSamples * NumChannels * BitsPerSample/8);
    ByteRate  =  (unsigned int )(SampleRate * NumChannels * BitsPerSample/8);
    BlockAlign = ( short int) (NumChannels * BitsPerSample/8);
    ChunkSize=  ( unsigned int)(36 + Subchunk2Size);
    //.........................................



    cout<<"Nombre de canaux="<<NumChannels<<endl;  // 1 ou 2
    cout<<"Frequence d'echantillonage="<<SampleRate<<" Hz"<<endl;  // 8000 ou 44100 ou...
    cout<<"Bits par echantillon="<<BitsPerSample<<endl;  // 8 ou 16 ou ...

    cout<<"Nombre d'octets par echantillon="<<BlockAlign<<endl; // se deduit de cidessus

    cout<<"Nombre d' echantillons ="<<NumSamples<<endl;  //
    cout<<"Pas de temps dt ="<<Dt<<" s."<<endl;
    cout<<"Durée T ="<<T<<" s."<<endl;



//......Ecrit l'entete du fichier wav.....................................;
    pfich.write(ChunkID,4);
    pfich.write((const char*)&ChunkSize,4);
    pfich.write(Format,4);
    pfich.write(Subchunk1ID,4);
    pfich.write((const char*)&Subchunk1Size,4);
    pfich.write((const char*)&AudioFormat,2);
    pfich.write((const char*)&NumChannels,2);
    pfich.write((const char*)&SampleRate,4);
    pfich.write((const char*)&ByteRate,4);
    pfich.write((const char*)&BlockAlign,2);
    pfich.write((const char*)&BitsPerSample,2);
    pfich.write((const char*)Subchunk2ID,4);
    pfich.write((const char*)&Subchunk2Size,4);


//........... normalise le signal ..............

// double Smax=32767; // le maxi voulu
    double Smax=16000; // le maxi voulu
    double smax=0;
    for(int is=0; is<signal.size(); is++) // cherche max
        if(signal(is)>smax)
            smax=signal(is);

    double Coef=Smax/smax;

//...............; Ecrit les données .............


    for(int is=0; is<signal.size(); is++)
    {
        short int  cx=(short int )(signal(is)*Coef);
        pfich.write((const char*)&cx,2);

    }// for is

//....... ferme fichier ....................;
    pfich.close();

}