backup before big change

7 files changed, 122 insertions, 1 deletions

diff --git a/.gitignore b/.gitignore
index 155bfa1..4cfe30a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,6 @@ neurobench_testing/
 *.hdf5
 #AI code assistant
 .aider*
+
+#spires reservoir
+spires/
diff --git a/MIL_testing/hp_memristor.cir b/MIL_testing/hp_memristor.cir
new file mode 100644
index 0000000..1865a7a
--- /dev/null
+++ b/MIL_testing/hp_memristor.cir
@@ -0,0 +1,45 @@
+* HP Memristor SPICE Model
+* For Transient Analysis only
+* created by Zdenek and Dalibor Biolek
+**************************
+* Ron, Roff - Resistance in ON / OFF States
+* Rinit - Resistance at T=0
+* D - Width of the thin film
+* uv - Migration coefficient
+* p - Parameter of the WINDOW-function
+* for modeling nonlinear boundary conditions
+* x - W/D Ratio, W is the actual width
+* of the doped area (from 0 to D)
+*
+.SUBCKT memristor Plus Minus PARAMS:
++ Ron=1K Roff=100K Rinit=80K D=10N uv=10F p=1
+***********************************************
+* DIFFERENTIAL EQUATION MODELING *
+***********************************************
+Gx 0 x value={ I(Emem)*uv*Ron/D^2*f(V(x),p)}
+Cx x 0 1 IC={(Roff-Rinit)/(Roff-Ron)}
+Raux x 0 1T
+
+* RESISTIVE PORT OF THE MEMRISTOR *
+*******************************
+Emem plus aux value={-I(Emem)*V(x)*(Roff-Ron)}
+Roff aux minus {Roff}
+***********************************************
+*Flux computation*
+***********************************************
+*does not work with ngspice
+*Eflux flux 0 value={SDT(V(plus,minus))}
+***********************************************
+*Charge computation*
+***********************************************
+*does not work with ngspice
+*Echarge charge 0 value={SDT(I(Emem))}
+***********************************************
+* WINDOW FUNCTIONS
+* FOR NONLINEAR DRIFT MODELING *
+***********************************************
+*window function, according to Joglekar
+.func f(x,p)={1-(2*x-1)^(2*p)}
+*proposed window function
+;.func f(x,i,p)={1-(x-stp(-i))^(2*p)}
+.ENDS memristor
diff --git a/MIL_testing/memristor_waveforms.png b/MIL_testing/memristor_waveforms.png
new file mode 100644
index 0000000..54cd69e
--- /dev/null
+++ b/MIL_testing/memristor_waveforms.png
diff --git a/MIL_testing/reading_output.py b/MIL_testing/reading_output.py
new file mode 100644
index 0000000..3e37dfe
--- /dev/null
+++ b/MIL_testing/reading_output.py
@@ -0,0 +1,36 @@
+from spicelib import RawRead
+import matplotlib.pyplot as plt
+
+raw = RawRead("simulation.raw")
+print(raw.get_trace_names())
+
+time = raw.get_trace("time").get_wave()
+current = raw.get_trace("i(v.xmem.v_emem)").get_wave()
+voltage = raw.get_trace("v(in)").get_wave()
+state_variable = raw.get_trace("v(xmem.x)").get_wave()
+
+print(f"final current: {current[-1]}")
+
+print(f"state variable : {state_variable[-1]}")
+
+# Create a figure with two stacked subplots
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
+
+# 1. Plot Current vs Time
+ax1.plot(time, current, color='blue', linewidth=2, label='Memristor Current')
+ax1.set_ylabel('Current (A)')
+ax1.set_title('Memristor Simulation Results')
+ax1.grid(True, linestyle='--', alpha=0.7)
+ax1.legend()
+
+# 2. Plot Voltage vs Time
+ax2.plot(time, voltage, color='red', linewidth=2, label='Input Voltage v(in)')
+ax2.set_xlabel('Time (s)')
+ax2.set_ylabel('Voltage (V)')
+ax2.grid(True, linestyle='--', alpha=0.7)
+ax2.legend()
+
+# Adjust layout to prevent overlap and display the plot
+plt.tight_layout()
+plt.savefig("memristor_waveforms.png", dpi=300, bbox_inches='tight')
+print("Saved waveforms to memristor_waveforms.png")
diff --git a/MIL_testing/simulation.cir b/MIL_testing/simulation.cir
new file mode 100644
index 0000000..8b0701d
--- /dev/null
+++ b/MIL_testing/simulation.cir
@@ -0,0 +1,34 @@
+* simulation.cir
+* Testbench for HP Memristor
+.options savecurrents
+.save all
+
+.include "hp_memristor.cir"
+
+* Voltage source
+.param Vin = ${VOLTAGE_IN}
+.param Vlast = ${VOLTAGE_LAST}
+V1 in 0 PULSE(VOLTAGE_LAST VOLTAGE_IN 0 1u 1u 1m 3m)
+
+* Memristor
+XMEM in out memristor 
++ Ron=1k 
++ Roff=100k 
++ Rinit=80k 
++ D=10n 
++ uv=10f 
++ p=1
+
+* Load resistor
+Rload out 0 1k
+
+.tran 1u 100m uic
+
+.control
+
+run
+write simulation.raw
+quit
+
+.endc
+.end
diff --git a/MIL_testing/simulation.raw b/MIL_testing/simulation.raw
new file mode 100644
index 0000000..0c51710
--- /dev/null
+++ b/MIL_testing/simulation.raw
diff --git a/MIL_testing/spice_backend.py b/MIL_testing/spice_backend.py
index 674ae30..c6be8a1 100644
--- a/MIL_testing/spice_backend.py
+++ b/MIL_testing/spice_backend.py
@@ -1,5 +1,6 @@
 import subprocess
 import re
+from spicelib import RawRead
 
 class SpiceBackend:
     """This is what actually does all of the spice stuff"""
@@ -8,6 +9,8 @@ class SpiceBackend:
         self.voltage = 0
         self.current = 0
         self.time = 0
+        self.state = None
+        self.dt = 1e-9
 
     def simulate(self, voltage_signal: float, return_current: bool = False) -> float:
         """runs the simulation for one time step"""
@@ -43,7 +46,7 @@ class SpiceBackend:
         #need to find a way to pass voltage signal to ngspice here
         #also maybe a time series resolution??
         result = subprocess.run(
-                ["ngspice" "-b" "spicefile.cir"],
+                ["ngspice" "-b" "simulation.cir"],
                 capture_output=True,
                 text=True
         )