/* The functions in this file allow to perform four-vector calculations with 
   complex vectors */

#include <four_vector.h>


/* the function V4ini() initializes a contravariant four-vector */

fvec V4ini(dcomplex t, vec r)
{
     fvec a;

     a.fv[0].re = Cre(t);
     a.fv[0].im = Cim(t);

     a.fv[1].re = r.v[0].re;
     a.fv[1].im = r.v[0].im;

     a.fv[2].re = r.v[1].re;
     a.fv[2].im = r.v[1].im;
		  
     a.fv[3].re = r.v[2].re;
     a.fv[3].im = r.v[2].im;

     return a;
}


/* the function V4zero() returns a zero four-vector */

fvec V4zero()
{
     fvec a;

     a.fv[0].re = 0.;
     a.fv[0].im = 0.;
     	  
     a.fv[1].re = 0.;
     a.fv[1].im = 0.;
		  
     a.fv[2].re = 0.;
     a.fv[2].im = 0.;
		  
     a.fv[3].re = 0.;
     a.fv[3].im = 0.;

     return a;
}


/* the function V4echo() echoes a four-vector to the screen */

void V4echo( fvec a)
{
     printf("[");Cecho(a.fv[0]);
     printf(",  \n");Cecho(a.fv[1]);
     printf(",  \n");Cecho(a.fv[2]);
     printf(",  \n");Cecho(a.fv[3]);
     printf("]");
}


/* the function V4add() adds two contravariant four-vectors */

fvec V4add(fvec a, fvec b)
{
     fvec t;
     int mu;

     for(mu = 0; mu < 4; mu++)
     {
	t.fv[mu] = Cadd(a.fv[mu], b.fv[mu]);
     }
     return t;
}


/* the function V4sub() subtracts two contravariant four-vectors */

