a working otto cycle

.gitignore

@@ -506,3 +506,4 @@ results/**.npy
 /python/energy_flow_proper/11_new_ho_comparison/ho_data_local/
 /python/energy_flow_proper/11_new_ho_comparison/local_ho_data/
 /python/energy_flow_proper/10_antizeno_engine/taurus/
+/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/taurus/

python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.envrc

@@ -1,2 +1,2 @@
-use_flake
+use_flake --impure
 eval "$shellHook"

python/energy_flow_proper/09_dynamic_two_bath_one_qubit/first_experiments.org

@@ -12,6 +12,7 @@ with the markovian result.
   import qutip as qt
   import utilities as ut
   import stocproc
+  import matplotlib.pyplot as plt
 #+end_src
 
 Init ray and silence stocproc.
@@ -23,7 +24,7 @@ Init ray and silence stocproc.
 #+end_src
 
 #+RESULTS:
-: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '192.168.100.170', 'raylet_ip_address': '192.168.100.170', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-15_13-45-28_971048_59572/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-15_13-45-28_971048_59572/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-15_13-45-28_971048_59572', 'metrics_export_port': 53281, 'gcs_address': '192.168.100.170:61047', 'address': '192.168.100.170:61047', 'node_id': '19344eac29e15d6404c14920662bd8259981954f0a9fde73ef5f8b4b'})
+: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '192.168.100.170', 'raylet_ip_address': '192.168.100.170', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-15_16-33-19_035369_76479/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-15_16-33-19_035369_76479/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-15_16-33-19_035369_76479', 'metrics_export_port': 56879, 'gcs_address': '192.168.100.170:64999', 'address': '192.168.100.170:64999', 'node_id': 'cd3c35a56063d205e53b29934473b89b3baea526868d0ff8e806b08e'})
 
 #+begin_src jupyter-python :results none
   from hops.util.logging_setup import logging_setup
@@ -55,7 +56,7 @@ Let's lay down some basic functionality.
 #+begin_src jupyter-python :results none
   t_compression = .1
   t_expansion = .1
-  comp_ratio = .25
+  comp_ratio = .5
   hot_switch_ratio = .8
   cold_switch_ratio = .8
 
@@ -75,8 +76,8 @@ Let's lay down some basic functionality.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fbff84a3340> |
+| <matplotlib.lines.Line2D | at | 0x7f7678510df0> |
-[[file:./.ob-jupyter/3cc22c46295a7b959768d60bfea0e4c0142e2b43.svg]]
+[[file:./.ob-jupyter/8326a30a0b67dac1fcbfa726d3677ad081f50681.svg]]
 :END:
 
 ** Hot Thermalization
@@ -89,8 +90,8 @@ Let's lay down some basic functionality.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fbff840e1c0> |
+| <matplotlib.lines.Line2D | at | 0x7f7669dd3b50> |
-[[file:./.ob-jupyter/98aec1955458432d8e29d21e1675b698b47b50fb.svg]]
+[[file:./.ob-jupyter/beef0051823ad8d69a03749d66188767c06f00d6.svg]]
 :END:
 
 ** Compression
@@ -101,8 +102,8 @@ Let's lay down some basic functionality.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fbff83eeb20> |
+| <matplotlib.lines.Line2D | at | 0x7f7669d4e100> |
-[[file:./.ob-jupyter/8ebe76752707b723232c25684baa8b2b3a0af60e.svg]]
+[[file:./.ob-jupyter/2c1c45c42776db17721492efacd6f4856f8fe3b3.svg]]
 :END:
 
 ** Cold Thermalization
@@ -115,8 +116,8 @@ Let's lay down some basic functionality.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fbfddd4deb0> |
+| <matplotlib.lines.Line2D | at | 0x7f7669d33610> |
-[[file:./.ob-jupyter/ffdfdc73f371b6485c38c7714e3d6f2075350606.svg]]
+[[file:./.ob-jupyter/a425435a0d34373dd92cf2a415ceec18a1fd69af.svg]]
 :END:
 
 ** Full Cycle
@@ -125,29 +126,73 @@ Now we turn the system around after each fill-cycle.
   H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_expansion, 1]), 1)
   L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)]
   period = 1
-  scale = .1
+  scale = .05
   H_cyc, *L = [ScaleTime(op, scale) for op in [H_cyc, *L]]
 
   ω_mod = 2*np.pi / (period) * scale
-  plt.plot(t, H_cyc.operator_norm(t))
-  plt.plot(t, L[0].operator_norm(t))
-  plt.plot(t, L[1].operator_norm(t))
   ω_mod
   t_cycle = period/scale
-  periods = 20
+  periods = 10
   t_max = t_cycle * periods
-  dt = .01
+  dt = .1
+
+  t = np.linspace(0, t_cycle, 1000)
+
+  plt.plot(t, H_cyc.operator_norm(t), label=r"$||H_\mathrm{S}||$")
+  plt.plot(t, L[0].operator_norm(t), label=r"$||L_\mathrm{c}||$")
+  plt.plot(t, L[1].operator_norm(t), label=r"$||L_\mathrm{h}||$")
+  plt.xlabel(r"$\tau$")
+  plt.legend()
+  fs.export_fig("modulation")
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/c6fcbae6c5d076e45f48b0d4e003265146c48409.svg]]
+[[file:./.ob-jupyter/2b79808a3127f0a93ce6cbde00ae4888c83b2696.svg]]
+
+** Shifts
+#+begin_src jupyter-python
+  from scipy.optimize import minimize_scalar
+
+  ref_model = QubitModelMutliBath(
+      ω_c=[1] * 2,
+      T=[1, 20],
+  )
+
+  shifts = []
+  scales = []
+
+  for i, t_exp in enumerate([0, t_expansion / scale]):
+      ω_exp = H_cyc.operator_norm(i)
+
+      def objective(ω_s):
+          ref_model.ω_s[i] = ω_s
+
+          return -ref_model.full_thermal_spectral_density(i)(ω_exp)
+
+      res = minimize_scalar(
+          objective, 1, method="bounded", bounds=(0.01, ω_exp), options=dict(maxiter=100)
+      )
+      shifts.append(res.x)
+      scales.append(-res.fun)
+      ω = np.linspace(0.1, 10, 1000)
+      plt.plot(ω, ref_model.full_thermal_spectral_density(i)(ω) / -res.fun)
+
+  print(shifts, scales)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+: [0.010005778765750251, 0.5000003902186526] [0.3656483223378975, 3.198956858189882]
+[[file:./.ob-jupyter/97d3960ba8288d117b1bf5774e57ac8bb1500b83.svg]]
+:END:
+
 
 * Model
 #+begin_src jupyter-python :results none
   model = QubitModelMutliBath(
-      δ=[.05*4, .01*4],
+      δ=list(1/np.array(scales) * .1/2),#[.05*4, .01*4],
-      ω_c=[1] * 2,
+      ω_c=ref_model.ω_c,
-      ω_s=[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1],
+      ω_s=shifts,#[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1],
       t=ut.linspace_with_strobe(0, t_max, int(t_max // dt), ω_mod),
       ψ_0=(qt.basis([2], [0])*0  + qt.basis([2], [1])*1),
       description=f"A simple otto cycle.",
@@ -156,95 +201,60 @@ Now we turn the system around after each fill-cycle.
       truncation_scheme="simplex",
       driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
       thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
-      T = [1, 20],
+      T = ref_model.T, #= [1, 20],
       L = L,
       H = H_cyc,
-      therm_methods=["tanhsinh", "fft"],
+      therm_methods=["fft", "fft"],
   )
 #+end_src
 
 
 #+begin_src jupyter-python
-  aux.integrate(model, 1000)
+  aux.integrate(model, 10000)
 #+end_src
 
 #+RESULTS:
-: [INFO    hops.core.integration     59572] Choosing the nonlinear integrator.
+: [INFO    hops.core.integration     76479] Choosing the nonlinear integrator.
-: [INFO    hops.core.integration     59572] Using 8 integrators.
+: [INFO    hops.core.integration     76479] Using 8 integrators.
-: [INFO    hops.core.integration     59572] Some 0 trajectories have to be integrated.
+: [INFO    hops.core.integration     76479] Some 30 trajectories have to be integrated.
-: [INFO    hops.core.integration     59572] Using 1820 hierarchy states.
+: [INFO    hops.core.integration     76479] Using 1820 hierarchy states.
-: 0it [00:00, ?it/s]
+: 100% 30/30 [00:22<00:00,  1.32it/s]
 
 
 #+begin_src jupyter-python
-  f, a = pu.plot_energy_overview(model)
+  f, a = pu.plot_energy_overview(model, strobe_frequency=ω_mod)
-  a.plot(model.t, model.H.operator_norm(model.t))
+  # a.plot(model.t, model.H.operator_norm(model.t))
-  a.plot(model.t, model.L[0].operator_norm(model.t))
+  # a.plot(model.t, model.L[0].operator_norm(model.t))
-  a.plot(model.t, model.L[1].operator_norm(model.t))
+  # a.plot(model.t, model.L[1].operator_norm(model.t))
+  ax.set_xlabel(r"$\tau$")
   a.legend()
+  #fs.export_fig("energy_strobe", y_scaling=.6)
 #+end_src
 
 #+RESULTS:
-:RESULTS:
+[[file:./.ob-jupyter/5d6c6959a02462528cdb3d541f8dd5ecf303f362.svg]]
-#+begin_example
-  Loading: 100% 8/8 [00:00<00:00, 192.76it/s]
-  [INFO    root                      59572] Writing cache to: results/5e8c8c2b39bd3b8ad75fb0da753ca3989cec2ca758e5b60d542cfb0c39e6009d_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_op_exp_task_33_None_7b4b9e9cfa38e0c95465b4087684ba9103d527ed24c440718c98f01d3693d041.npy
-  Loading: 100% 16/16 [00:00<00:00, 50.04it/s]
-  [INFO    root                      59572] Writing cache to: results/flow_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_flow_worker_33_None_c15d067b912279587be3d3939223173fe5368154834a5a60648fc0a46950f8cf.npy
-  Loading: 100% 16/16 [00:00<00:00, 86.42it/s]
-  [INFO    root                      59572] Writing cache to: results/interaction_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_interaction_energy_task_33_None_58455bf5ce6e35b96d04516e2a659a907e0c44b588167e00faa3fdf3e226bb11.npy
-  Loading: 100% 16/16 [00:00<00:00, 79.44it/s]
-  [INFO    root                      59572] Writing cache to: results/interaction_power_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_interaction_energy_task_33_None_58455bf5ce6e35b96d04516e2a659a907e0c44b588167e00faa3fdf3e226bb11.npy
-  Loading: 100% 8/8 [00:00<00:00, 245.72it/s]
-  [INFO    root                      59572] Writing cache to: results/143e514d782a0bf3369be86475c494b0df4baa9952cf06bbf576d262f2dd27a4_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_op_exp_task_33_None_7b4b9e9cfa38e0c95465b4087684ba9103d527ed24c440718c98f01d3693d041.npy
-#+end_example
-: <matplotlib.legend.Legend at 0x7fbfc42bd580>
-[[file:./.ob-jupyter/229ace105ac0a848d7b968ebdfd79d9fd38dbe9e.svg]]
-:END:
 
 #+begin_src jupyter-python
   fig, ax = plt.subplots()
   with aux.get_data(model) as data:
-        power = model.interaction_power(data).for_bath(0)
+        power = model.interaction_power(data).sum_baths()
         power_sys = model.system_power(data).sum_baths()
         energy = model.total_energy_from_power(data)
-        pu.plot_with_σ(model.t, power, ax=ax)
+        pu.plot_with_σ(model.t, power, ax=ax, label="Interaction")
-        pu.plot_with_σ(model.t, power_sys, ax=ax, label="sys")
+        pu.plot_with_σ(model.t, power_sys, ax=ax, label="System")
-        pu.plot_with_σ(model.t, energy, ax=ax)
 
   ax.legend()
+  ax.set_title("Power")
+  ax.set_xlabel(r"$\tau$")
+  fs.export_fig("power", y_scaling=.4)
 #+end_src
 
 #+RESULTS:
-:RESULTS:
+[[file:./.ob-jupyter/4938909f0a3699e7c9912ab201d2f04199b355e7.svg]]
-# [goto error]
-#+begin_example
-  ---------------------------------------------------------------------------
-  RuntimeError                              Traceback (most recent call last)
-  Input In [97], in <cell line: 2>()
-        1 fig, ax = plt.subplots()
-  ----> 2 with aux.get_data(model) as data:
-        3       power = model.interaction_power(data).for_bath(0)
-        4       power_sys = model.system_power(data).sum_baths()
-
-  File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:146, in get_data(model, data_path, read_only, **kwargs)
-      135         return HIData(
-      136             path,
-      137             hi_key=model.hops_config,
-     (...)
-      142             **kwargs,
-      143         )
-      145 else:
-  --> 146     raise RuntimeError(f"No data found for model with hash '{hexhash}'.")
-
-  RuntimeError: No data found for model with hash '2676c736eb2bf1d2a1f37e21cda1823a3dc6b3f2cedb82509cfa6231a5232b2f'.
-#+end_example
-[[file:./.ob-jupyter/5285f6b51867d45830747c7d1501193929fe8fd1.svg]]
-:END:
 
 
 #+begin_src jupyter-python
-  steady_index = 6
+  steady_index = 3
   cycle_times, cycle_indices = ut.strobe_times(model.t, ω_mod)
 
 
@@ -257,7 +267,7 @@ Now we turn the system around after each fill-cycle.
       )
 
   e_diff = total_energies.slice(0) - total_energies.slice(-1)
-  mean_power = e_diff * (1 / (model.t[-1] - model.t[steady_index]))
+  mean_power = e_diff * (1 / (model.t[-1] - cycle_times[steady_index]))
   efficiency = e_diff * (
       1 / (hot_bath_energies.slice(0) - hot_bath_energies.slice(-1)).value
   )
@@ -265,13 +275,13 @@ Now we turn the system around after each fill-cycle.
 #+end_src
 
 #+RESULTS:
-| -0.011704395396932077 | 8.941753951852788e-05 | -0.6487259446355118 |
+| 0.0019806021074448648 | 2.3085106663538818e-05 | 0.1975630167587971 |
 
 #+begin_src jupyter-python
   with aux.get_data(model) as data:
       ρ = data.rho_t_accum.mean[:]
       σ_ρ = data.rho_t_accum.ensemble_std[:]
-      plt.plot(model.t, ρ[:, 1, 1])
+      plt.plot(model.t, ρ[:, 0, 0])
       xs = np.einsum("tij,ji->t", ρ, qt.sigmax().full())
       ys = np.einsum("tij,ji->t", ρ, qt.sigmay().full())
       zs = np.einsum("tij,ji->t", ρ, qt.sigmaz().full())
@@ -279,9 +289,9 @@ Now we turn the system around after each fill-cycle.
 
 #+RESULTS:
 :RESULTS:
-: /nix/store/dbz51ji4p5pcppra5h7yqhwcrbqgq40v-python3-3.9.13-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part
+: /nix/store/yf7vf1pic6xfwr035xywcczjm08psvpa-python3-3.9.13-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part
 :   return np.asarray(x, float)
-[[file:./.ob-jupyter/9b6dbfd5d863d78fc20ce574cf7e282bc8341b0c.svg]]
+[[file:./.ob-jupyter/69c7da5a4fa9adc0fb49b0ff097ea6fbc28c58ca.svg]]
 :END:
 
 
@@ -299,9 +309,9 @@ Now we turn the system around after each fill-cycle.
 
 #+RESULTS:
 :RESULTS:
-: /nix/store/dbz51ji4p5pcppra5h7yqhwcrbqgq40v-python3-3.9.13-env/lib/python3.9/site-packages/qutip/bloch.py:639: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword
+: /nix/store/yf7vf1pic6xfwr035xywcczjm08psvpa-python3-3.9.13-env/lib/python3.9/site-packages/qutip/bloch.py:639: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword
 :   self.fig.canvas.draw()
-: /nix/store/dbz51ji4p5pcppra5h7yqhwcrbqgq40v-python3-3.9.13-env/lib/python3.9/site-packages/IPython/core/pylabtools.py:151: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword
+: /nix/store/yf7vf1pic6xfwr035xywcczjm08psvpa-python3-3.9.13-env/lib/python3.9/site-packages/IPython/core/pylabtools.py:151: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword
 :   fig.canvas.print_figure(bytes_io, **kw)
-[[file:./.ob-jupyter/7806ea3a984dc21f51c7f09f8720ae8e361465a5.svg]]
+[[file:./.ob-jupyter/633b45a276b59c015a575e42b611b52283311ebe.svg]]
 :END:

python/energy_flow_proper/09_dynamic_two_bath_one_qubit/mount_taurus.sh

@@ -0,0 +1,7 @@
+#!/usr/bin/env bash
+sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus  s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_11/project/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.data .data
+sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus  s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_11/project/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/results results
+
+
+
+sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus  s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_11/project/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/  taurus

python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor

@@ -0,0 +1,125 @@
+import figsaver as fs
+import plot_utils as pu
+from hiro_models.one_qubit_model import QubitModelMutliBath, StocProcTolerances
+import hiro_models.model_auxiliary as aux
+import numpy as np
+import qutip as qt
+import utilities as ut
+import stocproc
+
+import ray
+ray.shutdown()
+#ray.init(address='auto')
+ray.init()
+
+from hops.util.logging_setup import logging_setup
+import logging
+logging_setup(logging.INFO)
+
+from hops.util.dynamic_matrix import SmoothStep, Periodic, ConstantMatrix, ScaleTime, Shift, Piecewise
+
+H_op = 2*(1/2 * (qt.sigmaz() + qt.identity(2))).full()
+L_op = (1/2 * qt.sigmax()).full()
+L_op_2 = (1/2 * qt.sigmax()).full()
+print(H_op, L_op, sep="\n")
+
+t_compression = .1
+t_expansion = .1
+comp_ratio = .5
+hot_switch_ratio = .8
+cold_switch_ratio = .8
+
+t_hot = (1 - t_expansion - t_compression) / 2
+t_hot_switch = t_hot * hot_switch_ratio / 2
+t_cold = (1 - t_expansion - t_compression) / 2
+t_cold_switch = t_cold * cold_switch_ratio / 2
+t = np.linspace(0, 10, 1000)
+
+%matplotlib inline
+H_exp = SmoothStep(comp_ratio * H_op, 0, t_expansion, 2) + H_op * (1 - comp_ratio)
+tt = np.linspace(0, 1, 1000)
+plt.plot(tt, H_exp.operator_norm(tt))
+
+L_hot = SmoothStep(L_op, t_expansion, t_expansion + t_hot_switch, 2) - SmoothStep(
+      L_op, t_expansion + t_hot - t_hot_switch, t_expansion + t_hot, 2
+  )
+plt.plot(tt, L_hot.operator_norm(tt))
+
+H_comp = ConstantMatrix(H_op) - SmoothStep(comp_ratio * H_op, t_expansion + t_hot, t_expansion + t_hot + t_compression, 2)
+plt.plot(tt, H_comp.operator_norm(tt))
+
+L_cold = SmoothStep(L_op_2, t_expansion + t_hot + t_compression, t_expansion + t_hot + t_compression + t_cold_switch, 3) - SmoothStep(
+    L_op_2, t_expansion + t_hot + t_compression + t_cold - t_cold_switch, 1, 2
+)
+plt.plot(tt, L_cold.operator_norm(tt))
+
+H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_expansion, 1]), 1)
+L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)]
+period = 1
+scale = .1
+H_cyc, *L = [ScaleTime(op, scale) for op in [H_cyc, *L]]
+
+ω_mod = 2*np.pi / (period) * scale
+plt.plot(t, H_cyc.operator_norm(t))
+plt.plot(t, L[0].operator_norm(t))
+plt.plot(t, L[1].operator_norm(t))
+ω_mod
+t_cycle = period/scale
+periods = 20
+t_max = t_cycle * periods
+dt = .01
+
+model = QubitModelMutliBath(
+    δ=[.5, .1],
+    ω_c=[1] * 2,
+    ω_s=[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1],
+    t=ut.linspace_with_strobe(0, t_max, int(t_max // dt), ω_mod),
+    ψ_0=(qt.basis([2], [0])*0  + qt.basis([2], [1])*1),
+    description=f"A simple otto cycle.",
+    k_max=4,
+    bcf_terms=[6]*2,
+    truncation_scheme="simplex",
+    driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
+    thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
+    T = [0, 20],
+    L = L,
+    H = H_cyc,
+    therm_methods=["tanhsinh", "fft"],
+)
+
+aux.integrate(model, 200)
+
+%matplotlib inline
+f, a = pu.plot_energy_overview(model)
+a.plot(model.t, model.H.operator_norm(model.t))
+a.plot(model.t, model.L[0].operator_norm(model.t))
+a.plot(model.t, model.L[1].operator_norm(model.t))
+a.legend()
+
+fig, ax = plt.subplots()
+with aux.get_data(model) as data:
+      power = model.interaction_power(data).sum_baths()
+      power_sys = model.system_power(data).sum_baths()
+      energy = model.total_energy_from_power(data).sum_baths()
+      pu.plot_with_σ(model.t, power, ax=ax)
+      pu.plot_with_σ(model.t, power_sys, ax=ax, label="sys")
+      pu.plot_with_σ(model.t, energy, ax=ax)
+
+ax.legend()
+
+with aux.get_data(model) as data:
+    ρ = data.rho_t_accum.mean[:]
+    σ_ρ = data.rho_t_accum.ensemble_std[:]
+    plt.plot(model.t, ρ[:, 1, 1])
+    xs = np.einsum("tij,ji->t", ρ, qt.sigmax().full())
+    ys = np.einsum("tij,ji->t", ρ, qt.sigmay().full())
+    zs = np.einsum("tij,ji->t", ρ, qt.sigmaz().full())
+
+b = qt.Bloch()
+b.add_points([xs, ys, zs])
+b.view = [20, 20]
+b.point_size = [0.1]
+b.sphere_alpha = 0.1
+b.size = [20, 20]
+b.render()
+b.fig

python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor.py

@@ -6,6 +6,7 @@ import numpy as np
 import qutip as qt
 import utilities as ut
 import stocproc
+import matplotlib.pyplot as plt
 
 import ray
 ray.shutdown()
@@ -25,7 +26,7 @@ print(H_op, L_op, sep="\n")
 
 t_compression = .1
 t_expansion = .1
-comp_ratio = .25
+comp_ratio = .5
 hot_switch_ratio = .8
 cold_switch_ratio = .8
 
@@ -55,23 +56,57 @@ plt.plot(tt, L_cold.operator_norm(tt))
 H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_expansion, 1]), 1)
 L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)]
 period = 1
-scale = .1
+scale = .05
 H_cyc, *L = [ScaleTime(op, scale) for op in [H_cyc, *L]]
 
 ω_mod = 2*np.pi / (period) * scale
-plt.plot(t, H_cyc.operator_norm(t))
-plt.plot(t, L[0].operator_norm(t))
-plt.plot(t, L[1].operator_norm(t))
 ω_mod
 t_cycle = period/scale
-periods = 20
+periods = 10
 t_max = t_cycle * periods
-dt = .01
+dt = .1
+
+t = np.linspace(0, t_cycle, 1000)
+
+plt.plot(t, H_cyc.operator_norm(t), label=r"$||H_\mathrm{S}||$")
+plt.plot(t, L[0].operator_norm(t), label=r"$||L_\mathrm{c}||$")
+plt.plot(t, L[1].operator_norm(t), label=r"$||L_\mathrm{h}||$")
+plt.xlabel(r"$\tau$")
+plt.legend()
+fs.export_fig("modulation")
+
+from scipy.optimize import minimize_scalar
+
+ref_model = QubitModelMutliBath(
+    ω_c=[1] * 2,
+    T=[1, 20],
+)
+
+shifts = []
+scales = []
+
+for i, t_exp in enumerate([0, t_expansion / scale]):
+    ω_exp = H_cyc.operator_norm(i)
+
+    def objective(ω_s):
+        ref_model.ω_s[i] = ω_s
+
+        return -ref_model.full_thermal_spectral_density(i)(ω_exp)
+
+    res = minimize_scalar(
+        objective, 1, method="bounded", bounds=(0.01, ω_exp), options=dict(maxiter=100)
+    )
+    shifts.append(res.x)
+    scales.append(-res.fun)
+    ω = np.linspace(0.1, 10, 1000)
+    plt.plot(ω, ref_model.full_thermal_spectral_density(i)(ω) / -res.fun)
+
+print(shifts, scales)
 
 model = QubitModelMutliBath(
-    δ=[.05*4, .01*4],
+    δ=list(1/np.array(scales) * .1/2),#[.05*4, .01*4],
-    ω_c=[1] * 2,
+    ω_c=ref_model.ω_c,
-    ω_s=[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1],
+    ω_s=shifts,#[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1],
     t=ut.linspace_with_strobe(0, t_max, int(t_max // dt), ω_mod),
     ψ_0=(qt.basis([2], [0])*0  + qt.basis([2], [1])*1),
     description=f"A simple otto cycle.",
@@ -80,32 +115,36 @@ model = QubitModelMutliBath(
     truncation_scheme="simplex",
     driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
     thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
-    T = [1, 20],
+    T = ref_model.T, #= [1, 20],
     L = L,
     H = H_cyc,
-    therm_methods=["tanhsinh", "fft"],
+    therm_methods=["fft", "fft"],
 )
 
-aux.integrate(model, 1000)
+aux.integrate(model, 10000)
 
-f, a = pu.plot_energy_overview(model)
+f, a = pu.plot_energy_overview(model, strobe_frequency=ω_mod)
-a.plot(model.t, model.H.operator_norm(model.t))
+# a.plot(model.t, model.H.operator_norm(model.t))
-a.plot(model.t, model.L[0].operator_norm(model.t))
+# a.plot(model.t, model.L[0].operator_norm(model.t))
-a.plot(model.t, model.L[1].operator_norm(model.t))
+# a.plot(model.t, model.L[1].operator_norm(model.t))
+ax.set_xlabel(r"$\tau$")
 a.legend()
+#fs.export_fig("energy_strobe", y_scaling=.6)
 
 fig, ax = plt.subplots()
 with aux.get_data(model) as data:
-      power = model.interaction_power(data).for_bath(0)
+      power = model.interaction_power(data).sum_baths()
       power_sys = model.system_power(data).sum_baths()
       energy = model.total_energy_from_power(data)
-      pu.plot_with_σ(model.t, power, ax=ax)
+      pu.plot_with_σ(model.t, power, ax=ax, label="Interaction")
-      pu.plot_with_σ(model.t, power_sys, ax=ax, label="sys")
+      pu.plot_with_σ(model.t, power_sys, ax=ax, label="System")
-      pu.plot_with_σ(model.t, energy, ax=ax)
 
 ax.legend()
+ax.set_title("Power")
+ax.set_xlabel(r"$\tau$")
+fs.export_fig("power", y_scaling=.4)
 
-steady_index = 6
+steady_index = 3
 cycle_times, cycle_indices = ut.strobe_times(model.t, ω_mod)
 
 
@@ -118,7 +157,7 @@ with aux.get_data(model) as data:
     )
 
 e_diff = total_energies.slice(0) - total_energies.slice(-1)
-mean_power = e_diff * (1 / (model.t[-1] - model.t[steady_index]))
+mean_power = e_diff * (1 / (model.t[-1] - cycle_times[steady_index]))
 efficiency = e_diff * (
     1 / (hot_bath_energies.slice(0) - hot_bath_energies.slice(-1)).value
 )
@@ -127,7 +166,7 @@ mean_power.value, mean_power.σ, efficiency.value
 with aux.get_data(model) as data:
     ρ = data.rho_t_accum.mean[:]
     σ_ρ = data.rho_t_accum.ensemble_std[:]
-    plt.plot(model.t, ρ[:, 1, 1])
+    plt.plot(model.t, ρ[:, 0, 0])
     xs = np.einsum("tij,ji->t", ρ, qt.sigmax().full())
     ys = np.einsum("tij,ji->t", ρ, qt.sigmay().full())
     zs = np.einsum("tij,ji->t", ρ, qt.sigmaz().full())