added ZF, LF geometry, dynamic approximation width/hopp << 1.

src/external/LF_GL/main.cpp

@@ -13,10 +13,12 @@
 void lf_gl_syntax()
 {
     std::cout << std::endl;
-    std::cout << "usage lg_gl [-p field width hopp [-t tmax] -g [1|0] -o flnOut] | [-h]" << std::endl;
+    std::cout << "usage lg_gl [[-p field width hopp | -a 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 << " -a: 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;
@@ -25,11 +27,45 @@ void lf_gl_syntax()
 }
 
 // ----------------------------------------------------------------------------
-int lf_gl_write(const std::string fln, const bool gaussian_only, const std::vector<double> &param,
+/**
+ * <p>ZF Approximation for dynamic (LF geometry).
+ * Gaussian approximation:   Abragam function, A. Yaouanc and P. Dalmas, p. 123, Eq. (6.62)
+ * Lorentzian approximation: A. Yaouanc and P. Dalmas, p. 130, Eq. (6.89)
+ *
+ * @param gaussian if true, Gaussian approximation, otherwise Lorentzian
+ * @param param [0]: field = 0, [1]: width (1/us), [2]: hopp (1/us)
+ * @param tmax calculate upto tmax
+ * @param tt time vector
+ * @param pol polarization vector
+ */
+void lf_zf_approx(const bool gaussian, const std::vector<double> &param,
+                  const double tmax, std::vector<double> &tt, std::vector<double> &pol)
+{
+  tt.clear();
+  pol.clear();
+  double t=0, dval;
+  do {
+    tt.push_back(t);
+    t += 0.001;
+    if (gaussian) {
+      dval = exp(-2.0*pow(param[1]/param[2], 2.0)*(exp(-param[2]*t)-1.0+param[2]*t));
+    } else {
+      dval = exp(-sqrt(4.0*pow(param[1]/param[2], 2.0)*(exp(-param[2]*t)-1.0+param[2]*t)));
+    }
+    pol.push_back(dval);
+  } while (t <= tmax);
+}
+
+// ----------------------------------------------------------------------------
+int lf_gl_write(const std::string fln, const bool zf_approx,
+                const bool gaussian_only,
+                const std::vector<double> &param,
                 const std::vector<double> &tt, const std::vector<double> &pol)
 {
   std::ofstream fout(fln, std::ofstream::out);
 
+  if (zf_approx)
+    fout << "# ZF Approximation." << std::endl;
   if (gaussian_only)
     fout << "# Gaussian only" << std::endl;
   else
@@ -50,9 +86,10 @@ int main(int argc, char *argv[])
   std::vector<double> param{0.0, 0.0, 0.0};
   double tmax{10.0};
   bool gaussian_only = false;
+  bool zf_approx = false;
   std::string flnOut="";
 
-  if (argc < 9) {
+  if (argc < 8) {
     lf_gl_syntax();
     return 0;
   }
@@ -62,7 +99,7 @@ int main(int argc, char *argv[])
     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;
+        std::cout << "**ERROR** in handling parameters '-p'. Not enough input present." << std::endl;
         lf_gl_syntax();
         return 1;
       }
@@ -112,6 +149,46 @@ int main(int argc, char *argv[])
       }
       param[2] = dval;
       i += 3;
+    } else if (!strcmp(argv[i], "-a")) {
+      if (i+2 >= argc) {
+        std::cout << std::endl;
+        std::cout << "**ERROR** in handling parameters for '-a'. Not enough input present." << std::endl;
+        lf_gl_syntax();
+        return 1;
+      }
+      param[0] = 0; // ZF approximation
+      try {
+          dval = std::stod(argv[i+1], &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+1])) {
+          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+2], &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+2])) {
+          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 += 2;
+      zf_approx = true;
     } else if (!strcmp(argv[i], "-t")) {
       if (i+1 >= argc) {
         std::cout << std::endl;
@@ -176,10 +253,14 @@ int main(int argc, char *argv[])
     std::cout << "Gaussian LF only" << std::endl;
   else
     std::cout << "Gaussian averaged -> Lorentz LF" << std::endl;
+  std::cout << "ZF Approx: " << zf_approx << std::endl;
   std::cout << "flnOut: " << flnOut << std::endl;
 
   std::vector<double> tt;
   std::vector<double> pol;
+  if (zf_approx) {
+    lf_zf_approx(gaussian_only, param, tmax, tt, pol);
+  } else { // full calculation
     if (gaussian_only) {
       PGKT_LF gkt_lf(param, tmax);
       if (!gkt_lf.IsValid()) {
@@ -197,8 +278,9 @@ int main(int argc, char *argv[])
       tt  = lgkt_lf.GetTime();
       pol = lgkt_lf.GetPol();
     }
+  }
 
-  lf_gl_write(flnOut, gaussian_only, param, tt, pol);
+  lf_gl_write(flnOut, zf_approx, gaussian_only, param, tt, pol);
 
   return 0;
 }