Add hamilton function

6 files changed, 101 insertions, 38 deletions

diff --git a/Cargo.lock b/Cargo.lock
index 61a4282..b8033a2 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -31,6 +31,16 @@ version = "1.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 
 [[package]]
+name = "blas"
+version = "0.19.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "blas-sys 0.7.1 (registry+https://github.com/rust-lang/crates.io-index)",
+ "libc 0.2.66 (registry+https://github.com/rust-lang/crates.io-index)",
+ "num-complex 0.1.43 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
+[[package]]
 name = "blas-src"
 version = "0.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -41,6 +51,14 @@ version = "0.4.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 
 [[package]]
+name = "blas-sys"
+version = "0.7.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "libc 0.2.66 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
+[[package]]
 name = "byteorder"
 version = "1.3.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -178,11 +196,13 @@ dependencies = [
 name = "dequantify"
 version = "0.1.0"
 dependencies = [
+ "blas 0.19.1 (registry+https://github.com/rust-lang/crates.io-index)",
  "expm 0.1.1 (registry+https://github.com/rust-lang/crates.io-index)",
  "ndarray 0.13.0 (registry+https://github.com/rust-lang/crates.io-index)",
  "ndarray-linalg 0.12.0 (registry+https://github.com/rust-lang/crates.io-index)",
  "num-complex 0.2.3 (registry+https://github.com/rust-lang/crates.io-index)",
  "ocl 0.19.3 (registry+https://github.com/rust-lang/crates.io-index)",
+ "openblas-src 0.5.6 (registry+https://github.com/rust-lang/crates.io-index)",
 ]
 
 [[package]]
@@ -486,6 +506,11 @@ dependencies = [
 ]
 
 [[package]]
+name = "openblas-src"
+version = "0.5.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+
+[[package]]
 name = "ordered-float"
 version = "1.0.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
@@ -777,8 +802,10 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 "checksum backtrace 0.3.40 (registry+https://github.com/rust-lang/crates.io-index)" = "924c76597f0d9ca25d762c25a4d369d51267536465dc5064bdf0eb073ed477ea"
 "checksum backtrace-sys 0.1.32 (registry+https://github.com/rust-lang/crates.io-index)" = "5d6575f128516de27e3ce99689419835fce9643a9b215a14d2b5b685be018491"
 "checksum bitflags 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)" = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693"
+"checksum blas 0.19.1 (registry+https://github.com/rust-lang/crates.io-index)" = "cd6a9d7b0bfc50e540a216db5517248163572d9f9a502e47b64d5a97236e4314"
 "checksum blas-src 0.2.1 (registry+https://github.com/rust-lang/crates.io-index)" = "a7b894d1cc14af76750a80f64387f7986fe93a9d3786a9a8926a340ae50b2c51"
 "checksum blas-src 0.4.0 (registry+https://github.com/rust-lang/crates.io-index)" = "1c18550bb5de993fd6ba28b74b05b0aadd693002e566c90cec65636635f570de"
+"checksum blas-sys 0.7.1 (registry+https://github.com/rust-lang/crates.io-index)" = "13b1b279ceb25d7c4faaea95a5f7addbe7d8c34f9462044bd8e630cebcfc2440"
 "checksum byteorder 1.3.2 (registry+https://github.com/rust-lang/crates.io-index)" = "a7c3dd8985a7111efc5c80b44e23ecdd8c007de8ade3b96595387e812b957cf5"
 "checksum cauchy 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)" = "9c505375a294ec7cf79430c4ccb91ee1fab5f6390da3a264932fd303832a2de7"
 "checksum cblas 0.2.0 (registry+https://github.com/rust-lang/crates.io-index)" = "d82f331add33eceb4c41cb28d878049b96f56577016daf190831e94e4aece5db"
@@ -827,6 +854,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 "checksum ocl 0.19.3 (registry+https://github.com/rust-lang/crates.io-index)" = "fa54742f44f5813a4cb4c8c16f7b4069cb3a17a88a354754a50d6b951ae32e31"
 "checksum ocl-core 0.11.2 (registry+https://github.com/rust-lang/crates.io-index)" = "81bc628faf959b5e07b1251252926dfe0dd1b3f2709cef8998c97936ddbdaa74"
 "checksum ocl-core-vector 0.1.0 (registry+https://github.com/rust-lang/crates.io-index)" = "b4072920739958adeec5abedec51af70febc58f7fff0601aaa0827c1f3c8fefd"
+"checksum openblas-src 0.5.6 (registry+https://github.com/rust-lang/crates.io-index)" = "68d293fca3c73ad377ddd2236d32c828a50a611a5b472bf6a884b9b60a4acd97"
 "checksum ordered-float 1.0.2 (registry+https://github.com/rust-lang/crates.io-index)" = "18869315e81473c951eb56ad5558bbc56978562d3ecfb87abb7a1e944cea4518"
 "checksum proc-macro2 1.0.6 (registry+https://github.com/rust-lang/crates.io-index)" = "9c9e470a8dc4aeae2dee2f335e8f533e2d4b347e1434e5671afc49b054592f27"
 "checksum quote 1.0.2 (registry+https://github.com/rust-lang/crates.io-index)" = "053a8c8bcc71fcce321828dc897a98ab9760bef03a4fc36693c231e5b3216cfe"
diff --git a/Cargo.toml b/Cargo.toml
index 2b90ded..51776da 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -6,6 +6,8 @@ authors = ["Dennis Kobert <dennis@kobert.dev>",
 edition = "2018"
 
 [dependencies]
+blas = "0.19"
+openblas-src = "0.5"
 ocl = "0.19"
 ndarray-linalg = "0.12"
 num-complex = "0.2"
diff --git a/src/main.rs b/src/main.rs
index ab88302..773ea72 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,44 +1,42 @@
 pub mod plot;
 mod simulation;
+use num_complex::Complex;
 use plot::Graph;
+use simulation::time_evolution::*;
+use simulation::*;
+
+struct StaticHamiltonion;
+
+impl MatrixGen for StaticHamiltonion {
+    fn gen(_t: f64) -> ndarray::Array2<c64> {
+        let mut b = unsafe { ndarray::Array2::<c64>::uninitialized((4, 4)) };
+        expm::expm(&simulation::hamiltonians::X0(), &mut b);
+        b
+    }
+}
 
 fn main() {
-    println!("Hello, world!");
-    plot::Gnuplot::plot3d(
-        &[
-            [
-                0.14285714285714296,
-                0.5714285714285717,
-                0.0,
-                0.5714285714285717,
-                0.5714285714285717,
-                0.0,
-            ],
-            [
-                0.11910591478146304,
-                0.5652016093730936,
-                -0.01707887122942119,
-                0.6464319082591804,
-                0.498191752207936,
-                0.0002,
-            ],
-            [
-                0.21910591478146304,
-                0.7652016093730936,
-                -0.01707887122942119,
-                0.6464319082591804,
-                0.498191752207936,
-                0.0002,
-            ],
-            [
-                0.01910591478146304,
-                0.8652016093730936,
-                -0.01707887122942119,
-                0.6464319082591804,
-                0.498191752207936,
-                0.0002,
-            ],
-        ],
-        1,
-    )
+    let dt = 1.0;
+    let iterations = 100;
+
+    let state = simulation::twoxtwo_level::TwoXTwoLevel::new(
+        Complex::new(0.0, 0.0),
+        Complex::new(0.0, 0.0),
+        Complex::new(0.0, 0.0),
+        Complex::new(0.0, 0.0),
+    );
+    let trace = Vec::new();
+    for i in 0..iterations {
+        state.evolve::<StaticHamiltonion>(i as f64 * dt);
+        let fibre = state.fibrate();
+        trace.push([
+            fibre[0],
+            fibre[1],
+            fibre[2],
+            0.0,
+            0.0,
+            i as f64 / iterations as f64,
+        ]);
+    }
+    plot::Gnuplot::plot3d(&trace, 1)
 }
diff --git a/src/simulation/hamiltonians.rs b/src/simulation/hamiltonians.rs
new file mode 100644
index 0000000..03d8d1f
--- /dev/null
+++ b/src/simulation/hamiltonians.rs
@@ -0,0 +1,31 @@
+use super::c64;
+use ndarray::prelude::*;
+
+pub fn X0() -> Array2<c64> {
+    array![
+        [
+            c64::new(1.0, 0.0),
+            c64::new(0.0, 0.0),
+            c64::new(0.0, 0.0),
+            c64::new(0.0, 0.0),
+        ],
+        [
+            c64::new(0.0, 0.0),
+            c64::new(1.0, 0.0),
+            c64::new(0.0, 0.0),
+            c64::new(0.0, 0.0),
+        ],
+        [
+            c64::new(0.0, 0.0),
+            c64::new(0.0, 0.0),
+            c64::new(-1.0, 0.0),
+            c64::new(0.0, 0.0),
+        ],
+        [
+            c64::new(0.0, 0.0),
+            c64::new(0.0, 0.0),
+            c64::new(0.0, 0.0),
+            c64::new(-1.0, 0.0),
+        ]
+    ]
+}
diff --git a/src/simulation/mod.rs b/src/simulation/mod.rs
index 641444e..c72f275 100644
--- a/src/simulation/mod.rs
+++ b/src/simulation/mod.rs
@@ -1,3 +1,4 @@
+pub mod hamiltonians;
 pub mod s_5_fibration;
 pub mod time_evolution;
 pub mod two_level;
diff --git a/src/simulation/twoxtwo_level.rs b/src/simulation/twoxtwo_level.rs
index 8a34167..62425c8 100644
--- a/src/simulation/twoxtwo_level.rs
+++ b/src/simulation/twoxtwo_level.rs
@@ -25,11 +25,14 @@ impl State for TwoXTwoLevel {
         for val in conjugate.iter_mut() {
             *val = val.conj();
         }
+        let v = conjugate.dot(&X3().dot(&self.state))[(0, 0)];
 
         vec![
             conjugate.dot(&X0().dot(&self.state))[(0, 0)].re,
             conjugate.dot(&X1().dot(&self.state))[(0, 0)].re,
             conjugate.dot(&X2().dot(&self.state))[(0, 0)].re,
+            v.re,
+            v.im,
         ]
     }