(* **************************************************
Kuramoto reduction of a 2-dimensional system.
Can be used only after f[x,y,z], g[x,y,z],
yss[x], zss[x], and dcoup have been defined. 
************************************************** *)
   
reduction[f_,g_,yss_,zss_,dcoup_,xguess_] := (

vars = {x,y};
jac := Outer[D,{f[x,y,z],g[x,y,z]},vars];
m = Outer[D,jac,vars];
n = Outer[D,m,vars];

trjac = jac[[1,1]] + jac[[2,2]] /. {y->yss[x],
                                    z->zss[x]};
x0 = x /. FindRoot[ trjac==0, {x,xguess} ];
y0 = yss[x0];
z0 = zss[x0];

l0 = jac /. {x->x0,y->y0,z->z0};
m0 = m/2 /. {x->x0,y->y0,z->z0};
n0 = n/6 /. {x->x0,y->y0,z->z0};

rhs = { -D[f[x,y,z],z], -D[g[x,y,z],z] };
{dxdz,dydz} = LinearSolve[l0,rhs];
l1 = D[jac,x] dxdz + 
     D[jac,y] dydz + 
     D[jac,z] /. {x->x0,y->y0,z->z0};

esystem1 = Eigensystem[l0];
esystem2 = Eigensystem[-Transpose[l0]];
omega = Im[ esystem1[[1,1]] ];
u = esystem1[[2,1]];
ustar = esystem2[[2,2]];
ustar /= ustar.u;

vplus = -LinearSolve[(l0-2 omega I IdentityMatrix[2]), 
                     m0.u.u];
vminus = Conjugate[vplus];
vknot = -2 LinearSolve[l0,m0.u.Conjugate[u]];

ab = ustar.l1.u;
pq = -( 2 ustar.m0.u.vknot + 
        2 ustar.m0.Conjugate[u].vplus + 
        3 ustar.n0.u.u.Conjugate[u] );
d1d2 = ustar.dcoup.u;

alpha = Im[d1d2] / Re[d1d2];
beta = Im[ab] / Re[ab];
eta = Im[pq] / Re[pq];

eta1 = - (5 + alpha^2) / (4 alpha);
eta2 = -1 / alpha;
eta3 = (alpha^2 - 1) / (2 alpha);
If[ alpha < 0,
    If[ eta > eta1, region = 1, 
    If[ eta > eta2, region = 2, 
    If[ eta > eta3, region = 3, region = 4 ]]],
    If[ eta < eta1, region = 1,
    If[ eta < eta2, region = 2,
    If[ eta < eta3, region = 3, region = 4 ]]] ];
   
matrixout = { {x0,y0,z0}, 
              {Re[ab],Re[pq],Re[d1d2]}, 
              {Im[ab],Im[pq],Im[d1d2]},
              {alpha,beta,eta,region} } )

(* **************************************************
Reduction of the non-dimensional generic beta
cell model from the BMB paper by Artie. 
************************************************** *)
   
reductionBMB[lambda_,gca_,xguess_] := 
Module[{minf,ninf,f,g,yss,zss,dcoup},(

(* These are the equations for the first 
   nondimensionalization, resulting in the 
   upside-down bifurcation diagram.

minf[x_] := 1/(1+Exp[(25x-45)/12]);
ninf[x_] := 1/(1+Exp[(25x-42)/5.6]);

f[x_,y_,z_] := -gca minf[x] x/lambda -
               10 y (x-4)/lambda -
               z (x-4);
g[x_,y_,z_] := ninf[x] - y;

yss[x_] := ninf[x];
zss[x_] := ( -gca minf[x] x/lambda - 
             10 ninf[x] (x-4)/lambda ) / (x-4); *)
             
(* And these are the equations for the second
   nondimensionalization, resulting in a right-side-up
   bifurcation diagram. *)
   
minf[x_] := 1/(1+Exp[(205-25x)/12]);
ninf[x_] := 1/(1+Exp[(208-25x)/5.6]);

f[x_,y_,z_] := -gca minf[x] (x-10)/lambda -
               10 y (x-6)/lambda -
               z (x-6);
g[x_,y_,z_] := ninf[x] - y;

yss[x_] := ninf[x];
zss[x_] := ( -gca minf[x] (x-10)/lambda - 
             10 ninf[x] (x-6)/lambda ) / (x-6);
             
dcoup = {{1,0},{0,0}};

reduction[f,g,yss,zss,dcoup,xguess] )]

