mat_SU_N.cpp

Go to the documentation of this file.

#include "mat_SU_N.h"
using namespace SU_N;

#include <cassert>

#include "randomNumbers.h"

//==========================================================
Mat_SU_N& Mat_SU_N::reunit()
{
  // This implementation is applicable only to SU(3).
  //                             H.Matsufuru [5 Feb 2012]

  assert(m_Nc == 3);

  Mat_SU_N rhs(m_Nc);
  for (int i = 0; i < 2 * m_Nc * m_Nc; ++i) {
    rhs.va[i] = va[i];
  }

  double sn1 = 1.0 / (rhs.va[0] * rhs.va[0] + rhs.va[1] * rhs.va[1]
                      + rhs.va[2] * rhs.va[2] + rhs.va[3] * rhs.va[3]
                      + rhs.va[4] * rhs.va[4] + rhs.va[5] * rhs.va[5]);
  sn1 = sqrt(sn1);

  va[0] = rhs.va[0] * sn1;
  va[1] = rhs.va[1] * sn1;
  va[2] = rhs.va[2] * sn1;
  va[3] = rhs.va[3] * sn1;
  va[4] = rhs.va[4] * sn1;
  va[5] = rhs.va[5] * sn1;

  double sp1r = va[0] * rhs.va[6] + va[1] * rhs.va[7]
                + va[2] * rhs.va[8] + va[3] * rhs.va[9]
                + va[4] * rhs.va[10] + va[5] * rhs.va[11];
  double sp1i = va[0] * rhs.va[7] - va[1] * rhs.va[6]
                + va[2] * rhs.va[9] - va[3] * rhs.va[8]
                + va[4] * rhs.va[11] - va[5] * rhs.va[10];

  va[6]  = rhs.va[6] - sp1r * va[0] + sp1i * va[1];
  va[7]  = rhs.va[7] - sp1r * va[1] - sp1i * va[0];
  va[8]  = rhs.va[8] - sp1r * va[2] + sp1i * va[3];
  va[9]  = rhs.va[9] - sp1r * va[3] - sp1i * va[2];
  va[10] = rhs.va[10] - sp1r * va[4] + sp1i * va[5];
  va[11] = rhs.va[11] - sp1r * va[5] - sp1i * va[4];

  double sn2 = 1.0 / (va[6] * va[6] + va[7] * va[7]
                      + va[8] * va[8] + va[9] * va[9]
                      + va[10] * va[10] + va[11] * va[11]);
  sn2 = sqrt(sn2);

  va[6]  = va[6] * sn2;
  va[7]  = va[7] * sn2;
  va[8]  = va[8] * sn2;
  va[9]  = va[9] * sn2;
  va[10] = va[10] * sn2;
  va[11] = va[11] * sn2;

  va[12] = va[2] * va[10] - va[3] * va[11]
           - va[4] * va[8] + va[5] * va[9];

  va[13] = -va[2] * va[11] - va[3] * va[10]
           + va[4] * va[9] + va[5] * va[8];

  va[14] = va[4] * va[6] - va[5] * va[7]
           - va[0] * va[10] + va[1] * va[11];

  va[15] = -va[4] * va[7] - va[5] * va[6]
           + va[0] * va[11] + va[1] * va[10];

  va[16] = va[0] * va[8] - va[1] * va[9]
           - va[2] * va[6] + va[3] * va[7];

  va[17] = -va[0] * va[9] - va[1] * va[8]
           + va[2] * va[7] + va[3] * va[6];

  return *this;
}


//==========================================================
Mat_SU_N& Mat_SU_N::set_random(RandomNumbers *rand)
{
  // temporary implementation: only applicable to SU(3).
  //                            H.Matsufuru [5 Feb 2012]

  assert(m_Nc == 3);

  static const double PI  = 4.0 * atan(1.0);
  double              PI2 = 2.0 * PI;

  for (int j = 0; j < m_Nc; ++j) {
    int    j2  = j * 2 * m_Nc;
    double rn1 = rand->get();
    double rn2 = rand->get();
    double rn3 = rand->get();
    double rn4 = rand->get();
    double rn5 = rand->get();

    double c1   = 1.0 - 2.0 * rn1;
    double s1   = sqrt(1.0 - c1 * c1);
    double v1_2 = s1 * cos(PI2 * rn2);
    double v1_3 = s1 * sin(PI2 * rn2);

    va[j2 + 0] = c1 * cos(PI2 * rn3);
    va[j2 + 1] = c1 * sin(PI2 * rn3);
    va[j2 + 2] = v1_2 * cos(PI2 * rn4);
    va[j2 + 3] = v1_2 * sin(PI2 * rn4);
    va[j2 + 4] = v1_3 * cos(PI2 * rn5);
    va[j2 + 5] = v1_3 * sin(PI2 * rn5);
  }

  reunit();

  return *this;
}


//==========================================================
//==================================================END=====