fvec V4sub(fvec a, fvec b)
{
     fvec t;
     int i;

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


/* the function V4mul() gives the contraction of two four-vectors */

dcomplex V4mul(fvec a, fvec b)
{
     dcomplex c;
     int mu1, mu2;

     c = Complex(0., 0.);
     
     for(mu1 = 0; mu1 <= 3; mu1++)
     {
       for(mu2 = 0; mu2 <= 3; mu2++)
       {
       c = Cadd(c, CRmul(metric(mu1, mu2),
                         Cmul(a.fv[mu1], b.fv[mu2]) ) 
               );
       }
     }
 
     return c;  
}


/* the function V4conjg() returns the complex conjugate of a four-vector */

fvec V4conjg( fvec a)
{
     fvec t;
     int i;

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


/* the function V4Rmul() gives the multiplication of a 
   four-vector with a real number */

fvec V4Rmul( fvec a, double b)
{
     fvec t;

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

     return t;
}

/* the function V4Cmul() gives the multiplication of a four-vector 
   with a complex number */

fvec V4Cmul( fvec a, dcomplex b)
{
     fvec t;

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

     return t;
}


double metric(int mu, int nu)
{
     double g;
  
     if(mu == nu)
     {
        if (mu == 0)
	  g = 1.;
        else
          g = -1.;	
     }
     else g = 0.;
     
     return g;
}  


double antisymm_epsilon(int ka, int la, int mu, int nu)
{
     double g;
  
     if( (ka == la) || (ka == mu) || (ka == nu) || 
         (la == mu) || (la == nu) || (mu == nu) )
     {
       g = 0.;	
     }
     else 
     {
       switch(ka)
       {
       case 0:
             {
             switch(la)              
             {
	     case 1:
                   {
                   if ((mu == 2) && (nu == 3))
                      g = +1; 
                   if ((mu == 3) && (nu == 2))
                      g = -1;
 
                   break;
                   }
	     case 2:
                   {
                   if ((mu == 1) && (nu == 3))
                      g = -1; 
                   if ((mu == 3) && (nu == 1))
                      g = +1;

                   break;
                   }
	     case 3:
                   {
                   if ((mu == 1) && (nu == 2))
                      g = +1; 
                   if ((mu == 2) && (nu == 1))
                      g = -1;

                   break;
                   }
             }  
	     break;
	     }        
       case 1:
             {
             switch(la)              
             {
	     case 0:
                   {
                   if ((mu == 2) && (nu == 3))
                      g = -1; 
                   if ((mu == 3) && (nu == 2))
                      g = +1;
 
                   break;
                   }
	     case 2:
                   {
                   if ((mu == 0) && (nu == 3))
                      g = +1; 
                   if ((mu == 3) && (nu == 0))
                      g = -1;

                   break;
                   }
	     case 3:
                   {
                   if ((mu == 0) && (nu == 2))
                      g = -1; 
                   if ((mu == 2) && (nu == 0))
                      g = +1;

                   break;
                   }
             }  
	     break;
	     }        
       case 2:
             {
             switch(la)              
             {
	     case 0:
                   {
                   if ((mu == 1) && (nu == 3))
                      g = +1; 
                   if ((mu == 3) && (nu == 1))
                      g = -1;
 
                   break;
                   }
	     case 1:
                   {
                   if ((mu == 0) && (nu == 3))
                      g = -1; 
                   if ((mu == 3) && (nu == 0))
                      g = +1;

                   break;
                   }
	     case 3:
                   {
                   if ((mu == 0) && (nu == 1))
                      g = +1; 
                   if ((mu == 1) && (nu == 0))
                      g = -1;

                   break;
                   }
             }  
             
	     break;
	     }        
       case 3:
             {
             switch(la)              
             {
	     case 0:
                   {
                   if ((mu == 1) && (nu == 2))
                      g = -1; 
                   if ((mu == 2) && (nu == 1))
                      g = +1;
 
                   break;
                   }
	     case 1:
                   {
                   if ((mu == 0) && (nu == 2))
                      g = +1; 
                   if ((mu == 2) && (nu == 0))
                      g = -1;

                   break;
                   }
	     case 2:
                   {
                   if ((mu == 0) && (nu == 1))
                      g = -1; 
                   if ((mu == 1) && (nu == 0))
                      g = +1;

                   break;
                   }
             }  
             
	     break;
	     }        
     }
     }
     
     return g;
}  



spinor_mat SMadd(spinor_mat a, spinor_mat b)
{
     spinor_mat c;
     int i, j;

     for(i = 0; i < 4; i++)
     {
        for(j = 0; j < 4; j++)
        {
	   c.sm[i][j] = Cadd(a.sm[i][j], b.sm[i][j]);
	 }
     }
     return c;
}


spinor_mat SMsub(spinor_mat a, spinor_mat b)
{
     spinor_mat c;
     int i, j;

     for(i = 0; i < 4; i++)
     {
        for(j = 0; j < 4; j++)
        {
	   c.sm[i][j] = Csub(a.sm[i][j], b.sm[i][j]);
	 }
     }
     return c;
}


spinor_mat SMRmul(spinor_mat a, double z)
{
     spinor_mat b;
     int i, j;

     for(i = 0; i < 4; i++)
     {
        for(j = 0; j < 4; j++)
        {
	   b.sm[i][j] = CRmul(z, a.sm[i][j]);
	 }
     }
     return b;
}


spinor_mat SMCmul(spinor_mat a, dcomplex z)
{
     spinor_mat b;
     int i, j;

     for(i = 0; i < 4; i++)
     {
        for(j = 0; j < 4; j++)
        {
	   b.sm[i][j] = Cmul(a.sm[i][j], z);
	 }
     }
     return b;
}


spinor_mat SMmul(spinor_mat a, spinor_mat b)
{
     spinor_mat c;
     int i, j;

     for(i = 0; i < 4; i++)
     {
        for(j = 0; j < 4; j++)
        {
	   int k;
	   dcomplex temp = Complex(0., 0.);
 
           for(k = 0; k < 4; k++)
           {	   
           temp = Cadd(temp, Cmul(a.sm[i][k], b.sm[k][j]) );
	   }
	   
           c.sm[i][j] = temp;
	 }
     }
     return c;
}


spinor_mat SMmul3(spinor_mat a, spinor_mat b, spinor_mat c)
{     
     return SMmul(a, SMmul(b,c));
}


spinor_mat SMmul4(spinor_mat a, spinor_mat b, spinor_mat c, spinor_mat d)
{     
     return SMmul( a, SMmul( b, SMmul(c, d) ) );
}


spinor_mat SMmul5(spinor_mat a, spinor_mat b, spinor_mat c, spinor_mat d, spinor_mat e)
{     
     return SMmul( a, SMmul( b, SMmul(c, SMmul(d, e) ) ) );
}


/* The function dirac_gamma() gives the contravariant Dirac matrices 
   as 4 by 4 matrices in spinor space  (convention Bjorken and Drell 1964). */


spinor_mat dirac_gamma(int mu)
{
     spinor_mat gamma;

     if( mu == 0)
     {
     gamma.sm[0][0] = Complex(1., 0.);
     gamma.sm[0][1] = Complex(0., 0.);
     gamma.sm[0][2] = Complex(0., 0.);
     gamma.sm[0][3] = Complex(0., 0.);
     gamma.sm[1][0] = Complex(0., 0.);
     gamma.sm[1][1] = Complex(1., 0.);
     gamma.sm[1][2] = Complex(0., 0.);
     gamma.sm[1][3] = Complex(0., 0.);
     gamma.sm[2][0] = Complex(0., 0.);
     gamma.sm[2][1] = Complex(0., 0.);
     gamma.sm[2][2] = Complex(-1., 0.);
     gamma.sm[2][3] = Complex(0., 0.);
     gamma.sm[3][0] = Complex(0., 0.);
     gamma.sm[3][1] = Complex(0., 0.);
     gamma.sm[3][2] = Complex(0., 0.);
     gamma.sm[3][3] = Complex(-1., 0.);
     }
	   
     if( mu == 1)
     {
     gamma.sm[0][0] = Complex(0., 0.);
     gamma.sm[0][1] = Complex(0., 0.);
     gamma.sm[0][2] = Complex(0., 0.);
     gamma.sm[0][3] = Complex(1., 0.);
     gamma.sm[1][0] = Complex(0., 0.);
     gamma.sm[1][1] = Complex(0., 0.);
     gamma.sm[1][2] = Complex(1., 0.);
     gamma.sm[1][3] = Complex(0., 0.);
     gamma.sm[2][0] = Complex(0., 0.);
     gamma.sm[2][1] = Complex(-1., 0.);
     gamma.sm[2][2] = Complex(0., 0.);
     gamma.sm[2][3] = Complex(0., 0.);
     gamma.sm[3][0] = Complex(-1., 0.);
     gamma.sm[3][1] = Complex(0., 0.);
     gamma.sm[3][2] = Complex(0., 0.);
     gamma.sm[3][3] = Complex(0., 0.);
     }
	   
     if( mu == 2)
     {
     gamma.sm[0][0] = Complex(0., 0.);
     gamma.sm[0][1] = Complex(0., 0.);
     gamma.sm[0][2] = Complex(0., 0.);
     gamma.sm[0][3] = Complex(0., -1.);
     gamma.sm[1][0] = Complex(0., 0.);
     gamma.sm[1][1] = Complex(0., 0.);
     gamma.sm[1][2] = Complex(0., 1.);
     gamma.sm[1][3] = Complex(0., 0.);
     gamma.sm[2][0] = Complex(0., 0.);
     gamma.sm[2][1] = Complex(0., 1.);
     gamma.sm[2][2] = Complex(0., 0.);
     gamma.sm[2][3] = Complex(0., 0.);
     gamma.sm[3][0] = Complex(0., -1.);
     gamma.sm[3][1] = Complex(0., 0.);
     gamma.sm[3][2] = Complex(0., 0.);
     gamma.sm[3][3] = Complex(0., 0.);
     }

     if( mu == 3)
     {
     gamma.sm[0][0] = Complex(0., 0.);
     gamma.sm[0][1] = Complex(0., 0.);
     gamma.sm[0][2] = Complex(1., 0.);
     gamma.sm[0][3] = Complex(0., 0.);
     gamma.sm[1][0] = Complex(0., 0.);
     gamma.sm[1][1] = Complex(0., 0.);
     gamma.sm[1][2] = Complex(0., 0.);
     gamma.sm[1][3] = Complex(-1., 0.);
     gamma.sm[2][0] = Complex(-1., 0.);
     gamma.sm[2][1] = Complex(0., 0.);
     gamma.sm[2][2] = Complex(0., 0.);
     gamma.sm[2][3] = Complex(0., 0.);
     gamma.sm[3][0] = Complex(0., 0.);
     gamma.sm[3][1] = Complex(1., 0.);
     gamma.sm[3][2] = Complex(0., 0.);
     gamma.sm[3][3] = Complex(0., 0.);
     }
     
     return gamma;
}  


spinor_mat dirac_gamma5()
{
     spinor_mat a;

     a.sm[0][0] = Complex(0., 0.);
     a.sm[0][1] = Complex(0., 0.);
     a.sm[0][2] = Complex(1., 0.);
     a.sm[0][3] = Complex(0., 0.);
     a.sm[1][0] = Complex(0., 0.);
     a.sm[1][1] = Complex(0., 0.);
     a.sm[1][2] = Complex(0., 0.);
     a.sm[1][3] = Complex(1., 0.);
     a.sm[2][0] = Complex(1., 0.);
     a.sm[2][1] = Complex(0., 0.);
     a.sm[2][2] = Complex(0., 0.);
     a.sm[2][3] = Complex(0., 0.);
     a.sm[3][0] = Complex(0., 0.);
     a.sm[3][1] = Complex(1., 0.);
     a.sm[3][2] = Complex(0., 0.);
     a.sm[3][3] = Complex(0., 0.);

     return a;       
}


/* The function dirac_sigma gives the contravariant tensor sigma 
   in spinor space (convention Bjorken and Drell 1964). */

spinor_mat dirac_sigma(int mu, int nu)
{
     spinor_mat a;

     a = SMCmul( SMsub( SMmul(dirac_gamma(mu), dirac_gamma(nu)),
                        SMmul(dirac_gamma(nu), dirac_gamma(mu))
                      ),
                 Complex(0.0, 0.5) 
               );

     return a;  
}


spinor_mat unit_spinor_mat()
{
     spinor_mat a;

     a.sm[0][0] = Complex(1., 0.);
     a.sm[0][1] = Complex(0., 0.);
     a.sm[0][2] = Complex(0., 0.);
     a.sm[0][3] = Complex(0., 0.);
     a.sm[1][0] = Complex(0., 0.);
     a.sm[1][1] = Complex(1., 0.);
     a.sm[1][2] = Complex(0., 0.);
     a.sm[1][3] = Complex(0., 0.);
     a.sm[2][0] = Complex(0., 0.);
     a.sm[2][1] = Complex(0., 0.);
     a.sm[2][2] = Complex(1., 0.);
     a.sm[2][3] = Complex(0., 0.);
     a.sm[3][0] = Complex(0., 0.);
     a.sm[3][1] = Complex(0., 0.);
     a.sm[3][2] = Complex(0., 0.);
     a.sm[3][3] = Complex(1., 0.);

     return a;       
}


spinor_mat zero_spinor_mat()
{
     spinor_mat a;

     a.sm[0][0] = Complex(0., 0.);
     a.sm[0][1] = Complex(0., 0.);
     a.sm[0][2] = Complex(0., 0.);
     a.sm[0][3] = Complex(0., 0.);
     a.sm[1][0] = Complex(0., 0.);
     a.sm[1][1] = Complex(0., 0.);
     a.sm[1][2] = Complex(0., 0.);
     a.sm[1][3] = Complex(0., 0.);
     a.sm[2][0] = Complex(0., 0.);
     a.sm[2][1] = Complex(0., 0.);
     a.sm[2][2] = Complex(0., 0.);
     a.sm[2][3] = Complex(0., 0.);
     a.sm[3][0] = Complex(0., 0.);
     a.sm[3][1] = Complex(0., 0.);
     a.sm[3][2] = Complex(0., 0.);
     a.sm[3][3] = Complex(0., 0.);

     return a;       
}


/* The function fvec_slash gives the contraction of a four-vector a
   with the dirac gamma matrices (convention Bjorken and Drell 1964). */

spinor_mat fvec_slash(fvec a)
{
     spinor_mat b;
     int mu1, mu2;

     b = zero_spinor_mat();     

     for(mu1 = 0; mu1 <= 3; mu1++)
     {
       for(mu2 = 0; mu2 <= 3; mu2++)
       {
       b = SMadd(b, SMCmul(dirac_gamma(mu1), 
                           CRmul(metric(mu1, mu2), a.fv[mu2]) ) 
                );
       }
     }
 
     return b;  
}



/* The function spinor_pos() gives the positive-energy Dirac spinor
   for a particle with momentum "p", mass "mass" 
   and spin projection "sproj" : in the restframe of the particle this 
   is an eigenstate of the z-component of the spin. */ 

spinor spinor_pos(vec p, double mass, double sproj)
{
     double en;
     spinor a;

     en = sqrt( pow(Vnorm(p), 2.) + pow(mass, 2.) );
     
     if ( sproj == 0.5)
     {
     a.s[0] = Complex(1., 0.);
     a.s[1] = Complex(0., 0.);
     a.s[2] = CRmul(1. / (en + mass),
		    Vscalar(p, pauli(0.5, 0.5))
                   );
     a.s[3] = CRmul(1. / (en + mass),
		    Vscalar(p, pauli(- 0.5, 0.5))
                   );
     }

     if ( sproj == - 0.5)
     {
     a.s[0] = Complex(0., 0.);
     a.s[1] = Complex(1., 0.);
     a.s[2] = CRmul(1. / (en + mass),
                    Vscalar(p, pauli(0.5, - 0.5))
		   );
     a.s[3] = CRmul(1. / (en + mass),
                    Vscalar(p, pauli(- 0.5, - 0.5))
		   );
     }
      
     return a;
}


/* The function spinor_pos_xaxis() gives the positive-energy Dirac spinor
   for a particle with momentum "p", mass "mass" 
   and spin projection "sproj" in the x direction. */ 

spinor spinor_pos_xaxis(vec p, double mass, double sproj)
{
     double en;
     spinor a;

     en = sqrt( pow(Vnorm(p), 2.) + pow(mass, 2.) );
     
     if ( sproj == 0.5)
     {
     a.s[0] = Complex(1./sqrt(2.), 0.);
     a.s[1] = Complex(1./sqrt(2.), 0.);
     a.s[2] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(p.v[0], Cmul(p.v[1], Complex(0., - 1.)) ),
		         p.v[2]
		        ) 
                   );
     a.s[3] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Csub(Cadd(p.v[0], Cmul(p.v[1], Complex(0., 1.)) ),
		         p.v[2]
		        ) 
                   );
     }

     if ( sproj == - 0.5)
     {
     a.s[0] = Complex(1./sqrt(2.), 0.);
     a.s[1] = Complex(- 1./sqrt(2.), 0.);
     a.s[2] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(CRmul(-1.,p.v[0]), Cmul(p.v[1],Complex(0., 1.)) ),
		         p.v[2]
		        ) 
                   );
     a.s[3] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(p.v[0], Cmul(p.v[1], Complex(0., 1.)) ),
		         p.v[2]
		        ) 
                   );
     }
      
     return a;
}


