/* The functions in this file allow to perform vector calculus with */
/* complex vectors */

#include <vector.h>


/* the function Vini() initializes a vector */

vec Vini(dcomplex x, dcomplex y, dcomplex z)
{
     vec t;

     t.v[0].re = Cre(x);
     t.v[0].im = Cim(x);

     t.v[1].re = Cre(y);
     t.v[1].im = Cim(y);

     t.v[2].re = Cre(z);
     t.v[2].im = Cim(z);

     return t;
}


/* the function Vpolar() constructs a complex vector */
/* starting from the magnitude and the polar angles */

vec Vpolar(double magn, double th, double phi)
{
     vec t;

     t.v[0].re = magn * sin(th) * cos(phi);
     t.v[0].im = 0.0;

     t.v[1].re = magn * sin(th) * sin(phi);
     t.v[1].im = 0.0;

     t.v[2].re = magn * cos(th);
     t.v[2].im = 0.0;

     return t;
}


/* The function polar_angle() returns the polar angle(in radians)
/* of a real vector. */
/* Range of polar angle : [0, PI]. */

double polar_angle(vec t)
{
     return acos( (t.v[2].re) / Vnorm(t) );
}


/* The function azimut_angle() returns the azimuthal angle(in radians)
/* of a real vector. */
/* Range of azimuthal angle : [0, 2 * PI]. */

double azimut_angle(vec t)
{
     if (t.v[0].re == 0.)
     {
	if (t.v[1].re >= 0.)
	   return .5 * PI;
	else
	   return 1.5 * PI;
     }
     else
     {
	if (t.v[0].re > 0.)
	{
	   if(t.v[1].re >= 0.)
	      return atan(t.v[1].re / t.v[0].re);
	   else
	      return (2 * PI - atan(fabs(t.v[1].re / t.v[0].re)) );
	}
	else
	{
	   if(t.v[1].re >= 0.)
	      return (PI - atan(fabs(t.v[1].re / t.v[0].re)) );
	   else
	      return (PI + atan(fabs(t.v[1].re / t.v[0].re)) );
	}
     }
}


/* the function Vunit() returns the cartesian unit vectors */

vec Vunit( int i)
{
     vec t;

     switch (i)
     {
     case 1 :
	    {
	    t = Vini( Complex(1.0, 0.0),
		      Complex(0.0, 0.0),
		      Complex(0.0, 0.0) );
	    return t;
	    break;
	    }
     case 2 :
	    {
	    t = Vini( Complex(0.0, 0.0),
		      Complex(1.0, 0.0),
		      Complex(0.0, 0.0) );
	    return t;
	    break;
	    }
     case 3 :
	    {
	    t = Vini(Complex(0.0, 0.0),
                     Complex(0.0, 0.0),
                     Complex(1.0, 0.0) );
	    return t;
	    break;
	    }
     }

}

/* the function Vspherunit() returns the spherical unit vectors */

vec Vspherunit( int lam)
{
     vec t;

     switch (lam)
     {
     case 1  :
	     {
	     t = Vini( Complex( - 1./sqrt(2.), 0.0 ),
		       Complex(0.0 , - 1./sqrt(2.) ),
		       Complex(0.0, 0.0) );
	     return t;
	     break;
	     }
     case -1 :
	     {
	     t = Vini( Complex( + 1./sqrt(2.), 0.0 ),
		       Complex(0.0 , - 1./sqrt(2.) ),
		       Complex(0.0,0.0) );
	     return t;
	     break;
	     }
     case 0  :
	     {
	     t = Vini(Complex(0.0,0.0),Complex(0.0,0.0),Complex(1.0,0.0) );
	     return t;
	     break;
	     }
     }

}


/* the function Vspher() returns the spherical components of a vector */

dcomplex Vspher(vec a, int lam)
{
     return Vscalar( a, Vspherunit(lam) );

}


/* the function Vecho() echoes a vector to the screen */

void Vecho( vec a)
{
     printf("[");Cecho(a.v[0]);
     printf(",  \n");Cecho(a.v[1]);
     printf(",  \n");Cecho(a.v[2]);
     printf("]");
}


/* the function Vadd() adds two vectors */

vec Vadd(vec a, vec b)
{
     vec t;
     int i;

     for(i = 0; i < 3; i++)
     {
	t.v[i] = Cadd(a.v[i], b.v[i]);
     }
     return t;
}


/* the function Vsub() subtracts two vectors */

