allow to define Tmax from the cmd line.

src/external/LF_GL/PGKT_LF.cpp

@@ -38,10 +38,9 @@
  * <p>CTOR
  *
  * @param param param[0]=field (G), param[1]=width (1/us), param[2]=hopp (1/us)
- * @param tt time vector in (ns)
+ * @param tmax maxium of time to be used in (us).
- * @param pol polarization vector
  */
-PGKT_LF::PGKT_LF(std::vector<double> &param) : fParam(param)
+PGKT_LF::PGKT_LF(std::vector<double> &param, const double tmax) : fParam(param), fTmax(tmax)
 {
   fParam[0] *= fGammaMu;
 
@@ -62,11 +61,10 @@ int PGKT_LF::DynamicGaussKTLF()
   // check if there is an empty time vector, if yes create one, otherwise use the given one
   if (fTime.size() == 0) {
     double t = 0.0;
-    unsigned int i=0;
     do {
-      t = i++ * 1.0e-4; // 100ps steps
+      t += 1.0e-4; // 100ps steps
       fTime.push_back(t);
-    } while (t < 10.0);
+    } while (t < fTmax);
   }
   fPol.resize(fTime.size());
 
@@ -101,15 +99,17 @@ int PGKT_LF::DynamicGaussKTLF()
  */
 void PGKT_LF::CalculateDynKTLF()
 {
-  const double Tmax = 20.0; // 20 us
+  double tmax = 20.0; // 20 us
-  unsigned int N = static_cast<unsigned int>(16.0*Tmax*fParam[0]);
+  if (tmax < fTmax)
+    tmax = fTmax;
+  unsigned int N = static_cast<unsigned int>(16.0*tmax*fParam[0]);
 
   // check if width is very high
   if (fParam[1] > 0.1) {
     double tmin = 20.0;
     tmin = fabs(sqrt(3.0)/fParam[1]);
 
-    unsigned int Nrate = static_cast<unsigned int>(25.0 * Tmax / tmin);
+    unsigned int Nrate = static_cast<unsigned int>(25.0 * tmax / tmin);
     if (Nrate > N) {
       N = Nrate;
     }
@@ -128,7 +128,7 @@ void PGKT_LF::CalculateDynKTLF()
   // calculate the P^(0)(t) exp(-nu t) vector
   std::vector<double> p0exp(N);
   double t = 0.0;
-  double dt = Tmax/N;
+  double dt = tmax/N;
   double nu = fParam[0] / (2.0*std::numbers::pi_v<double>);
   fDynLFdt = dt; // keep it since it is needed in GetDynKTLFValue()
   for (unsigned int i=0; i<N; i++) {

src/external/LF_GL/PGKT_LF.h

@@ -35,7 +35,7 @@
 class PGKT_LF
 {
   public:
-    PGKT_LF(std::vector<double> &param);
+    PGKT_LF(std::vector<double> &param, const double tmax);
 
     int DynamicGaussKTLF();
     bool IsValid() { return fValid; }
@@ -51,6 +51,7 @@ class PGKT_LF
     double fDynLFdt{0.0001};
 
     bool fValid{false};
+    double fTmax{10.0};
     std::vector<double> fParam;
     std::vector<double> fTime;
     std::vector<double> fPol;

src/external/LF_GL/PLGKT_LF.cpp

@@ -36,7 +36,7 @@
 #include "PGKT_LF.h"
 
 //-----------------------------------------------------------------------------
-PLGKT_LF::PLGKT_LF(std::vector<double> &param) : fParam(param)
+PLGKT_LF::PLGKT_LF(std::vector<double> &param, const double tmax) : fParam(param), fTmax(tmax)
 {
   if (DynamicLGKTLF() == 0)
     fValid = true;
@@ -70,8 +70,10 @@ PLGKT_LF::PLGKT_LF(std::vector<double> &param) : fParam(param)
  *                                  \right].
  * \f}
  *
- * <p>where f(r) = P_z(t, \Delta_{\rm G} = r\cdot\Delta_{\rm L}, B, \nu)
+ * <p>where
- *
+ * \f[
+ *   f(r) = P_z(t, \Delta_{\rm G} = r\cdot\Delta_{\rm L}, B, \nu) \cdot \sqrt{\frac{2}{\pi}} \frac{1}{r^2} \, \exp(-\frac{1}{2 r^2})
+ * \f]
  * @return 0 on success, >0 otherwise
  */
 int PLGKT_LF::DynamicLGKTLF()
@@ -84,9 +86,26 @@ int PLGKT_LF::DynamicLGKTLF()
   double scale, rr2;
   fTime.clear();
   fPol.clear();
+
+  if ((fParam[0] < 0.5) && (fParam[2] == 0.0)) { // very close to ZF -> use ZF
+    // generate time vector
+    double t = 0.0;
+    do {
+      t += 1.0e-4; // 100ps steps
+      fTime.push_back(t);
+    } while (t < fTmax);
+
+    // generate polarization vector
+    for (auto tt : fTime) {
+      fPol.push_back(0.333333333333 + 0.666666666667 * (1.0 - fParam[1]*tt) * exp(-fParam[1]*tt));
+    }
+
+    return 0;
+  }
+
   for (unsigned int i=0; i<rr.size(); i++) {
     pp[1] = rr[i] * fParam[1];
-    PGKT_LF gkt_lf(pp);
+    PGKT_LF gkt_lf(pp, fTmax);
     if (!gkt_lf.IsValid()) {
       std::cout << "**ERROR** in Gaussian LF calculation" << std::endl;
       return 2;

src/external/LF_GL/PLGKT_LF.h

@@ -34,7 +34,7 @@
 
 class PLGKT_LF {
   public:
-    PLGKT_LF(std::vector<double> &param);
+    PLGKT_LF(std::vector<double> &param, const double tmax);
 
     int DynamicLGKTLF();
     bool IsValid() { return fValid; }
@@ -44,6 +44,7 @@ class PLGKT_LF {
 
   private:
     bool fValid{false};
+    double fTmax{10.0};
     std::vector<double> fParam;
     std::vector<double> fTime;
     std::vector<double> fPol;

src/external/LF_GL/main.cpp

@@ -13,10 +13,11 @@
 void lf_gl_syntax()
 {
     std::cout << std::endl;
-    std::cout << "usage lg_gl [-p field width hopp -g [1,0] -o flnOut] | [-h]" << std::endl;
+    std::cout << "usage lg_gl [-p field width hopp [-t tmax] -g [1|0] -o flnOut] | [-h]" << std::endl;
     std::cout << " -p: field in (G)" << std::endl;
     std::cout << "     width in (1/us)" << std::endl;
     std::cout << "     hopp  in (1/us)" << std::endl;
+    std::cout << " -t: tmax in (us). Default: tmax=10 (us)." << std::endl;
     std::cout << " -g: 1=Gaussian LF data, 0=Gaussian averaged -> Lorentz LF data" << std::endl;
     std::cout << " -o: flnOut = output file name." << std::endl;
     std::cout << " -h: this help" << std::endl;
@@ -47,7 +48,7 @@ int lf_gl_write(const std::string fln, const bool gaussian_only, const std::vect
 int main(int argc, char *argv[])
 {
   std::vector<double> param{0.0, 0.0, 0.0};
-  double field{0.0};
+  double tmax{10.0};
   bool gaussian_only = false;
   std::string flnOut="";
 
@@ -56,37 +57,91 @@ int main(int argc, char *argv[])
     return 0;
   }
   double dval;
+  std::size_t pos{0};
   for (int i=1; i<argc; i++) {
     if (!strcmp(argv[i], "-p")) {
+      if (i+3 >= argc) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling parameters. Not enough input present." << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
       try {
-        dval = std::stod(argv[i+1]);
+        dval = std::stod(argv[i+1], &pos);
       } catch (...) {
         std::cout << std::endl;
-        std::cout << "**ERROR** in handling of field." << std::endl;
+        std::cout << "**ERROR** in handling of field: '" << argv[i+1] << "'" << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
+      if (pos != strlen(argv[i+1])) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling of field: '" << argv[i+1] << "'" << std::endl;
         lf_gl_syntax();
         return 1;
       }
       param[0] = dval;
       try {
-        dval = std::stod(argv[i+2]);
+        dval = std::stod(argv[i+2], &pos);
       } catch (...) {
         std::cout << std::endl;
         std::cout << "**ERROR** in handling of width." << std::endl;
         lf_gl_syntax();
         return 1;
       }
+      if (pos != strlen(argv[i+2])) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling of field: '" << argv[i+1] << "'" << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
       param[1] = dval;
       try {
-        dval = std::stod(argv[i+3]);
+        dval = std::stod(argv[i+3], &pos);
       } catch (...) {
         std::cout << std::endl;
         std::cout << "**ERROR** in handling of hopp." << std::endl;
         lf_gl_syntax();
         return 1;
       }
+      if (pos != strlen(argv[i+3])) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling of field: '" << argv[i+1] << "'" << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
       param[2] = dval;
       i += 3;
+    } else if (!strcmp(argv[i], "-t")) {
+      if (i+1 >= argc) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling option -g. Not enough input present." << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
+      try {
+        dval = std::stod(argv[i+1], &pos);
+      } catch (...) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling of tmax." << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
+      if (pos != strlen(argv[i+1])) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling of field: '" << argv[i+1] << "'" << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
+      tmax = dval;
+      i++;
     } else if (!strcmp(argv[i], "-g")) {
+      if (i+1 >= argc) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling option -g. Not enough input present." << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
       if (!strcmp(argv[i+1], "1")) {
         gaussian_only = true;
       } else if (!strcmp(argv[i+1], "0")) {
@@ -99,6 +154,12 @@ int main(int argc, char *argv[])
       }
       i++;
     } else if (!strcmp(argv[i], "-o")) {
+      if (i+1 >= argc) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling option -o. file name missing?" << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
       flnOut = std::string(argv[i+1]);
     }
   }
@@ -120,7 +181,7 @@ int main(int argc, char *argv[])
   std::vector<double> tt;
   std::vector<double> pol;
   if (gaussian_only) {
-    PGKT_LF gkt_lf(param);
+    PGKT_LF gkt_lf(param, tmax);
     if (!gkt_lf.IsValid()) {
       std::cout << "**ERROR** in Gaussian LF calculation" << std::endl;
       return 2;
@@ -128,7 +189,7 @@ int main(int argc, char *argv[])
     tt  = gkt_lf.GetTime();
     pol = gkt_lf.GetPol();
   } else {
-    PLGKT_LF lgkt_lf(param);
+    PLGKT_LF lgkt_lf(param, tmax);
     if (!lgkt_lf.IsValid()) {
       std::cout << "**ERROR** in Gaussian/Lorentzian LF calculation" << std::endl;
       return 2;