flyopunta/src/src_uti/gcomplx.cpp

/***************************************************************************
                          gcomplx.cpp  -  description
                             -------------------
    begin                : Wed Aug 15 2001
    copyright            : (C) 2001 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *

***************************************************************************/
#define GCOMPLX
#include "parm.h"
#include "gcomplx.h"

   Gcomplx Rgc[100];
   Gcomplx *Prgc;
   int Ireg=0;


// -------


 Gcomplx::Gcomplx(double parm)
  { real=parm;
    imm=parm;
    norma=normaq=0.0;
  }

 Gcomplx::Gcomplx(double a, double b)
  {
     real=a;
     imm=b;
     Norma();
  }

 Gcomplx::~Gcomplx(){}
//-----------------------------------


 double  Gcomplx::Norma()
 { normaq=real*real + imm*imm;
   norma=sqrt(normaq);

   return norma;
 }

  double& Gcomplx::operator[](int i) {
          if(i==0)return real; else return imm;}

// asssegnamento
  Gcomplx & Gcomplx::operator=(Gcomplx& a)                //  ritorna il vettore a..
    { real=a.real;imm=a.imm;     //   si puo' scrivere A=B=C..
      normaq=a.normaq;norma=a.norma; return a;
    }

  Gcomplx & Gcomplx::operator=(Gcomplx * a)
    {   *this=*a;  return *a; }

  Gcomplx &  Gcomplx::operator=(double *a)
    {real=a[0];imm=a[1];
     normaq=real*real + imm*imm;
     norma=sqrt(normaq);
     return *this;
    }
  Gcomplx &  Gcomplx::operator=(double a)
    {real=a;imm=0.;
     normaq=a*a;
     norma=fabs(a);
     return *this;
    }




// incrementi

  Gcomplx & Gcomplx::operator+=(Gcomplx& a)  // attenzione non calcola
    { real+=a.real;imm+=a.imm; return *this;}

  Gcomplx & Gcomplx::operator+=(Gcomplx * a)  // attenzione non calcol
    {  *this+=*a; return *this;}

  Gcomplx & Gcomplx::operator+=(double a)  // attenzione non calcol
    {  real+=a; return *this;}


  Gcomplx & Gcomplx::operator-=(Gcomplx * a) // attenzione non calcol
    {  *this-=*a; return *this;}

  Gcomplx & Gcomplx::operator-=(Gcomplx& a)                 // attenzione non calcola no
    { real-=a.real;imm-=a.imm; return *this;}

  Gcomplx & Gcomplx::operator-=(double a)  // attenzione non calcol
    {  real-=a; return *this;}



 // operatore unario di cambio segno

  Gcomplx & Gcomplx::operator-()
     { real=-real;imm=-imm;return *this;}



// prodotto di due complessi

 Gcomplx & Gcomplx::operator*(Gcomplx &A)
   { Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);

     Prgc->real=real*A.real-imm*A.imm;
     Prgc->imm=imm*A.real+real*A.imm;
     Prgc->Norma();
     return *Prgc;
    }

 Gcomplx & Gcomplx::operator*(Gcomplx *A)
  { return *this*(*A);}


 // prodotto per uno scalare

 Gcomplx & Gcomplx::operator*(double b)
   { Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
     Prgc->real=real*b;
     Prgc->imm=imm*b;
     Prgc->norma=norma*fabs(b);
     Prgc->normaq=normaq*b*b;
     return *Prgc;
    }

 Gcomplx & operator*(double c,Gcomplx & A)  // friend function for symmetry
    { return A*c;}

 Gcomplx & Gcomplx::operator*=(Gcomplx &A)
   {
     double re=real*A.real-imm*A.imm;
     double im =imm*A.real+real*A.imm;
     real=re; imm=im;
     Norma();
     return *this;
   }