vec Vsub(vec a, vec b)
{
     vec t;
     int i;

     for(i = 0; i < 3; i++)
     {
	t.v[i] = Csub(a.v[i], b.v[i]);
     }
     return t;
}


/* the function Vconjg() returns the complex conjugate of a vector */

vec Vconjg( vec a)
{
     vec t;
     int i;

     for (i = 0; i < 3; i++)
     {
	t.v[i] = Conjg(a.v[i]);
     }
     return t;
}


/* the function Vscalar() return the scalar product of two vectors */

dcomplex Vscalar( vec a, vec b)
{
     dcomplex s;
     int i;

     s = Complex(0.0,0.0);
     for(i = 0; i < 3; i++)
     {
	s = Cadd( s, Cmul(a.v[i], b.v[i]) );
     }
     return s;
}


/* the function Vnorm() return the norm of a vector */

double Vnorm( vec a)
{
      return sqrt( Cre( Vscalar(a, Vconjg(a)) ) );
}


/* the function Vcross() returns the cross product of two vectors */

vec Vcross( vec a, vec b)
{
     vec t;

     t.v[0] = Csub(Cmul(a.v[1],b.v[2]),Cmul(a.v[2],b.v[1]));
     t.v[1] = Csub(Cmul(a.v[2],b.v[0]),Cmul(a.v[0],b.v[2]));
     t.v[2] = Csub(Cmul(a.v[0],b.v[1]),Cmul(a.v[1],b.v[0]));

     return t;
}


/* the function VRmul() gives the multiplication of a vector with a 
   real number */

vec VRmul( vec a, double b)
{
     vec t;

     t.v[0] = CRmul(b, a.v[0]);
     t.v[1] = CRmul(b, a.v[1]);
     t.v[2] = CRmul(b, a.v[2]);

     return t;
}

/* the function VCmul() gives the multiplication of a vector with a 
   complex number */

vec VCmul( vec a, dcomplex b)
{
     vec t;

     t.v[0] = Cmul(a.v[0],b);
     t.v[1] = Cmul(a.v[1],b);
     t.v[2] = Cmul(a.v[2],b);

     return t;
}


/* the function dsign() returns the sign of a double */

double dsign(double a)
{
  if (a < 0.) return -1.;
  else return 1.;
}


/* The function pauli() returns a vector */
/* consisting of the three pauli matrices. */
/* The row and column indices take the values + 1/2 or - 1/2  . */
/* The 4 different elements form a pauli matrix as follows :
/*
		   |  ++  +-  |
		   |          |
		   |  -+  --  |

   < b | s_x | a > = 1 - d(ab)
   < b | s_y | a > = i sign(a) (1 - d(ab))
   < b | s_z | a > = sign(a) d(ab)

   where d(ab) is the delta function of ab i.e. d(ab) = 1 if a == b */

vec pauli(double row, double col)
{
     vec t;
     double delta;

     delta = fabs(row + col);

     t.v[0] = Complex(1. - delta,0.);
     t.v[1] = Complex(0., (dsign(col)) * (1. - delta));
     t.v[2] = Complex( (dsign(col)) * delta,0.);

     return t;
}  


/* the function unit() returns the matrixelements 
   of the 2 x 2 unit matrix */

double unitm(double s1, double s2)
{
     if (s1 == s2)
	return 1.0;
     else
	return 0.0;
}


/* The function antisymm3_epsilon() returns the symbol epsilon  
   antisymmetric in 3 indices */
   
double antisymm3_epsilon(int i, int j, int k)
{
     double eps;
  
     if( (i == j) || (i == k) || (j == k) )
     {
       eps = 0.;	
     }
     else 
     {
       switch(i)
       {
       case 1:
             {
             switch(j)              
             {
	     case 2:
                   {
                   if (k == 3)
                      eps = +1; 
                      
                   break;
                   }
	     case 3:
                   {
                   if (k == 2)
                      eps = -1; 

                   break;
                   }
             }  
	     break;
	     }        
       case 2:
             {
             switch(j)              
             {
	     case 1:
                   {
                   if (k == 3)
                      eps = -1; 
                      
                   break;
                   }
	     case 3:
                   {
                   if (k == 1)
                      eps = +1; 

                   break;
                   }
             }  
	     break;
	     }        
       case 3:
             {
             switch(j)              
             {
	     case 1:
                   {
                   if (k == 2)
                      eps = +1; 
                      
                   break;
                   }
	     case 2:
                   {
                   if (k == 1)
                      eps = -1; 

                   break;
                   }
             }  
	     break;
	     }        

       }
     }
     
     return eps;
}