MatrixForm[ reductionBMB[1.1,3.6,7.8] ]

MatrixForm[ reductionBMB[0.9,3.6,7.9] ]

MatrixForm[ reductionBMB[1.1,3.6,2.2] ]

MatrixForm[ reductionBMB[0.7,3.6,2.2] ]

varylambda[lmin_,lmax_,lstep_,gca_,xguess_] := (

x0 = xguess;
temp = {};
Do[ reduced = reductionBMB[lambda,gca,x0];
    x0 = reduced[[1]][[1]];
    z0 = reduced[[1]][[3]];
    p = reduced[[2]][[2]];
    alpha = reduced[[4]][[1]];
    eta = reduced[[4]][[3]];
    region = reduced[[4]][[4]];
    temp = Append[temp,{lambda,gca,x0,z0,p,alpha,eta,region}],
    
    {lambda,lmin,lmax,lstep} ]; )

varygca[lambda_,gmin_,gmax_,gstep_,xguess_] := (

x0 = xguess;
temp = {};
Do[ reduced = reductionBMB[lambda,gca,x0];
    x0 = reduced[[1]][[1]];
    z0 = reduced[[1]][[3]];
    p = reduced[[2]][[2]];
    alpha = reduced[[4]][[1]];
    eta = reduced[[4]][[3]];
    region = reduced[[4]][[4]];
    temp = Append[temp,{lambda,gca,x0,z0,p,alpha,eta,region}],
    
    {gca,gmin,gmax,gstep} ]; )

PrintTable[table_,filename_] := (
out = OpenWrite[filename];
Do[ WriteString[out,ToString[Part[Part[table,i],1]], "  ",
                    ToString[Part[Part[table,i],2]], "  ",
                    ToString[Part[Part[table,i],3]], "  ",
                    ToString[Part[Part[table,i],4]], "  ",
                    ToString[Part[Part[table,i],5]], "  ",
                    ToString[Part[Part[table,i],6]], "  ",
                    ToString[Part[Part[table,i],7]], "  ",
                    ToString[Part[Part[table,i],8]], "\n"],
    {i,1,Length[table]} ];
Close[out]; )

PrintTable[temp,"lambda.out"]

PrintTable[temp,"gca.out"]

varylambda[1.1,0.05,-0.05,3.6,2.25]

varygca[1.1,3.6,5.5,0.1,2.25]

MatrixForm[temp]

pars[delta_,gamma_,z1_] := (

gcoup = gamma Re[ab] (z1-z0) / Re[d1d2];
z2 = z0 + delta (z1-z0);

{gcoup,z2}

)

pars[0.98,0.5,0.56]

(* **************************************************
Reduction of the Morris-Lecar equations.
Using the nondimensionalized version of Han,
Kurrer, and Kuramoto.
************************************************** *)
   
reductionML[xguess_] := 
Module[{minf,yinf,tauy,f,g,yss,zss,dcoup},(

minf[x_] := (1+Tanh[(x+0.01)/0.15]) / 2;
yinf[x_] := (1+Tanh[(x-0.1)/0.145]) / 2;
tauy[x_] := 1 / Cosh[(x-0.1)/(2 0.145)];

f[x_,y_,z_] := z - 
               1 minf[x] (x-1) - 
               2 y (x+0.7) -
               0.5 (x+0.5);
g[x_,y_,z_] := 1.15 (yinf[x] - y) / tauy[x];

yss[x_] := yinf[x];

zss[x_] := 1 minf[x] (x-1) +
           2 yinf[x] (x+0.7) +
           0.5 (x+0.5);

dcoup = {{1,0},{0,0}};

reduction[f,g,yss,zss,dcoup,xguess] )]

