/***************************************************************************

  PLGKT_LF.cpp

  Author: Andreas Suter
  e-mail: andreas.suter@psi.ch

***************************************************************************/

/***************************************************************************
 *   Copyright (C) 2007-2025 by Andreas Suter                              *
 *   andreas.suter@psi.ch                                                  *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <iostream>
#include <cmath>
#include <algorithm>
#include <numbers>

#include "PLGKT_LF.h"
#include "PGKT_LF.h"

//-----------------------------------------------------------------------------
PLGKT_LF::PLGKT_LF(std::vector<double> &param) : fParam(param)
{
  if (DynamicLGKTLF() == 0)
    fValid = true;
}

//-----------------------------------------------------------------------------
/**
 * <p>The weight density \f$\rho_{\Delta_{\rm L}}(\Delta_{\rm G}) dG_{\rm G} =
 * \sqrt{\frac{2}{\pi}} \frac{1}{r^2} \exp\left(-\frac{1}{2 r^2}\right) dr\f$, with
 * \f$r = \frac{\Delta_{\rm G}}{\Delta_{\rm L}}\f$
 * (see A. Yaouanc and P. Dalmas de Reotier, ``Muon Spin Rotation, Relaxation, and Resonance'', p.129)
 *
 * <p>Integration method: Simpson on the \f$r-\f$ intervals:
 * \f$[0.2, 0.6], [0.6, 1.0], [1.0, 3.0], [3.0, 5.0], [5.0, 7.5], [7.5, 10.0], [10.0, 55.0], [55.0, 100.0]\f$
 *
 * <p>This leads to
 * \f{eqnarray*}{
 *   P_z(t, \Delta_{\rm L}, B, \nu) &=& \frac{0.4}{6} \left[f(0.2) + 4 f(0.4) + f(0.6) \right] +
 *                                      \frac{0.4}{6} \left[f(0.6) + 4 f(0.8) + f(1.0) \right] \\
 *                                  & & \frac{2}{6} \left[f(1.0) + 4 f(2.0) + f(3.0) \right] +
 *                                      \frac{2}{6} \left[f(3.0) + 4 f(4.0) + f(5.0) \right] \\
 *                                  & & \frac{2.5}{6} \left[f(5.0) + 4 f(6.25) + f(7.5) \right] +
 *                                      \frac{2.5}{6} \left[f(7.5) + 4 f(8.75) + f(10.0) \right] \\
 *                                  & & \frac{45}{6} \left[f(10.0) + 4 f(32.5) + f(55.0) \right] +
 *                                      \frac{45}{6} \left[f(55.0) + 4 f(77.5) + f(100.0) \right] \\
 *                                  &=& \frac{1}{6} \left[
 *                                      0.4 f(0.2) + 1.6 f(0.4) + 0.8 f(0.6) + 1.6 f(0.8) + 2.4 f(1.0) +\\
 *                                  & & 8.0 f(2.0) + 4 f(3.0) + 8 f(4.0) + 4.5 f(5.0) + 10.0 f(6.25) +\\
 *                                  & & 5 f(7.5) + 10.0 f(8.75) + 47.5 f(10.0) + 180.0 f(32.5) + 90 f(55.0) +\\
 *                                  & & 180 f(77.5) + 45 f(100.0)
 *                                  \right].
 * \f}
 *
 * <p>where f(r) = P_z(t, \Delta_{\rm G} = r\cdot\Delta_{\rm L}, B, \nu)
 *
 * @return 0 on success, >0 otherwise
 */
int PLGKT_LF::DynamicLGKTLF()
{
  std::vector<double> rr={0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 6.25, 7.5, 8.75, 10.0, 32.5, 55.0, 77.5, 100.0};
  std::vector<double> ww={0.4, 1.6, 0.8, 1.6, 2.4, 8.0, 4.0, 8.0, 4.5, 10.0, 5.0, 10.0, 47.5, 180.0, 90.0, 180.0, 45.0};

  std::vector<double> pp={fParam[0], 0.0, fParam[2]};
  std::vector<double> pol;
  double scale, rr2;
  fTime.clear();
  fPol.clear();
  for (unsigned int i=0; i<rr.size(); i++) {
    pp[1] = rr[i] * fParam[1];
    PGKT_LF gkt_lf(pp);
    if (!gkt_lf.IsValid()) {
      std::cout << "**ERROR** in Gaussian LF calculation" << std::endl;
      return 2;
    }
    if (fTime.size()==0) {
      fTime = gkt_lf.GetTime();
      fPol.resize(fTime.size());
    }
    rr2 = pow(rr[i],2.0);
    scale = ww[i]*exp(-0.5/rr2)/rr2;
    pol.clear();
    pol = gkt_lf.GetPol();
    for (unsigned int j=0; j<fPol.size(); j++)
      fPol[j] += scale*pol[j];
  }

  scale = pow(2.0/std::numbers::pi_v<double>, 0.5)/6.0;
  std::transform(fPol.begin(), fPol.end(), fPol.begin(), [&scale](double el) { return el *= scale;});

  return 0;
}