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

+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

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.

scripts/profstats

-#!env python
+#!/usr/bin/env python2
 #coding: utf-8
 
 import json

scripts/tsplot

-#!env python
+#!/usr/bin/env python2
 #coding: utf-8
 
 import argparse