/* The function spinor_pos_yaxis() gives the positive-energy Dirac spinor
   for a particle with momentum "p", mass "mass" 
   and spin projection "sproj" in the y direction. */ 

spinor spinor_pos_yaxis(vec p, double mass, double sproj)
{
     double en;
     spinor a;

     en = sqrt( pow(Vnorm(p), 2.) + pow(mass, 2.) );
     
     if ( sproj == 0.5)
     {
     a.s[0] = Complex(1./sqrt(2.), 0.);
     a.s[1] = Complex(0., 1./sqrt(2.));
     a.s[2] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(p.v[1], Cmul(p.v[0], Complex(0., 1.)) ),
		         p.v[2]
		        ) 
                   );
     a.s[3] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(p.v[0], Cmul(p.v[1], Complex(0., 1.)) ),
		         Cmul(p.v[2], Complex(0., -1.) )
		        ) 
                   );
     }

     if ( sproj == - 0.5)
     {
     a.s[0] = Complex(1./sqrt(2.), 0.);
     a.s[1] = Complex(0., - 1./sqrt(2.));
     a.s[2] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(CRmul(-1.,p.v[1]), Cmul(p.v[0],Complex(0.,-1.)) ),
		         p.v[2]
		        ) 
                   );
     a.s[3] = CRmul(1. / (en + mass) * 1./sqrt(2.),
		    Cadd(Cadd(p.v[0], Cmul(p.v[1], Complex(0., 1.)) ),
		         Cmul(p.v[2], Complex(0., 1.) )
		        ) 
                   );
     }
      
     return a;
}


