#! /usr/bin/gawk -f
# Last edited on 2019-11-06 17:40:46 by jstolfi

function F(z,kappa,lambda,dead,per,  p,Fz) {
  # Return the {per} iterate of {F} at {z}.
  Fz = z;
  for(p = 1; p <= per; p++) { 
    # At this point, Fz is the {p-1} iterate of {F} at {z}:
    Fz = F1(Fz,kappa,lambda,dead)
  }
  return Fz;
}

function DF(z,kappa,lambda,dead,per, p,Fz,DFz) {
  # Return the derivative of the {per} iterate of {F} at {z}.
  Fz = z;
  DFz = 1;
  for(p = 1; p <= per; p++) { 
    # At this point, {Fz} and {DFz} are the value and derivative of
    # the {p-1} iterate of {F} at {z}:
    DFz = DFz*DF1z(Fz,kappa,lambda,dead)
    Fz = (1-Fz)*phi(kappa*Fz+lambda);
  }
  return DFz;
}

function F1(z,kappa,lambda,dead) {
  # Returns the value of {F} at {z}.
  if (dead == 0) {
    return phi(kappa*z+lambda)
  } else {
    return (1-z)*phi(kappa*z+lambda)
  }
}

function DF1z(z,kappa,lambda,dead) {
  # Returns the derivative of {F} at {z}.
  if (dead == 0){
    return Dphi(kappa*z+lambda)
  } else {
    return -phi(kappa*z + lambda) + (1 - z)*Dphi(kappa*z+lambda)*kappa
  }
}

function phi(u) {
  # Piecewise-linear firing function with normalized threshhold 
  # potential {0} and normalized saturation potential {1},
  # evaluated at {u}.
  if (u <= 0.0) { return 0.0; }
  else if (u >= 1.0){ return 1.0; }
  else { return u; }
}

function Dphi(u) {
  # Derivative of {phi} at {u}.
  if (u <= 0.0) { return 0.0; }
  else if (u >= 1.0){ return 0.0; }
  else { return 1.0; }
}

function mknum(x) {
  # Converts {x} to number, stripping any quotes.
  gsub(/[\"\']/, "", x);
  return x + 0;
}

function min(x,y) {
  return (x <= y ? x : y);
}

function max(x,y) {
  return (x >= y ? x : y);
}

function abs(x) {
  return (x >= 0 ? x : -x);
}

function error(msg) {
  printf "** %s\n", msg > "/dev/stderr";
  exit(1);
}