Working multithreaded version backup

bachelorarbeit/Cargo.lock

@@ -456,6 +456,7 @@ name = "bachelorarbeit"
 version = "0.1.0"
 dependencies = [
  "ahash",
+ "dashmap",
  "eframe",
  "egui",
  "egui_extras",
@@ -465,6 +466,9 @@ dependencies = [
  "indicatif",
  "nohash-hasher",
  "num",
+ "ray",
+ "rayon",
+ "scc",
  "wasm-bindgen-futures",
 ]
 
@@ -849,6 +853,25 @@ dependencies = [
  "cfg-if",
 ]
 
+[[package]]
+name = "crossbeam-deque"
+version = "0.8.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "613f8cc01fe9cf1a3eb3d7f488fd2fa8388403e97039e2f73692932e291a770d"
+dependencies = [
+ "crossbeam-epoch",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.9.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
+dependencies = [
+ "crossbeam-utils",
+]
+
 [[package]]
 name = "crossbeam-utils"
 version = "0.8.19"
@@ -871,6 +894,20 @@ version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "96a6ac251f4a2aca6b3f91340350eab87ae57c3f127ffeb585e92bd336717991"
 
+[[package]]
+name = "dashmap"
+version = "6.0.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "804c8821570c3f8b70230c2ba75ffa5c0f9a4189b9a432b6656c536712acae28"
+dependencies = [
+ "cfg-if",
+ "crossbeam-utils",
+ "hashbrown",
+ "lock_api",
+ "once_cell",
+ "parking_lot_core",
+]
+
 [[package]]
 name = "derivative"
 version = "2.2.0"
@@ -1053,6 +1090,12 @@ dependencies = [
  "egui",
 ]
 
+[[package]]
+name = "either"
+version = "1.13.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0"
+
 [[package]]
 name = "emath"
 version = "0.27.2"
@@ -2560,6 +2603,32 @@ version = "0.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
 
+[[package]]
+name = "ray"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4706e01695c83e94a0a61af525170acefa9061efc4567dc15ef183bea3b328a8"
+
+[[package]]
+name = "rayon"
+version = "1.10.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa"
+dependencies = [
+ "either",
+ "rayon-core",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.12.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2"
+dependencies = [
+ "crossbeam-deque",
+ "crossbeam-utils",
+]
+
 [[package]]
 name = "redox_syscall"
 version = "0.3.5"
@@ -2679,6 +2748,15 @@ dependencies = [
  "winapi-util",
 ]
 
+[[package]]
+name = "scc"
+version = "2.1.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "05ccfb12511cdb770157ace92d7dda771e498445b78f9886e8cdbc5140a4eced"
+dependencies = [
+ "sdd",
+]
+
 [[package]]
 name = "scoped-tls"
 version = "1.0.1"
@@ -2704,6 +2782,12 @@ dependencies = [
  "tiny-skia",
 ]
 
+[[package]]
+name = "sdd"
+version = "2.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "177258b64c0faaa9ffd3c65cd3262c2bc7e2588dbbd9c1641d0346145c1bbda8"
+
 [[package]]
 name = "serde"
 version = "1.0.199"

bachelorarbeit/Cargo.toml

@@ -15,6 +15,10 @@ egui_plot = "0.27.2"
 indicatif = "0.17.8"
 nohash-hasher = "0.2.0"
 ahash = "0.8.11"
+ray = "0.2.0"
+rayon = "1.10.0"
+scc = "2.1.6"
+dashmap = "6.0.1"
 
 [target.'cfg(target_arch = "wasm32")'.dependencies]
 wasm-bindgen-futures = "0.4"

bachelorarbeit/plotter.py

@@ -0,0 +1,125 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scienceplots
+
+i = 1
+
+names = [
+    "omnes_integrand",
+    "omnes_integrand_zoom",
+    "omnes_integrand_tan",
+    "phi0_integrand",
+    "phi0_integrand_zoom",
+    "phi0_integrand_tan",
+    "delta",
+    "omnes",
+    "phi0",
+    "phi0_delta_rel_diff",
+]
+
+files = [open("exports/" + n + ".txt") for n in names]
+
+export_filenames = ["figures/" + n + ".png" for n in names]
+
+titles = [
+    r"$\frac{\delta_1^1(s') - \delta_1^1(s)}{s'(s'-s)}$, $\sqrt{s}=1$ GeV",
+    r"$\frac{\delta_1^1(s') - \delta_1^1(s)}{s'(s'-s)}$, $\sqrt{s}=1$ GeV",
+    r"$\frac{1}{\cos^2(u)}\frac{\delta_1^1(\tan{u}) - \delta_1^1(s)}{\tan{u}(\tan{u}-s)}$, $\sqrt{s}=1$ GeV",
+    r"$\frac{\ln{|\frac{F(s')}{F(s)}|^2}}{(x^2+s_0)(x^2+s_0-s)}$, $\sqrt{s}=1$ GeV",
+    r"$\frac{\ln{|\frac{F(s')}{F(s)}|^2}}{(x^2+s_0)(x^2+s_0-s)}$, $\sqrt{s}=1$ GeV",
+    r"$\frac{1}{\cos^2(u)}\frac{\ln{|\frac{F(\tan^2{u}+s_0)}{F(s)}|^2}}{(\tan^2{u}+s_0)(\tan^2{u}+s_0-s)}$, $\sqrt{s}=1$ GeV",
+    "",
+    "",
+    "",
+    "",
+]
+
+xlabels = [
+    r"$\sqrt{s}$ [GeV]",
+    r"$\sqrt{s}$ [GeV]",
+    r"$u$",
+    r"$x$",
+    r"$x$",
+    r"$u$",
+    r"$\sqrt{s}$ [GeV]",
+    r"$\sqrt{s}$ [GeV]",
+    r"$\sqrt{s}$ [GeV]",
+    r"$\sqrt{s}$ [GeV]",
+]
+
+ylabels = [
+    r"integrand",
+    r"integrand",
+    r"integrand",
+    r"integrand",
+    r"integrand",
+    r"integrand",
+    r"$\delta_1^1(s)$",
+    r"$|\Omega(s)|^2$",
+    r"$\phi_0(s)$",
+    r"$\frac{\phi_0(s)}{\delta_1^1(s)}$",
+]
+
+y_logarithmic = [
+    False,
+    False,
+    False,
+    False,
+    False,
+    False,
+    False,
+    True,
+    False,
+    False,
+]
+
+x_sqrt = [
+    True,
+    True,
+    False,
+    False,
+    False,
+    False,
+    True,
+    True,
+    True,
+    True,
+]
+
+for i in range(len(names)):
+    print(i)
+    plt.figure()
+    filecont = []
+    for line in files[i]:
+        filecont.append(line.split(','))
+
+    data = []
+    for curve in filecont:
+        curve_numbers = []
+        for val in curve:
+            try:
+                curve_numbers.append(float(val))
+            except ValueError:
+                continue
+        if x_sqrt[i]:
+            curve_numbers_paired = np.array([[x**0.5 for x in curve_numbers[0::2]], curve_numbers[1::2]])
+        else:
+            curve_numbers_paired = np.array([curve_numbers[0::2], curve_numbers[1::2]])
+        data.append(curve_numbers_paired)
+
+    # Plotting
+
+    plt.style.use(['science'])
+
+    plt.figure(figsize=(4,3))
+    plt.title(titles[i], fontsize = 15)
+    if y_logarithmic[i]:
+        plt.yscale('log')
+    for j in range(len(data)):
+        plt.plot(data[j][0], data[j][1])
+    plt.xlabel(xlabels[i], fontsize = 12)
+    plt.ylabel(ylabels[i], fontsize = 12)
+    plt.grid(which='major', color='#000000', linestyle='-', alpha = 0.3)
+    plt.grid(which='minor', color='#000000', linestyle='-', alpha=0.1)
+    plt.savefig(export_filenames[i], dpi=400)
+    plt.close()

bachelorarbeit/simpleplot.py

@@ -0,0 +1,133 @@
+import matplotlib.pyplot as plt
+import numpy as np
+data = [
+    [
+        0.00046921958193724134,
+        7.358616407326934e-5,
+        8.723023031975696e-5,
+        8.153727384807574e-5,
+        9.090152038249144e-5,
+        8.414070625384351e-5,
+        8.612395652507487e-5,
+        8.776483357408882e-5,
+        8.703793625031153e-5,
+        8.595498892316478e-5,
+    ],
+    [
+        0.00030936685630689187,
+        7.978835482247426e-5,
+        3.943681093321327e-5,
+        3.438668332045314e-5,
+        2.0189312051055452e-5,
+        2.874604966751626e-5,
+        2.700407617384215e-5,
+        2.4825341293444048e-5,
+        2.5628809087074877e-5,
+        2.678726585980673e-5,
+    ],
+    [
+        0.0003035764863147383,
+        0.00015586851063420504,
+        2.4490574906943507e-5,
+        1.7902688270066136e-5,
+        1.6650197575129866e-5,
+        7.59094323854459e-6,
+        3.342059771260786e-6,
+        2.8747982954158147e-6,
+        5.142103307598234e-7,
+        2.69394593033212e-6,
+    ],
+    [
+        0.00031160663209994777,
+        5.37083752694123e-5,
+        1.5304054572173875e-5,
+        9.256441461102938e-7,
+        6.723313460188507e-6,
+        8.254317026779034e-6,
+        1.4014951930096942e-6,
+        2.3551311784109075e-6,
+        5.980254742521396e-7,
+        3.597584800507647e-6,
+    ],
+    [
+        0.0003096809009284218,
+        7.174695932898878e-5,
+        1.8277500006158576e-5,
+        4.939461074715545e-6,
+        3.07613414929353e-6,
+        6.533050137313978e-6,
+        4.1030853825674285e-8,
+        2.2538858418519325e-6,
+        9.132272886791526e-7,
+        5.169395710824531e-6,
+    ],
+    [
+        0.0003078627597422612,
+        0.000106444698403374,
+        2.0882908283836876e-5,
+        1.0885016441308792e-5,
+        2.1767025349062052e-5,
+        3.2850613408408336e-6,
+        2.3900967810464024e-6,
+        5.198947351736649e-6,
+        2.426213355777307e-6,
+        8.623114753358863e-7,
+    ],
+    [
+        0.00030836014251955923,
+        9.61302294200106e-5,
+        2.0296397885921635e-5,
+        9.659277290574586e-6,
+        9.207541751266035e-6,
+        2.0492995951304493e-6,
+        2.5594252984761923e-6,
+        7.261965937477299e-7,
+        3.73917648510691e-7,
+        1.2425543701732877e-6,
+    ],
+    [
+        0.0003088364016595735,
+        8.721442142234359e-5,
+        1.962601701765987e-5,
+        8.169525936518873e-6,
+        3.1747640369950147e-6,
+        5.728087832945761e-7,
+        1.945445256779088e-6,
+        1.0661919567223066e-6,
+        1.1508338138010998e-6,
+        6.123502085397803e-7,
+    ],
+    [
+        0.0003086558304487008,
+        9.080158735530475e-5,
+        1.9872981781188237e-5,
+        8.736252950058976e-6,
+        6.22104723513317e-6,
+        3.8569119342746205e-7,
+        2.1402028302919263e-6,
+        9.468634343257065e-8,
+        6.590138900142151e-7,
+        2.6713974365932813e-9,
+    ],
+    [
+        0.0003083751385098976,
+        9.654880744447425e-5,
+        2.0235402455659468e-5,
+        9.544949382433998e-6,
+        1.7823923909920936e-5,
+        1.6227578926164554e-6,
+        2.2769853512683014e-6,
+        4.159153095706358e-6,
+        1.8888938431160796e-6,
+        5.116018929607336e-7,
+    ],
+]
+n = [100,500,1000,2000,5000,7500,10000,15000,25000,50000];
+
+
+for n_o in range(len(data)):
+    # plt.scatter([n_o] * len(data[n_o]), data[n_o], label=str(n[n_o]))
+    plt.plot(n, data[n_o], label=str(n[n_o]))
+    plt.plot
+plt.legend()
+plt.show()

bachelorarbeit/src/calc.rs

@@ -1,8 +1,11 @@
 use core::f64::consts::PI;
 use gauss_quad::GaussLegendre;
-use std::collections::HashMap;
 use std::io;
-use std::rc::Rc;
+use std::sync::Arc;
+use rayon::prelude::*;
+// use scc::DashMap;
+
+use dashmap::DashMap;
 
 use std::{fs::File, io::BufReader};
 
@@ -10,63 +13,68 @@ use crate::gq_storage;
 
 pub type Complex = num::complex::Complex<f64>;
 
-const M_ROH: f64 = 0.7736;
+pub const M_ROH: f64 = 0.7736;
-const M_PI: f64 = 0.13957;
+pub const M_PI: f64 = 0.13957;
 pub const S0: f64 = 4.0 * M_PI * M_PI;
-const M_K: f64 = 0.496;
+pub const M_K: f64 = 0.496;
-const B0: f64 = 1.055;
+pub const B0: f64 = 1.055;
-const B1: f64 = 0.15;
+pub const B1: f64 = 0.15;
-const LAMBDA1: f64 = 1.57;
+pub const LAMBDA1: f64 = 1.57;
-const LAMBDA2: f64 = -1.96;
+pub const LAMBDA2: f64 = -1.96;
-const LAMBDA_HIGH: f64 = 10.0;
+pub const LAMBDA_HIGH: f64 = 10.0;
 pub const EPSILON: f64 = 1e-15;
-const S_MID_CUTOFF: f64 = 1.42 * 1.42;
+pub const S_MID_CUTOFF: f64 = 1.42;
-const INF: f64 = 1e6;
+pub const INF: f64 = 1e6;
 
 #[derive(Default)]
 pub struct Cache {
-    pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>), nohash_hasher::BuildNoHashHasher<u32>>,
+    pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>), nohash_hasher::BuildNoHashHasher<u32>>,
-    pub delta_lut: HashMap<u64, f64, nohash_hasher::BuildNoHashHasher<u64>>,
+    pub delta_lut: DashMap<u64, f64, nohash_hasher::BuildNoHashHasher<u64>>,
-    pub omnes_lut: HashMap<
+    pub omnes_lut: DashMap<
         u32,
-        HashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>,
+        DashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>,
         nohash_hasher::BuildNoHashHasher<u32>,
     >, // First key is n, second is s
 }
 
 // pub struct Cache {
-//     pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>), RandomState>,
+//     pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>), RandomState>,
-//     pub delta_lut: HashMap<u64, f64, RandomState>,
+//     pub delta_lut: DashMap<u64, f64, RandomState>,
-//     pub omnes_lut: HashMap<
+//     pub omnes_lut: DashMap<
 //         u32,
-//         HashMap<u64, Complex, RandomState>,
+//         DashMap<u64, Complex, RandomState>,
 //         RandomState,
 //     >, // First key is n, second is s
 // }
 
 // pub struct Cache {
-//     pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>)>,
+//     pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>)>,
-//     pub delta_lut: HashMap<u64, f64>,
+//     pub delta_lut: DashMap<u64, f64>,
-//     pub omnes_lut: HashMap<u32, HashMap<u64, Complex>> // First key is n, second is s
+//     pub omnes_lut: DashMap<u32, DashMap<u64, Complex>> // First key is n, second is s
 // }
 
 impl Cache {
     pub fn from_file(filepath: &str) -> io::Result<Cache> {
+        // println!("Cache from_file");
         let file = File::open(filepath)?;
         let mut file = BufReader::new(file);
 
-        let me = Cache {
+        let gql = DashMap::with_hasher(Default::default());
-            gauss_quad_lut: gq_storage::deserialize(&mut file, |_| true)?
+        for ent in gq_storage::deserialize(&mut file, |_| true)?
                 .into_iter()
                 .map(|(roots, weights)| {
                     (
                         u32::try_from(roots.len()).unwrap(),
-                        (Rc::from(roots), Rc::from(weights)),
+                        (Arc::from(roots), Arc::from(weights)),
                     )
                 })
-                .collect(),
+        {
-            delta_lut: HashMap::with_hasher(Default::default()),
+            let _ = gql.insert(ent.0, ent.1);
-            omnes_lut: HashMap::with_hasher(Default::default()),
+        }
+        let me = Cache {
+            gauss_quad_lut: gql,
+            delta_lut: DashMap::with_hasher(Default::default()),
+            omnes_lut: DashMap::with_hasher(Default::default()),
         };
         Ok(me)
     }
@@ -74,18 +82,22 @@ impl Cache {
     #[allow(unused)]
     pub fn from_slice(data: &[u8]) -> io::Result<Cache> {
         let mut cursor = std::io::Cursor::new(data);
-        let me = Cache {
+        let mut gq = DashMap::with_hasher(Default::default());
-            gauss_quad_lut: gq_storage::deserialize(&mut cursor, |_| true)?
+        for ent in gq_storage::deserialize(&mut cursor, |_| true)?
                 .into_iter()
                 .map(|(roots, weights)| {
                     (
                         u32::try_from(roots.len()).unwrap(),
-                        (Rc::from(roots), Rc::from(weights)),
+                        (Arc::from(roots), Arc::from(weights)),
                     )
                 })
-                .collect(),
+        {
-            delta_lut: HashMap::with_hasher(Default::default()),
+            gq.insert(ent.0, ent.1);
-            omnes_lut: HashMap::with_hasher(Default::default()),
+        }
+        let me = Cache {
+            gauss_quad_lut: gq,
+            delta_lut: DashMap::with_hasher(Default::default()),
+            omnes_lut: DashMap::with_hasher(Default::default()),
         };
         Ok(me)
     }
@@ -93,12 +105,8 @@ impl Cache {
 
 // derivated values
 
-fn lambda_0() -> f64 {
+fn lambda_0(cache: &Cache) -> f64 {
-    delta((2.0 * M_K).powi(2))
+    delta_with_lut(cache, (2.0 * M_K).powi(2))
-}
-
-fn delta_mid_cutoff() -> f64 {
-    delta(S_MID_CUTOFF)
 }
 
 fn atan_shift(x: f64) -> f64 {
@@ -110,33 +118,31 @@ fn atan_shift(x: f64) -> f64 {
     }
 }
 
-fn integrate<G: FnMut(f64) -> f64>(
+fn integrate<G: Fn(f64) -> f64 + Send + Sync>(
     roots: &[f64],
     weights: &[f64],
     a: f64,
     b: f64,
-    mut f: G,
+    f: G,
     print_prog: bool,
 ) -> f64 {
-    let mut sum: f64 = 0.0;
-    let mut i: usize = 0;
     if roots.len() != weights.len() {
         panic!("Error: roots and weights are of different length");
     }
-    if print_prog {
+    let prog = indicatif::ProgressBar::new(roots.len() as u64);
-        println!("Calculating single integration: ");
+    let prog_ = prog.clone();
-        let prog = indicatif::ProgressBar::new(roots.len() as u64);
+    let sum:f64 = (weights.par_iter(), roots.par_iter()).into_par_iter().map(move |(weight,root)| {
-        while i < roots.len() {
+        let prog = &prog_;
-            sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
+        // println!("Orello");
-            i += 1;
+        let res = weight * f(root * (b - a) / 2.0 + (a + b) / 2.0);
+        if print_prog {
             prog.inc(1);
         }
+        res
+        // 0.0
+    }).sum();
+    if print_prog {
         prog.finish();
-    } else {
-        while i < roots.len() {
-            sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
-            i += 1;
-        }
     }
     sum * (b - a) / 2.0
 }
@@ -148,6 +154,7 @@ fn integrate_complex<G: FnMut(f64) -> Complex>(
     b: f64,
     mut f: G,
 ) -> Complex {
+    // println!("integrate_complex");
     let mut sum: Complex = Complex::new(0.0, 0.0);
     let mut i: usize = 0;
     if roots.len() != weights.len() {
@@ -170,77 +177,117 @@ fn k(s: f64) -> f64 {
     (s / 4.0 - M_PI.powi(2)).sqrt()
 }
 
-pub fn delta_with_lut(cache: &mut Cache, s: f64) -> f64 {
+pub fn delta_with_lut(cache: &Cache, s: f64) -> f64 {
     match cache.delta_lut.get(&s.to_bits()) {
         Some(val) => *val,
         None => {
-            let val = delta(s);
+            let val = delta(cache, s);
-            cache.delta_lut.insert(s.to_bits(), val);
+            let _ = cache.delta_lut.insert(s.to_bits(), val);
             val
         }
     }
 }
 
-pub fn delta(s: f64) -> f64 {
+
-    if s <= 2.0 * M_K {
+
+const BLEND_RANGE: f64 = 0.1;
+
+#[allow(unused)]
+pub fn delta_low(cache: &Cache, s: f64) -> f64 {
+    atan_shift(
+        (s.sqrt() / (2.0 * k(s).powi(3))
+            * (M_ROH.powi(2) - s)
+            * ((2.0 * M_PI.powi(3)) / (M_ROH.powi(2) * s.sqrt()) + B0 + B1 * omega(s)))
+        .powi(-1),
+    )
+}
+
+#[allow(unused)]
+pub fn delta_mid(cache: &Cache, s: f64) -> f64 {
+    let par = 0.5 * s.sqrt() / M_K - 1.0;
+    lambda_0(cache) + LAMBDA1 * par + LAMBDA2 * par.powi(2)
+}
+
+#[allow(unused)]
+pub fn delta_high(cache: &Cache, s: f64) -> f64 {
+    PI + (delta(cache, S_MID_CUTOFF.powi(2)) - PI) * (LAMBDA_HIGH * LAMBDA_HIGH + S_MID_CUTOFF)
+        / (LAMBDA_HIGH * LAMBDA_HIGH + s)
+}
+
+#[allow(unused)]
+pub fn delta_blend(cache: &Cache, s: f64) -> f64 {
+    if (s - 4.0*M_K.powi(2)).abs() < BLEND_RANGE {
+        let t = (s - 4.0*M_K.powi(2) + BLEND_RANGE) / (2.0*BLEND_RANGE);
+        (1.0 - t) * delta_low(cache, s) + t * delta_mid(cache, s)
+    } else if (s - S_MID_CUTOFF.powi(2)).abs() < BLEND_RANGE {
+        let t = (s - S_MID_CUTOFF.powi(2) + BLEND_RANGE) / (2.0*BLEND_RANGE);
+        (1.0 - t) * delta_mid(cache, s) + t * delta_high(cache, s)
+    } else {
+        delta(cache, s)
+    }
+}
+
+fn delta(cache: &Cache, s: f64) -> f64 {
+    let s_sqrt = s.sqrt();
+    if s <= 4.0 * M_K.powi(2) {
         atan_shift(
-            (s.sqrt() / (2.0 * k(s).powi(3))
+            (s_sqrt / (2.0 * k(s).powi(3))
                 * (M_ROH.powi(2) - s)
-                * ((2.0 * M_PI.powi(3)) / (M_ROH.powi(2) * s.sqrt()) + B0 + B1 * omega(s)))
+                * ((2.0 * M_PI.powi(3)) / (M_ROH.powi(2) * s_sqrt) + B0 + B1 * omega(s)))
             .powi(-1),
         )
-    } else if s <= S_MID_CUTOFF {
+    } else if s <= S_MID_CUTOFF.powi(2) {
-        let par = 0.5 * s.sqrt() / M_K - 1.0;
+        let par = 0.5 * s_sqrt / M_K - 1.0;
-        lambda_0() + LAMBDA1 * par + LAMBDA2 * par.powi(2)
+        lambda_0(cache) + LAMBDA1 * par + LAMBDA2 * par.powi(2)
     } else {
-        PI + (delta_mid_cutoff() - PI) * (LAMBDA_HIGH * LAMBDA_HIGH + S_MID_CUTOFF)
+        PI + (delta_with_lut(cache, S_MID_CUTOFF.powi(2)) - PI) * (LAMBDA_HIGH * LAMBDA_HIGH + S_MID_CUTOFF)
             / (LAMBDA_HIGH * LAMBDA_HIGH + s)
     }
 }
 
-pub fn omnes_integrand_tan(cache: &mut Cache, s_tick: f64, s: f64) -> Complex {
+pub fn omnes_integrand_tan(cache: &Cache, s_tick: f64, s: f64) -> Complex {
     let sub = s_tick.tan();
     (delta_with_lut(cache, sub) - delta_with_lut(cache, s))
         / (sub * (sub - s + Complex::new(0.0, EPSILON)) * s_tick.cos().powi(2))
 }
 
-pub fn omnes_integrand_notan(cache: &mut Cache, s_tick: f64, s: f64) -> Complex {
+pub fn omnes_integrand_notan(cache: &Cache, s_tick: f64, s: f64) -> Complex {
     (delta_with_lut(cache, s_tick) - delta_with_lut(cache, s))
         / (s_tick * (s_tick - s + Complex::new(0.0, EPSILON)))
 }
 
-pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
+pub fn omnes(cache: &Cache, s: f64, n: u32, use_tan: bool) -> Complex {
     if let Some(inner_lut) = cache.omnes_lut.get(&n) {
         if let Some(res) = inner_lut.get(&s.to_bits()) {
             return *res;
         }
     }
     if !cache.omnes_lut.contains_key(&n) {
-        cache
+        let _ = cache
             .omnes_lut
-            .insert(n, HashMap::with_hasher(Default::default()));
+            .insert(n, DashMap::with_hasher(Default::default()));
     }
 
-    let roots: Rc<[f64]>;
+    let roots: Arc<[f64]>;
-    let weights: Rc<[f64]>;
+    let weights: Arc<[f64]>;
-    match cache.gauss_quad_lut.get_mut(&n) {
+    match cache.gauss_quad_lut.get(&n) {
         Some(gq_values) => {
-            roots = Rc::clone(&gq_values.0);
+            roots = Arc::clone(&gq_values.0);
-            weights = Rc::clone(&gq_values.1);
+            weights = Arc::clone(&gq_values.1);
         }
         None => {
             let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
-            roots = Rc::from(gq_values.0.as_slice());
+            roots = Arc::from(gq_values.0.as_slice());
-            weights = Rc::from(gq_values.1.as_slice());
+            weights = Arc::from(gq_values.1.as_slice());
         }
     }
 
     if !cache.gauss_quad_lut.contains_key(&n) {
-        let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
-        cache
-            .gauss_quad_lut
-            .insert(n, (Rc::from(gq_values.0), Rc::from(gq_values.1)));
         println!("n_omnes={} is not included in the lookup file. Consider generating a new lookup file via the appropriate function in gq_storage.rs", n);
+        let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
+        let _ = cache
+            .gauss_quad_lut
+            .insert(n, (Arc::from(gq_values.0), Arc::from(gq_values.1)));
     }
 
     let intgrl = if use_tan {
@@ -254,9 +301,9 @@ pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
     };
     let omnes_res = (s / PI * (intgrl)).exp()
         * (S0 / Complex::new(S0 - s, -EPSILON)).powf(delta_with_lut(cache, s) / PI);
-    cache
+    let _ = cache
         .omnes_lut
-        .get_mut(&n)
+        .get(&n)
         .unwrap()
         .insert(s.to_bits(), omnes_res);
 
@@ -264,7 +311,7 @@ pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
 }
 
 pub fn phi0_integrand(
-    cache: &mut Cache,
+    cache: &Cache,
     s_tick: f64,
     s: f64,
     n_omnes: u32,
@@ -283,14 +330,8 @@ pub fn phi0_integrand(
         (false, true) => s_tick.powi(2) + S0,
         (false, false) => s_tick,
     };
-    // let omnes_s = if use_cutoff && s >= cutoff_s {
-    //     cutoff_factor * s.powf(cutoff_power)
-    // } else {
-    //     omnes(cache, s, n_omnes, use_tan_omnes)
-    // };
     let omnes_sub = if use_cutoff && sub >= cutoff_s {
         cutoff_factor * sub.powf(cutoff_power)
-        // Complex::new((0.00002f64).powf(cutoff_power), 0.0)
     } else {
         omnes(cache, sub, n_omnes, use_tan_omnes)
     };
@@ -306,7 +347,7 @@ pub fn phi0_integrand(
 }
 
 pub fn phi0(
-    cache: &mut Cache,
+    cache: &Cache,
     s: f64,
     n_omnes: u32,
     n_phi0: u32,
@@ -321,21 +362,23 @@ pub fn phi0(
 ) -> f64 {
     if !cache.gauss_quad_lut.contains_key(&n_phi0) {
         let gq_values = GaussLegendre::nodes_and_weights(n_phi0.try_into().unwrap());
-        cache
+        let _ = cache
             .gauss_quad_lut
-            .insert(n_phi0, (Rc::from(gq_values.0), Rc::from(gq_values.1)));
+            .insert(n_phi0, (Arc::from(gq_values.0), Arc::from(gq_values.1)));
         println!("n_omnes={} is not included in the lookup file. Consider generating a new lookup file via the appropriate function in gq_storage.rs", n_phi0);
     }
-    let (ref roots, ref weights) = cache.gauss_quad_lut.get(&n_phi0).unwrap();
+    let (roots, weights) = {
-    let (roots, weights) = (Rc::clone(roots), Rc::clone(weights));
+        let ent = cache.gauss_quad_lut.get(&n_phi0).unwrap();
+        (Arc::clone(&ent.0), Arc::clone(&ent.1))
+    };
 
     let analyt = if use_reduced_integrand {
-        omnes(cache, s, n_omnes, use_tan_omnes).norm().ln() / S0.sqrt()
-    } else {
         (omnes(cache, s, n_omnes, use_tan_omnes) / omnes(cache, S0, n_omnes, use_tan_omnes))
             .norm()
             .ln()
             / (s * (s - S0).powf(1.5))
+    } else {
+        omnes(cache, s, n_omnes, use_tan_omnes).norm().ln() / S0.sqrt()
     };
 
     let cutoff_factor = omnes(cache, cutoff_s, n_omnes, use_tan_omnes) * cutoff_s.powf(-cutoff_power);

bachelorarbeit/src/main.rs

@@ -3,17 +3,21 @@
 // - Reduziertes Integral aus Paper zu numerisch berechenbarer Form umformen
 
 use indicatif;
-use std::{collections::HashMap, time::Instant};
+use std::time::Instant;
+use std::fs;
+use core::f64::consts::PI;
 
 use egui_plot::{log_grid_spacer, AxisHints, GridInput, GridMark, Legend, Line, Plot};
 use num::{complex::ComplexFloat, Float};
+use eframe::Theme;
+// use scc::HashMap;
 
 mod calc;
 mod gq_storage;
 
 struct App {
     plot_data: Vec<PlotCurve>,
-    plots_available: [String; 7],
+    plots_available: [String; 9],
     calc_cache: calc::Cache,
     scaling_type: ScalingType,
     a: f64,
@@ -28,6 +32,7 @@ struct App {
     use_cutoff: bool,
     cutoff_power: f64,
     cutoff_s: f64,
+    export_filename: String,
 }
 
 #[derive(Clone)]
@@ -46,27 +51,31 @@ impl Default for App {
                 "Phi0".to_string(),
                 "Phi0 / delta".to_string(),
                 "s^cutoff_power".to_string(),
+                "Omnes integrand".to_string(),
                 "phi0 integrand".to_string(),
                 "Phi0_cut / Phi0".to_string(),
+                "Phi0 / delta".to_string(),
             ],
             calc_cache: calc::Cache::default(),
             scaling_type: ScalingType::default(),
             a: calc::S0,
             b: 9.0,
-            n_points: 500,
+            n_points:2500,
-            n_omnes: 500,
+            n_omnes:2500,
-            n_phi0: 500,
+            n_phi0:2500,
             
             // Substitutions params
             use_tan_omnes: true,
             use_tan_phi0: true,
             use_xsub_phi0: true,
-            use_reduced_integrand: true,
+            use_reduced_integrand: false,
 
             // Cutoff params
             use_cutoff: false,
             cutoff_power: -1.0,
             cutoff_s: 100.0, 
+
+            export_filename: "plot_data".to_string(),
         }
     }
 }
@@ -81,6 +90,7 @@ enum ScalingType {
 impl eframe::App for App {
     fn update(&mut self, ctx: &egui::Context, _frame: &mut eframe::Frame) {
         egui::CentralPanel::default().show(ctx, |ui| {
+            ui.set_min_width(ui.available_width());
             let cur_scaling_type = self.scaling_type;
 
             ui.vertical(|ui| {
@@ -95,12 +105,20 @@ impl eframe::App for App {
                                     "the lower bound of the interval that functions are calculated in",
                                 );
                                 float_text_edit_singleline(ui, &mut self.a, true);
+                                let seta_pi2 = ui.button("S0");
+                                if seta_pi2.clicked() {
+                                    self.a = calc::S0;
+                                }
                                 ui.end_row();
 
                                 ui.label("calc interval end:").on_hover_text(
                                     "the upper bound of the interval that functions are calculated in",
                                 );
                                 float_text_edit_singleline(ui, &mut self.b, true);
+                                let setb_pi2 = ui.button("pi/2");
+                                if setb_pi2.clicked() {
+                                    self.b = 1.5707963268;
+                                }
                                 ui.end_row();
 
                                 ui.label("N for omnes:").on_hover_text(
@@ -129,7 +147,7 @@ impl eframe::App for App {
                         let omnes_tan_cb =
                             ui.checkbox(&mut self.use_tan_omnes, "Use tan()-subst. for Omnes ('ot')");
                         if omnes_tan_cb.clicked() {
-                            self.calc_cache.omnes_lut = HashMap::with_hasher(Default::default()); // TODO split omnes_lut so that this is no longer necessary
+                            self.calc_cache.omnes_lut = dashmap::DashMap::with_hasher(Default::default()); // TODO split omnes_lut so that this is no longer necessary
                         }
                         ui.checkbox(&mut self.use_tan_phi0, "Use tan()-subst. for phi0 ('pt')");
                         ui.checkbox(&mut self.use_xsub_phi0, "Use second subst. for phi0 ('ss')");
@@ -178,6 +196,25 @@ impl eframe::App for App {
                             float_text_edit_singleline(ui, &mut self.cutoff_s, true);
                             ui.end_row();
                         });
+                        let textedit = egui::TextEdit::singleline(&mut self.export_filename);
+                        ui.add(textedit);
+                        let button_export = ui.button("Export data");
+                        if button_export.clicked() {
+                            let mut plot_data_str = String::new();
+                            for curve in &mut self.plot_data {
+                                for point in &mut curve.points {
+                                    plot_data_str.push_str(&point[0].to_string());
+                                    plot_data_str.push_str(",");
+                                    plot_data_str.push_str(&point[1].to_string());
+                                    plot_data_str.push_str(",");
+                                }
+                                plot_data_str.push_str("\n");
+                            }
+                            let mut filename = "./exports/".to_owned();
+                            filename.push_str(&self.export_filename);
+                            filename.push_str(".txt");
+                            fs::write(filename, plot_data_str).expect("Unable to write file");
+                        }
                     })
                 });
 
@@ -197,76 +234,52 @@ impl eframe::App for App {
                                 .collect();
                             match i {
                                 0 => {
-                                    let y_values_delta: Vec<f64> = if self.scaling_type
+                                    let y_values: Vec<f64> = x_values.iter().map(|&x| {
-                                        == ScalingType::Linear
+                                        let val = calc::delta_with_lut(&mut self.calc_cache, x);
-                                    {
+                                        match self.scaling_type {
-                                        x_values
+                                            ScalingType::Linear => {
-                                            .iter()
+                                                val
-                                            .map(|&x| calc::delta_with_lut(&mut self.calc_cache, x))
+                                            }
-                                            .collect()
+                                            ScalingType::Logarithmic => {
-                                    } else {
+                                                val.log10()
-                                        x_values
+                                            }
-                                            .iter()
+                                        }
-                                            .map(|&x| {
+                                    }).collect();
-                                                calc::delta_with_lut(&mut self.calc_cache, x)
-                                                    .log10()
-                                            })
-                                            .collect()
-                                    };
 
                                     self.plot_data.push(PlotCurve {
                                         points: x_values
                                             .iter()
-                                            .zip(y_values_delta.iter())
+                                            .zip(y_values.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name: format!("(#{}) Delta", self.plot_data.len() + 1),
-                                    })
+                                    });
                                 }
                                 1 => {
-                                    let y_values_omnes: Vec<f64> = if self.scaling_type
+                                    let cutoff_factor = calc::omnes(&mut self.calc_cache, self.cutoff_s, self.n_omnes, self.use_tan_omnes) * self.cutoff_s;
-                                        == ScalingType::Linear
+
-                                    {
+                                    let y_values: Vec<f64> = x_values.iter().map(|&x| {
-                                        print!("Calculating all plot points: ");
+                                       let val = calc::omnes(
-                                        let bar =
+                                            &mut self.calc_cache,
-                                            indicatif::ProgressBar::new(u64::from(self.n_points));
+                                            x,
-                                        let res = x_values
+                                            self.n_omnes,
-                                            .iter()
+                                            self.use_tan_omnes,
-                                            .map(|&x| {
+                                        ).abs().powi(2);
-                                                bar.inc(1);
+                                        match self.scaling_type {
-                                                calc::omnes(
+                                            ScalingType::Linear => {
-                                                    &mut self.calc_cache,
+                                                val
-                                                    x,
+                                            }
-                                                    self.n_omnes,
+                                            ScalingType::Logarithmic => {
-                                                    self.use_tan_omnes,
+                                                val.log10()
-                                                )
+                                            }
-                                                .abs()
+                                        }
-                                                .powi(2)
+                                    }).collect();
-                                            })
+
-                                            .collect();
-                                        bar.finish();
-                                        res
-                                    } else {
-                                        x_values
-                                            .iter()
-                                            .map(|&x| {
-                                                calc::omnes(
-                                                    &mut self.calc_cache,
-                                                    x,
-                                                    self.n_omnes,
-                                                    self.use_tan_omnes,
-                                                )
-                                                .abs()
-                                                .powi(2)
-                                                .log10()
-                                            })
-                                            .collect()
-                                    };
                                     self.plot_data.push(PlotCurve {
                                         points: x_values
                                             .iter()
-                                            .zip(y_values_omnes.iter())
+                                            .zip(y_values.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name: format!(
                                             "(#{}) Omnes, {} ot, n_o={}",
@@ -316,7 +329,7 @@ impl eframe::App for App {
                                         points: x_values
                                             .iter()
                                             .zip(y_values_phi0.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name: format!(
                                             "(#{}) phi0, {} ot, {} pt, {} ss, n_o={}, n_p={}",
@@ -331,10 +344,13 @@ impl eframe::App for App {
                                 }
                                 3 => {
                                     let y_values_phi0: Vec<f64> = {
+                                        let bar =
+                                            indicatif::ProgressBar::new(u64::from(self.n_points));
+                                        let t0 = Instant::now();
                                         x_values
                                             .iter()
                                             .map(|&x| {
-                                                let val = calc::phi0(
+                                                        let val = calc::phi0(
                                                     &mut self.calc_cache,
                                                     x,
                                                     self.n_omnes,
@@ -350,7 +366,7 @@ impl eframe::App for App {
                                                 ) / calc::delta_with_lut(
                                                     &mut self.calc_cache,
                                                     x,
-                                                );
+                                                ) - 1.0;
                                                 match self.scaling_type {
                                                     ScalingType::Linear => {
                                                         val
@@ -366,7 +382,7 @@ impl eframe::App for App {
                                         points: x_values
                                             .iter()
                                             .zip(y_values_phi0.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name: format!(
                                             "(#{}) phi0 / delta, {} ot, {} pt, {} ss, n_o={}, n_p={}",
@@ -398,21 +414,52 @@ impl eframe::App for App {
                                         points: x_values
                                             .iter()
                                             .zip(y_values_func.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name:"cont. polynomial".to_string()
                                     });
                                 }
                                 5 => {
-                                    let cutoff_factor = calc::omnes(&mut self.calc_cache, self.cutoff_s, self.n_omnes, self.use_tan_omnes) * self.cutoff_s;
-                                    println!("{:?}", cutoff_factor);
                                     let y_values_func: Vec<f64> = x_values
                                         .iter()
                                         .map(|&x| { 
+                                            let val = if self.use_tan_omnes {
+                                                calc::omnes_integrand_tan(&self.calc_cache, x, 1.2).abs()
+                                            } else {
+                                                calc::omnes_integrand_notan(&self.calc_cache, x, 1.2).abs()
+                                            };
+                                            match self.scaling_type {
+                                                ScalingType::Linear => {
+                                                    val
+                                                }
+                                                ScalingType::Logarithmic => {
+                                                    val.log10()
+                                                }
+                                            }
+                                        }).collect();
+                                    self.plot_data.push(PlotCurve {
+                                        points: x_values
+                                            .iter()
+                                            .zip(y_values_func.iter())
+                                            .map(|(&x, &y)| [x, y])
+                                            .collect(),
+                                        name: format!(
+                                            "(#{}) omnes integrand",
+                                            self.plot_data.len() + 1,
+                                        ),
+                                    })
+                                }
+                                6 => {
+                                    let bar = indicatif::ProgressBar::new(u64::from(self.n_points));
+                                    let cutoff_factor = calc::omnes(&mut self.calc_cache, self.cutoff_s, self.n_omnes, self.use_tan_omnes) * self.cutoff_s;
+                                    let y_values_func: Vec<f64> = x_values
+                                        .iter()
+                                        .map(|&x| { 
+                                            bar.inc(1);
                                             let val = calc::phi0_integrand(
                                                 &mut self.calc_cache, 
                                                 x, 
-                                                2.0, 
+                                                1.0, 
                                                 self.n_omnes, 
                                                 self.use_tan_omnes, 
                                                 self.use_tan_phi0, 
@@ -436,7 +483,7 @@ impl eframe::App for App {
                                         points: x_values
                                             .iter()
                                             .zip(y_values_func.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name: format!(
                                             "(#{}) phi0 integrand",
@@ -444,7 +491,7 @@ impl eframe::App for App {
                                         ),
                                     })
                                 }
-                                6 => {
+                                7 => {
                                     let bar = indicatif::ProgressBar::new(u64::from(self.n_points));
                                     let y_values: Vec<f64> =
                                         x_values.iter().map(|&x| {
@@ -490,7 +537,7 @@ impl eframe::App for App {
                                         points: x_values
                                             .iter()
                                             .zip(y_values.iter())
-                                            .map(|(&x, &y)| [x.powf(0.5), y])
+                                            .map(|(&x, &y)| [x, y])
                                             .collect(),
                                         name: format!(
                                             "(#{}) phi0_cut / phi0, cutoff power: {}",
@@ -499,7 +546,43 @@ impl eframe::App for App {
                                         ),
                                     });
                                 }
-                                7_usize.. => {
+                                8 => {
+                                    let bar = indicatif::ProgressBar::new(u64::from(self.n_points));
+                                    let y_values_func: Vec<f64> = x_values.iter().map(|&x| {
+                                        bar.inc(1);
+                                        let val = calc::phi0(
+                                            &mut self.calc_cache,
+                                            x,
+                                            self.n_omnes,
+                                            self.n_phi0,
+                                            self.use_tan_omnes,
+                                            self.use_tan_phi0,
+                                            self.use_xsub_phi0,
+                                            self.use_reduced_integrand,
+                                            x == self.a,
+                                            self.use_cutoff, 
+                                            self.cutoff_power,
+                                            self.cutoff_s
+                                        ) / calc::delta_with_lut(&mut self.calc_cache, x);
+                                        match self.scaling_type {
+                                            ScalingType::Linear => {
+                                                val
+                                            }
+                                            ScalingType::Logarithmic => {
+                                                val.log10()
+                                            }
+                                        }
+                                    }).collect();
+                                    self.plot_data.push(PlotCurve {
+                                        points: x_values
+                                            .iter()
+                                            .zip(y_values_func.iter())
+                                            .map(|(&x, &y)| [x, y])
+                                            .collect(),
+                                        name:"cont. polynomial".to_string()
+                                    });
+                                }
+                                9_usize.. => {
                                     panic!("This button has no assigned purpose!");
                                 }
                             }
@@ -522,7 +605,7 @@ impl eframe::App for App {
 
             let plot = Plot::new("my_plot")
                 .grid_spacing(5.0f32..=300.0f32)
-                .x_axis_label("sqrt(s) [GeV]")
+                .x_axis_label("s [GeV^2]")
                 .legend(Legend::default());
 
             let plot = match self.scaling_type {
@@ -543,7 +626,7 @@ impl eframe::App for App {
                             if v.value - v.value.floor() < 0.9 {
                                 log_gridmarks.push(GridMark {
                                     value: log_map(v.value),
-                                    step_size: v.step_size / 1.14,
+                                    step_size: v.step_size / 1.14, // This value just looks best
                                 });
                             }
                         }
@@ -561,7 +644,7 @@ impl eframe::App for App {
                     plot_ui.line(
                         Line::new(self.plot_data[i].points.clone())
                             .name(self.plot_data[i].name.clone()),
-                    ); // is clone necessary?
+                    );
                 }
             });
         });
@@ -572,6 +655,7 @@ impl eframe::App for App {
 fn main() {
     let mut app = App::default();
     app.calc_cache = calc::Cache::from_file("./gauss_quad_lut.morello").unwrap();
+    println!("{}", (2.0 * calc::M_K).powi(2));
 
     // let x_values: Vec<f64> = (0..app.n_points)
     //     .map(|x| -> f64 { f64::from(x) * ((app.b - app.a) / f64::from(app.n_points)) + app.a })
@@ -595,12 +679,36 @@ fn main() {
     //     })
     //     .collect();
 
-    let phi0 = calc::phi0(&mut app.calc_cache, 10000.0, 10000, 10000, true, true, true, true, true, false, 0.0, 0.0);
+
-    println!("relative difference: {}", phi0 / calc::delta(10000.0) - 1.0);
+    // let accuracies = [100,500,1000,2000,5000,7500,10000,15000,25000,50000];
+    // let mut results = [[0.0;10];10];
+    // for (o,n_o) in accuracies.iter().enumerate() {
+    //     println!("{:?}", n_o);
+    //     for (p,n_p) in accuracies.iter().enumerate() {
+    //         let phi0 = calc::phi0(&mut app.calc_cache, 10000.0, *n_o, *n_p , true, true, true, false, true, false, 0.0, 0.0);
+    //         results[o][p] = phi0 / calc::delta(10000.0) - 1.0;
+    //     }
+    // }
+
+    // println!("{:#?}", results);
+    
+    let phi0 = calc::phi0(&mut app.calc_cache, 2.0158, 100, 100, true, true, true, false, true, false, 0.0, 0.0);
+    println!("relative difference: {:e}", phi0 / calc::delta_with_lut(&mut app.calc_cache, 2.0158) - 1.0);
+
+    // let A = 1.42;
+    // let numerator = 4.0*A*A*calc::M_K * (calc::delta_with_lut(&mut app.calc_cache, A*A) - PI);
+    // let denominator = calc::LAMBDA1 + 2.0*calc::LAMBDA2 * (0.5*A/calc::M_K - 1.0);
+    // println!("{:?}", - numerator / denominator - A*A);
+
+    let options = eframe::NativeOptions {
+        viewport: egui::ViewportBuilder::default().with_inner_size([1500 as f32, 1000 as f32]),
+        default_theme: Theme::Dark,
+        ..Default::default()
+    };
 
     // eframe::run_native(
     //     "Omnes Calculator",
-    //     eframe::NativeOptions::default(),
+    //     options,
     //     Box::new(|_cc| Box::new(app)),
     // )
     // .unwrap();

ersterTest.py

@@ -23,8 +23,8 @@ lambda_0 = 0 # this is fixed by the low energy parametrization: lambda_0 = delta
 delta1_mid_high = 0 # again, this is fixed by the value of delta1_mid at 1.42 GeV (see later in the code)
 
 # Numerical constants
-gl_n_omnes = 100 # number of points used in the Gauss Legendre Integration
+gl_n_omnes = 10000 # number of points used in the Gauss Legendre Integration
-gl_n_phi0 = 100
+gl_n_phi0 = 10000
 
 gl_lookup = {}
 omnes_lookup = {}
@@ -195,7 +195,7 @@ def phi0(s, a=0, inf = 100000, optimized = True, n = gl_n_phi0):
 	c = np.pi / (s * np.sqrt(s_0)) * np.log(cmath.polar(omnes(s))[0]**2)
 	
 	if optimized == True:
-		integral = integral_gl_tan_reparam(phi0_integrand, s, a, inf, n, False)
+		integral = integral_gl_tan_reparam(phi0_integrand, s, a, inf, n, True)
 	else:
 		integral = integral_gl(phi0_integrand, s, a, inf, n)
 	return -s * np.sqrt(s-s_0) / (2*np.pi) * (integral - c)
@@ -269,7 +269,7 @@ ax.plot(x, y_delta, label = r"$\delta$")
 # print(time() - t0)
 
 t0 = time()
-print(phi0(10000, n=100))
+print(phi0(10000) / delta(10000))
 print(time() - t0)
 
 # ax.legend()
@@ -277,4 +277,4 @@ print(time() - t0)
 # ax.grid(which='major', color='#A4A4A4', linewidth=1)
 # ax.grid(which='minor', color='#B5B5B5', linestyle=':', linewidth=0.5)
 # ax.minorticks_on()
-# plt.show()
+# plt.show()