#include <stdlib.h>
#include <qdp.h>

extern int nt;

int
timeslice_func(int x[])
{
  return x[3];
}

void
w_meson2(QDP_DiracPropagator *src1, QDP_DiracPropagator *src2)
{
  QLA_Complex *prop0;
  QDP_DiracPropagator *localmat;
  QDP_Subset *timeslices;

  prop0 = (QLA_Complex *) malloc( nt * sizeof(QLA_Complex) );

  localmat = QDP_create_P();

  timeslices = QDP_create_subset(timeslice_func, nt);

  /* Multiply by gamma_0 for gamma_5 - gamma_5*gamma_0 correlator */
  QDP_P_eq_P_times_gamma( localmat, src1, 0, QDP_all );

  /* trace over propagators */
  QDP_c_eq_P_dot_P_multi( prop0, localmat, src2, timeslices, nt );

  QDP_destroy_subset(timeslices);
  QDP_destroy_P(localmat);
  free(prop0);
}