From 1cbc812e59ce51c24adc9e4c0798ff350e222045 Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Wed, 19 Oct 2016 16:46:31 +0200
Subject: [PATCH] Add analytic simple cable solver.

* New file (Julia module): PassiveCable.jl
* Update script documentation for profstats
* Add (too brief) documentation for PassiveCable.jl
* Fix `env`, `python` invocation on python scripts.
---
 scripts/PassiveCable.jl | 104 ++++++++++++++++++++++++++++++++++++++++
 scripts/README.md       |  56 ++++++++++++++++++++++
 scripts/profstats       |   2 +-
 scripts/tsplot          |   2 +-
 4 files changed, 162 insertions(+), 2 deletions(-)
 create mode 100644 scripts/PassiveCable.jl

diff --git a/scripts/PassiveCable.jl b/scripts/PassiveCable.jl
new file mode 100644
index 00000000..7c50a645
--- /dev/null
+++ b/scripts/PassiveCable.jl
@@ -0,0 +1,104 @@
+module PassiveCable
+
+export cable_normalize, cable, rallpack1
+
+# Compute solution g(x, t) to
+#
+#     ∂²g/∂x² - g - ∂g/∂t = 0
+#
+# on [0, L] × [0,∞), subject to:
+#
+#     g(x, 0) = 0
+#     ∂g/∂x (0, t) = 1
+#     ∂g/∂x (L, t) = 0
+#
+# (This implementation converges slowly for small t)
+
+function cable_normalized(x, t, L; tol=1e-8)
+    if t<=0
+        return 0.0
+    else
+        ginf = -cosh(L-x)/sinh(L)
+        sum = exp(-t/L)
+
+        for k = countfrom(1)
+            a = k*pi/L
+            e = exp(-t*(1+a^2))
+
+            sum += 2/L*e*cos(a*x)/(1+a^2)
+            resid_ub = e/(L*a^3*t)
+
+            if resid_ub<tol
+                break
+            end
+        end
+	return ginf+sum
+     end
+end
+
+
+# Compute solution f(x, t) to
+#
+#     λ²∂²f/∂x² - f - τ∂f/∂t = 0
+#
+# on [0, L] x [0, ∞), subject to:
+#
+#     f(x, 0) = V
+#     ∂f/∂x (0, t) = I·r
+#     ∂f/∂x (L, t) = 0
+#
+# where:
+#
+#     λ² = 1/(r·g)   length constant
+#     τ  = r·c       time constant
+#
+# In the physical model, the parameters correspond to the following:
+#
+#     L:  length of cable
+#     r:  linear axial resistivity
+#     g:  linear membrane conductance
+#     c:  linear membrane capacitance.
+#     V:  membrane reversal potential
+#     I:  injected current on the left end (x = 0) of the cable.
+
+function cable(x, t, L, lambda, tau, r, V, I; tol=1e-8)
+    scale = I*r*lambda;
+    if scale == 0
+	return V
+    else
+        tol_n = abs(tol/scale)
+	return scale*cable_normalized(x/lambda, t/tau, L/lambda, tol=tol_n) + V
+    end
+end
+
+
+# Rallpack 1 test
+#
+# One sided cable model with the following parameters:
+#
+#     RA = 1 Ω·m    bulk axial resistivity
+#     RM = 4 Ω·m²   areal membrane resistivity
+#     CM = 0.01 F/m²  areal membrane capacitance
+#     d  = 1 µm     cable diameter
+#     EM = -65 mV   reversal potential
+#     I  = 0.1 nA   injected current
+#     L  = 1 mm     cable length 
+#
+
+function rallpack1(x, t; tol=1e-8)
+    RA = 1
+    RM = 4
+    CM = 1e-2
+    d  = 1e-6
+    EM = -65e-3
+    I  = 0.1e-9
+    L  = 1e-3
+
+    r = 4*RA/(pi*d*d)
+    lambda = sqrt(d/4 * RM/RA)
+    tau = CM*RM
+
+    return cable(x, t, L, lambda, tau, r, EM, -I, tol=tol)
+end
+
+end #module
diff --git a/scripts/README.md b/scripts/README.md
index e645a762..230d5be4 100644
--- a/scripts/README.md
+++ b/scripts/README.md
@@ -78,4 +78,60 @@ of the timeseries data.
 Use the `-o` or `--output` option to save the plot as an image, instead of
 displaying it interactively.
 
+profstats
+---------
 
+`profstats` collects the profiling data output from multiple MPI ranks and performs
+a simple statistical summary.
+
+Input files are in the JSON format emitted by the profiling code.
+
+By default, `profstats` reports the quartiles of the times reported for each
+profiling region and subregion. With the `-r` option, the collated raw times
+are reported instead.
+
+Output is in CSV format.
+
+PassiveCable.jl
+---------------
+
+Compute analytic solutions to the simple passive cylindrical dendrite cable
+model with step current injection at one end from t=0.
+
+This is used to generate validation data for the first Rallpack test.
+
+Module exports the following functions:
+
+ * `cable_normalized(x, t, L; tol)`
+
+   Compute potential _V_ at position _x_ in [0, _L_] at time _t_ ≥ 0 according
+   to the normalized cable equation with unit length constant and time constant.
+
+   Neumann boundary conditions: _V'_(0) = 1; _V'_(L) = 0.
+   Initial conditions: _V_( _x_, 0) = 0.
+
+   Absolute tolerance `tol` defaults to 1e-8.
+
+ * `cable(x, t, L, lambda, tau, r, V, I; tol)`
+
+   Compute the potential given:
+      *  length constant `lambda`
+      *  time constant `tau`,
+      *  axial linear resistivity `r`
+      *  injected current of `I` at the origin
+      *  reversal potential `V`
+
+   Implied units must be compatible, e.g. SI units.
+
+   Absolute tolerance `tol` defaults to 1e-8.
+
+ * `rallpack1(x, t; tol)`
+
+   Compute the value of the potential in the Rallpack 1 test model at position
+   `x` and time `t`.
+
+   Parameters for the underlying cable equation calculation are taken from the
+   Rallpack model description in SI units; as the cable length is 1 mm in this
+   model, `x` can take values in [0, 0.001].
+
+   Absolute tolerance `tol` defaults to 1e-8.
diff --git a/scripts/profstats b/scripts/profstats
index 33f37356..88f68c72 100755
--- a/scripts/profstats
+++ b/scripts/profstats
@@ -1,4 +1,4 @@
-#!env python
+#!/usr/bin/env python2
 #coding: utf-8
 
 import json
diff --git a/scripts/tsplot b/scripts/tsplot
index 9ad6e2a5..8e89cbda 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -1,4 +1,4 @@
-#!env python
+#!/usr/bin/env python2
 #coding: utf-8
 
 import argparse
-- 
GitLab