//rapporto
 Gcomplx & Gcomplx::operator/(Gcomplx &A)
   {
     return *this * A.Inverso();
    }
 Gcomplx & Gcomplx::operator/(Gcomplx *A)
   { return *this / *A;     }

 Gcomplx & Gcomplx::operator/(double a)
   {
    a=1./a;
    return   *this*a;
   }
 Gcomplx & operator/(double c,Gcomplx & A)  //friend
    { return A.Inverso()*c;}

 Gcomplx & Gcomplx::operator/=(Gcomplx &A)
   {
     *this= *this * A.Inverso();
     Norma();
     return *this ;
    }


 // prodotto scalare
 double Gcomplx::operator%(Gcomplx & A)
   {  return(real*A.real-imm*A.imm);}

 double Gcomplx::operator%(Gcomplx * A)          // attenzione..funziona...........
   { return this[0]%A[0];}

 //somma

  Gcomplx&  Gcomplx::operator+(Gcomplx& A)
    {
      Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
      Prgc->real=real+A.real;
      Prgc->imm=imm  +A.imm;
      return *Prgc;}

  Gcomplx&  Gcomplx::operator+(Gcomplx* a)
     {  return *this + *a;}

  Gcomplx&  Gcomplx::operator+(double r)
     {Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
      Prgc->real=r+real;
      Prgc->imm=imm;
      return *Prgc;}
  Gcomplx&  operator+(double r,Gcomplx& A)
     {Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
      Prgc->real=r+A.real;
      Prgc->imm=A.imm;
      return *Prgc;}


// sottrazione

   Gcomplx&  Gcomplx::operator-(Gcomplx& A)
     {
       Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
       Prgc->real=real-A.real;
       Prgc->imm=imm  -A.imm;
       return *Prgc;}

   Gcomplx&  Gcomplx::operator-(Gcomplx* a)
      {  return *this - *a;}

   Gcomplx&  Gcomplx::operator-(double r)
      {Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
       Prgc->real=real-r;
       Prgc->imm=imm;
       return *Prgc;}
   Gcomplx&  operator-(double r,Gcomplx& A)
      {Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
       Prgc->real=r-A.real;
       Prgc->imm=A.imm;
       return *Prgc;}



 // coniugata

 Gcomplx & Gcomplx::operator&()
     {Ireg++; if(Ireg>99)Ireg=0; Prgc=(Rgc+Ireg);
      Prgc->real=real;
      Prgc->imm=-imm;
      return  *Prgc;}



// inverso
Gcomplx & Gcomplx::Inverso()
  { Ireg++; if(Ireg>99) Ireg=0; Prgc=(Rgc+Ireg);
   Prgc->real=real/normaq;
   Prgc->imm=-imm/normaq;
   Prgc->norma=1./norma;
   Prgc->normaq=1./normaq;
   return  *Prgc;

/*    Rgc[Ireg].real=real/normaq;
   Rgc[Ireg].imm=-imm/normaq;
   Rgc[Ireg].norma=1./norma;
   Rgc[Ireg].normaq=1./normaq;
   return  Rgc[Ireg];
*/

  }

 // ritorna il vettore normalizzato
 Gcomplx & Gcomplx::Verso()
   { Ireg++; if(Ireg>99) Ireg=0; Prgc=(Rgc+Ireg);
     Prgc->real=real/norma;
     Prgc->imm=imm/norma;
     Prgc->norma=Prgc->normaq=1.0;
     return *Prgc;
   }

 Gcomplx & Gcomplx::Set(double a,double b)
   { real=a;
     imm=b;
     normaq=real*real + imm*imm;
     norma=sqrt(normaq);
     return *this;
   }


  // output per debug...............
  void   Gcomplx::print(char *t)
    {
    
     using namespace std;                 
     Gout.precision(5);
     Gout<<std::setfill(' ')<<std::fixed;      
     Gout<<"\n"<<setw(10)<<t<<"=[r: "<<setw(12)<<real<<"  i:=["
     <<setw(12)<<imm<<"] |"<<setw(12)<<norma<<"| "<<endl;       
    }