-- Electrophysiology and common stuff for nerve models --

module(..., package.seeall);

-- useful constants
R = 8314.4        -- Universal gaz constant, mJ / mole*K
F = 96485.34      -- Faraday constant Coulomb / mole

-- Diffusion coefficients in aqueous solutions --
-- [cm^2/s]
-- (Hille, B., Ionic Channels of Excitable Membranes)
Daq = {K = 1.96e-5, Na = 1.33e-5, H = 9.31e-5, Ca = 0.79e-5}



local gauss_mem = {}
setmetatable(gauss_mem, {__mode = "v"})  -- make gauss_mem weak

function gaussian(t,p)
   --local key = t..'-'..p.t..'-'..p.w..'-'..p.a
   --if gauss_mem[key] then  
   --   return gauss_mem[key]
   --else
      local x,g,z
      z = 10
      if (t - p.t) * (t - p.t) < z * p.w then
	 g = (t - p.t)^2
	 g = -.5*g / p.w
	 x =  p.a * math.exp(g)
      else x = 0 end
      gauss_mem[key] = x
      return x
   --end
end

-- Some biophysics --

function NaKpump_flux(Ki, Ko, Nai, Nao, p)
   --                                mmole
   -- NOTE: should return fluxes in --------
   --                                cm^2 s  
   -- version 1.0 -- no optimization
   
   local jna, jk
   
   local Km_na = p.a_na + p.b_na * Ki
   local Km_k = p.a_k + p.b_k * Nao
   
   local mm_na = 1 / (1 + Km_na/Nai)^3 
   local mm_k =  1 / (1 + Km_k/Ko)^2

   jna = p.jnamax*mm_na*mm_k --+ arg.p.jrest
   jk = -jna / 1.5
   
   return {jna, jk}
end


-- Basic circuits --

function i2flux(I)
   return I/(1e3*F)
end

function flux2i(flux)
   return flux*1e3*F
end

function reciprocal(x)  return 1/x end

function rsumr(values)
   -- reciprocal of the summed reciprocals --
   local res = 0
   for i,g in ipairs(values) do
   --for i = 1,#values do 
      res = res + 1/g
   end
   return 1/res
end

function sum(values)
   local res = 0
   for _,v in ipairs(values) do res = res+v end
   return res
end

function series_resistances(rrr)
   return sum(rrr)
end

function series_conductances(ggg)
   return rsumr(ggg)
end

function parallel_resistances(rrr)
   return rsumr(rrr)
end

function parallel_conductances(ggg)
   return sum(ggg)
end

-- Electrochemistry functions --

function nernst(ci, co, T, z)
   if not z then z = 1 end -- default charge
   return R*T*math.log(co / ci)/(z*F)
end


-- GHK current equation --
function ghk(ci, co, V, T, z) 
   if not z then z = 1 end -- default charge
   local x = V*F*z / (R*T)
   if math.abs(V) < 1e-7 then
      return (ci - co)*R*T/z
   else
      local y = math.exp(x)
      return x*F*(co-ci*y)/(1 - y)
   end
end


local ghk_mem = {}
setmetatable(ghk_mem, {__mode = "v"})  -- make ghk_mem weak


function make_subkey(vector)
   local key = "-"
   for i,val in ipairs(vector) do
      key = key..'val'..'-'
   end
   return key
end

function ghk_rate_const(p,V)
      local x = V + p[1]
      local rc = (p[2]*x + p[3]) / (1 + p[4]*math.exp(x/p[5]))
      return rc 
end

function myelin_area(D, L)
    --return math.pi*(D+d)*(L + .5*(D-d))
    --return math.pi*(D*L + .5*(D^2 - d^2))
    return math.pi*D*L
end

function ring_area(D, d)
   return .25*math.pi*(D^2 - d^2)
end

function periaxonal_conductance(rho, d, h, l)
   return ring_area(d + 2*h, d)/ (rho*l)
end

function myelin_conductance(gs, D, d, L, N)
   -- two membranes per lamella
   --return .5 * math.pi * gs * L * D / N
   return myelin_area(D, L)*gs/(2*N)
end

function myelin_capacitance(cs, D, d,L, N)
   -- two membranes per lamella
   --return .5 * math.pi * cs * L * D / N 
   return cs *myelin_area(D, L) /(2*N) 
end