/* The function spinor_hel() gives the helicity spinor
   for a particle with momentum "p", mass 0 
   and helicity "hel" (+/- 1.). */ 

spinor spinor_hel(vec p, double hel)
{
     double en, th, phi;
     spinor a;

     th = polar_angle(p);
     phi = azimut_angle(p);
     en = Vnorm(p);
     
     if ( hel == +1.)
     {
     a.s[0] = Complex(cos(th/2.), 0.);
     a.s[1] = Complex(sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin(phi) );
     a.s[2] = Complex(cos(th/2.), 0.);
     a.s[3] = Complex(sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin(phi) );
     }


     if ( hel == -1.)
     {
     a.s[0] = Complex(- sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin (phi));
     a.s[1] = Complex(cos(th/2.), 0.0);
     a.s[2] = Complex(sin(th/2.) * cos(phi), 
                      - sin(th/2.) * sin (phi));
     a.s[3] = Complex(- cos(th/2.), 0.0);
     }

      
     return a;
}

/* The function spinor_anti_hel() gives the helicity spinor
   for an anti-particle with momentum "p", mass 0 
   and helicity "hel" (+/- 1.). */ 

spinor spinor_anti_hel(vec p, double hel)
{
     double en, th, phi;
     spinor a;

     th = polar_angle(p);
     phi = azimut_angle(p);
     en = Vnorm(p);
     
     if ( hel == +1.)
     {
     a.s[0] = Complex(sin(th/2.) * cos(phi), 
                      -sin(th/2.) * sin(phi) );
     a.s[1] = Complex(-cos(th/2.), 0.);
     a.s[2] = Complex(-sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin(phi) );
     a.s[3] = Complex(cos(th/2.), 0.);
     }


     if ( hel == -1.)
     {
     a.s[0] = Complex(- cos(th/2.), 0.0);
     a.s[1] = Complex(- sin(th/2.) * cos(phi), 
                      - sin(th/2.) * sin (phi));
     a.s[2] = Complex(- cos(th/2.), 0.0);
     a.s[3] = Complex(- sin(th/2.) * cos(phi), 
                      - sin(th/2.) * sin (phi));
     }

      
     return a;
}