MatrixForm[ reductionML[0] ]

(* **************************************************
Reduction of the Hindmarsh-Rose model.
************************************************** *)

reductionHR[a_,b_,c_,d_,input_,lambda_,xguess_] := 
Module[{f,g,yss,zss,dcoup},(

f[x_,y_,z_] := y - a x^3 + b x^2 + input - z;
g[x_,y_,z_] := lambda (c - d x^2 - y);

yss[x_] = c - d x^2;
zss[x_] = c - d x^2 - a x^3 + b x^2 + input;

dcoup = {{1,0},{0,0}};

reduction[f,g,yss,zss,dcoup,xguess] )]

MatrixForm[ reductionHR[1,3,1,5,2,1,5] ]

MatrixForm[ reductionHR[1,3,1,5,2,1,0] ]

MatrixForm[ reductionHR[1,3,1,5,2,2,5] ]

MatrixForm[ reductionHR[1,3,1,5,2,2,0] ]

MatrixForm[ reductionHR[1,3,1,5,2,2.5,5] ]

MatrixForm[ reductionHR[1,3,1,5,2,2.5,0] ]

MatrixForm[ reductionHR[1.1,2.8,1,4,2,2,3] ]

MatrixForm[ reductionHR[1,2.5,1,5,2,0.8,4] ]

(* **************************************************
Reduction of the model by Kawato, Sokabe and Suziki,
from 
Biol. Cybernetics, V34, 81-89, 1979.
************************************************** *)

reductionKSS[apar_,bpar_,cpar_,xguess_] := 
Module[{f,g,yss,zss,dcoup},(

f[x_,y_,z_] := cpar ( y + x - x^3/3 + z );
g[x_,y_,z_] := -( x - apar + bpar y ) / cpar;

yss[x_] := ( apar - x ) / bpar;
zss[x_] := x^3/3 - x - ( apar - x ) / bpar;

dcoup = {{1,0},{0,0}};

reduction[f,g,yss,zss,dcoup,xguess] )]

MatrixForm[ reductionKSS[0.7, 0.8, 1, 1 ] ]

(* **************************************************
Reduction of the Pernarowski model.  Passing the
coupling matrix as a parameter.
************************************************** *)

reductionMP[ua_,uhat_,ueta_,d_,xguess_] := 
Module[{f,g,yss,zss,dcoup},(

f[x_,y_,z_] := y;
g[x_,y_,z_] := -( ua y ((x-uhat)^2-ueta^2) +
                  z + x^3 - 3 (x+1) );

yss[x_] := 0;
zss[x_] := 3 (x+1) - x^3;

dcoup = d;

reduction[f,g,yss,zss,dcoup,xguess] )]

MatrixForm[ reductionMP[0.25,1.6,0.7,{{1,0},{1,1}},2] ]

(* **************************************************
Reduction of the Taylor-Kevrekidis model, from 
Physica D 51 (1991) 274-292.
Using B as the bifurcation parameter (z).
************************************************** *)

reductionTK[capd_,parbeta_,ysubc_,xguess_] := 
Module[{f,g,yss,zss,dcoup},(

f[x_,y_,z_] := -x + capd (1-x) Exp[y];
g[x_,y_,z_] := -y + z capd (1-x) Exp[y] - 
               parbeta (y-ysubc);

yss[x_] := Log[x/(capd (1-x))];
zss[x_] := (yss[x] (1+parbeta) - parbeta ysubc)/x;

dcoup = {{1,0},{0,1}};

reduction[f,g,yss,zss,dcoup,xguess] )]

MatrixForm[ reductionTK[0.15541882,1.05,0,0.5] ]

MatrixForm[ reductionTK[0.15541882,1.05,0,0.8] ]