/* The function spinor_part_hel() gives the helicity spinor
   for a spin 1/2 particle with momentum "p" 
   mass "mass" and helicity "hel" (+/- .5). */ 

spinor spinor_part_hel(vec p, double mass, double hel)
{
     double pmom, en, th, phi;
     spinor a;

     th = polar_angle(p);
     phi = azimut_angle(p);

     if(th > 3.14159265)
       phi = 0.;
      
     pmom = Vnorm(p);
     en = sqrt(pow(pmom,2.) + pow(mass,2.));
     
     if ( hel == .5)
     {
     a.s[0] = Complex(cos(th/2.), 0.);
     a.s[1] = Complex(sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin(phi) );
     a.s[2] = Complex(pmom/(en + mass) * cos(th/2.), 0.);
     a.s[3] = Complex(pmom/(en + mass) * sin(th/2.) * cos(phi), 
                      pmom/(en + mass) * sin(th/2.) * sin(phi) );
     }


     if ( hel == -.5)
     {
     a.s[0] = Complex(- sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin (phi));
     a.s[1] = Complex(cos(th/2.), 0.0);
     a.s[2] = Complex(pmom/(en + mass) * sin(th/2.) * cos(phi), 
                      - pmom/(en + mass) * sin(th/2.) * sin (phi));
     a.s[3] = Complex(- pmom/(en + mass) * cos(th/2.), 0.0);
     }
      
     return a;
}

/* The function spinor_antipart_hel() gives the helicity spinor
   for a spin 1/2 anti-particle with momentum "p" 
   mass "mass" and helicity "hel" (+/- .5). */ 

spinor spinor_antipart_hel(vec p, double mass, double hel)
{
     double pmom, en, th, phi;
     spinor a;

     th = polar_angle(p);
     phi = azimut_angle(p);

     if(th > 3.14159265)
       phi = 0.;

     pmom = Vnorm(p);
     en = sqrt(pow(pmom,2.) + pow(mass,2.));
     
     if ( hel == .5)
     {
     a.s[0] = Complex(pmom/(en + mass) * sin(th/2.) * cos(phi), 
                      - pmom/(en + mass) * sin(th/2.) * sin (phi));
     a.s[1] = Complex(- pmom/(en + mass) * cos(th/2.), 0.0);
     a.s[2] = Complex(- sin(th/2.) * cos(phi), 
                      sin(th/2.) * sin (phi));
     a.s[3] = Complex(cos(th/2.), 0.0);
     }

     if ( hel == -.5)
     {
     a.s[0] = Complex(- pmom/(en + mass) * cos(th/2.), 0.);
     a.s[1] = Complex(- pmom/(en + mass) * sin(th/2.) * cos(phi), 
                      - pmom/(en + mass) * sin(th/2.) * sin(phi) );
     a.s[2] = Complex(- cos(th/2.), 0.);
     a.s[3] = Complex(- sin(th/2.) * cos(phi), 
                      - sin(th/2.) * sin(phi) );
     }
      
     return a;
}


/* The function spinor_adj() gives the adjoint Dirac spinor */

spinor spinor_adj(spinor a)
{
     spinor b;

     b.s[0] = Conjg(a.s[0]); 
     b.s[1] = Conjg(a.s[1]);
     b.s[2] = CRmul(- 1., Conjg(a.s[2]));
     b.s[3] = CRmul(- 1., Conjg(a.s[3]));

     return b;
     
}


dcomplex spinleft_mat_spinright(spinor l, spinor_mat m, spinor r)
{
     int i, j;
     dcomplex sum = Complex(0., 0.);
     spinor l_adj;
      
     l_adj = spinor_adj(l);
     
     for(i = 0; i <= 3; i++)
     {
        for(j = 0; j <= 3; j++)
        {
	   sum = Cadd(sum, Cmul(l_adj.s[i], Cmul(m.sm[i][j], r.s[j])) );
	}
     }	
     
     return sum;
     
}