# Worksheet for the Large Amplitude Impedance for mechanisms with a
# single adsorbed inermediate.
#  Remember that I no longer means sqrt(-1)!
restart;with(linalg):with(plots):alias(I=BesselI,j=sqrt(-1)):
# Choose the order of harmonics we want to calculate to.  kmax is the
# maximum harmonic to be solved.
# e,g, kmax = 2 means up to second harmonic. mx is then the size of
# matrix
kmax:=2;mx:=2*kmax+1;
# Define a function for the (finite) series expansion for
# exp(b*sin(omega*t))
 actafel:=b->         
 I(0,b)+2*sum((-1)^'k'*I(2*'k'+1,b)*sin((2*'k'+1)*omega*t),'k'=0..kmax)
       +2*sum((-1)^'k'*I(2*'k',b)  *cos((2*'k') 
 *omega*t),'k'=1..kmax):         actafel(b):
# Fourier series for theta
 Theta:=theta[0]+sum(theta['k',`r`]*sin('k'*omega*t)+theta['k',`i`]*cos
 ('k'*omega*t),'k'=1..mx):
# Now get dTheta/dt
 dThetadt:=diff(Theta,t):
# The various harmonics of the derivative with respect to time are
 for k to mx do
  LHScos.k:=select(has,dThetadt,cos(k*omega*t))/cos(k*omega*t);
  LHSsin.k:=select(has,dThetadt,sin(k*omega*t))/sin(k*omega*t);
 od:
# DEFINE THE MECHANISM-SPECIFIC EQUATIONS HERE
# Define the rate law dtheta/dt = RHS, and the current function, each in
# the form of a Fourier series. Maple's combine
# function allows expression of products of Fourier series in a single
# Fourier series.
# Simplest to work in rate constants of units s^(-1) = rate constants in
# mol m^(-2) s^(-1) /Gamma in mol m^(-2).
# Then the current is also in units of per second.
# If the potential dependence is of interest, put it in the rate
# constants here.
# Example 1: Simple oxidative quasi-reversible adsorption, Langmuir
# isotherm, beta =1/2. V is a non-dim voltage
 v1:=k0*exp(V)*(1-Theta)*actafel(b)-k0*exp(-V)*Theta*actafel(-b):
 RHS:=combine(v1):
 current:=RHS:
# Example 2: The hydrogen evolution reaction, Langmuir isotherm, both
# betas=1/2. Not interested in potential dependence
 v1:=k1*(1-Theta)*actafel(-b)-km1*Theta*actafel(b):
 v2:=k2*Theta*actafel(-b)-km2*(1-Theta)*actafel(b):
 v3:=k3*Theta^2-km3*(1-Theta)^2:
 RHS:=combine(v1-v2-2*v3):
 current:=combine(-v1-v2):
# Collect terms so that we have a sum of ( harmonic * series of Bessel
# functions)
 for k to mx do
  RHS:=collect(RHS,cos(k*omega*t));
  RHS:=collect(RHS,sin(k*omega*t));
  current:=collect(current,cos(k*omega*t));
  current:=collect(current,sin(k*omega*t));
  od:RHS:current:
# Get rid of all the harmonics above mx 
 for k from mx+1 to 2*mx do
    RHS:=remove(has,RHS,cos(k*omega*t));
    RHS:=remove(has,RHS,sin(k*omega*t));
    current:=remove(has,current,cos(k*omega*t));
    current:=remove(has,current,sin(k*omega*t));
 od:RHS:current:
# Extract the coefficients of the harmonics we want. After these are
# removed the remaining parts of RHS and current are the dc parts
 for k to mx do
  RHScos.k:=select(has,RHS,cos(k*omega*t))/cos(k*omega*t);
  RHS:=remove(has,RHS,cos(k*omega*t));
  RHSsin.k:=select(has,RHS,sin(k*omega*t))/sin(k*omega*t);
  RHS:=remove(has,RHS,sin(k*omega*t));
  currentcos.k:=select(has,current,cos(k*omega*t))/cos(k*omega*t);
  current:=remove(has,current,cos(k*omega*t));
  currentsin.k:=select(has,current,sin(k*omega*t))/sin(k*omega*t);
  current:=remove(has,current,sin(k*omega*t));
 od:RHScos0:=RHS:currentcos0:=current:
# The equations to solve are (in real and imaginary theta variables) (sc
# means sin cos here)
 eqnssc:=NULL:varssc:=NULL:
 for k to mx do
   varssc:=varssc,theta[k,`r`],theta[k,`i`];
   eqns.k:=-RHSsin.k=-LHSsin.k;
   eqnc.k:=-RHScos.k=-LHScos.k;
  # print(eqns.k);print(eqnc.k);
   eqnssc:=eqnssc,eqns.k,eqnc.k:
 od:
 eqnssc:=[-RHScos0=0,eqnssc]:
 varssc:=[theta[0],varssc]:
# Now define some procedures to actually solve the equations
# These procedures take the equations or the matrix Asc, the list of
# variables, and a list of assigments for the parameters.
# The result is a list of values for the theta components. For the
# linear methods components higher than kmax 
# are set to zero. Uses the global kmax which must be consistent with
# the number of entries in varssc & eqnssc
# For examples of how to use these see below.
# 1. Numerical Solution of the matrix equations.
 flin:=proc(A,vars,params) local mx,ii,jj,B,bb,soln,ans,extras;
 global kmax;
 mx:=2*kmax+1;
 B:=submatrix(A,1..mx,1..mx):
 for ii to linalg[rowdim](B) do for jj to linalg[coldim](B) do
    B[ii,jj]:=subs(op(params),B[ii,jj])
 od; od;
 B:=map(evalf,B);
 bb:=convert(submatrix(A,1..mx,[2*mx+2]),vector):
 for ii to linalg[vectdim](bb) do
    bb[ii]:=subs(op(params),bb[ii]);
 od;
 soln:=convert(linsolve(B,bb),list):
 [op(zip((x,y)->x=y,varssc,soln)),
     op(map(x->x=0,[op(mx+1..-1,varssc)]))]:
 end:
# 2. Numerical Solution of the nonlinear equations
 fnonlin:=proc(eqns,vars,params) local intervals,k,feqns;
 intervals:={seq(varssc[k]=-1..1,k=1..nops(varssc))}: 
 feqns:={seq(subs(op(params),eqnssc[k])
         ,k=1..nops(eqnssc))}:
 fsolve(feqns,{op(varssc)},intervals);
 end:
# 3. Symbolic Solution of the Matrix equations (only useful for 5th or
# lower order matrices)
 lin:=proc(A,vars,params) local ii,jj,B,bb,mx,soln;
 global kmax;
 mx:=kmax*2+1;
 B:=submatrix(A,1..mx,1..mx):
 for ii to linalg[rowdim](B) do for jj to linalg[coldim](B) do
    B[ii,jj]:=subs(op(params),B[ii,jj])
 od; od;
 bb:=convert(submatrix(A,1..mx,[2*mx+2]),vector):
 for ii to linalg[vectdim](bb) do
    bb[ii]:=subs(op(params),bb[ii]);
 od;
 soln:=convert(linsolve(B,bb),list):
 [op(zip((x,y)->x=y,varssc,soln)),
     op(map(x->x=0,[op(mx+1..-1,varssc)]))]:
 end:
# Create the matrix of coefficients. Only required or works if the
# equations are linear
 Asc:=genmatrix(eqnssc, varssc,'flag'):
# EXAMPLE. Potential-dependence of impedance for the adsorption
# reaction. Must have hit the adsorption mechanism part above.
# Method is numerical solution of matrix equations.
# Case of b=1, omega=100, k0=2, V varies.
# First assemple all the fixed parameters into a set. Admittances are
# defined as current/b, so have units of s^(-1); similarly for
# impedances.
# YY will be used for the second harmonic current/b^2.
 rateconsts:={k0=2.};
   pars:={op(rateconsts),omega=100.,b=1.}:
   ans:=flin(Asc,varssc,pars):
   Yreal:=subs(op(pars),op(ans),currentsin1/b):
   Yim:=subs(op(pars),op(ans),currentcos1/b):
   Y:=Yreal+j*Yim:
   Z:=1/Y:Zreal:=evalc(Re(Z)):Zim:=evalc(Im(Z)):
   YYreal:=subs(op(pars),op(ans),currentsin2/b^2):
   YYim:=subs(op(pars),op(ans),currentcos2/b^2):
   YY:=YYreal+j*YYim:
   ZZ:=1/YY:ZZreal:=evalc(Re(ZZ)):ZZim:=evalc(Im(ZZ)):
# e.g., plot real part of second harmonic admittance vs potential
 plot(YYreal,V=-6..6);
# EXAMPLE. Is the 2nd harmonic response  zero at the standard potential
# for the simple adsorption reaction.
#  Must have hit adsorption in the mechanism part above.
# Use the symbolic solution.
# Want V=0,k0=1, arbitrary omega and b.
# Note that the harmonics higher than kmax have been arbitrarily set to
# zero and are not found by solving
# any equations (this is for compatibility with the fnonlin routine).
 pars:={V=0,k0=1};
 ans:=lin(Asc,varssc,pars);
 ans2:=simplify(ans,exp);
# Simplify uses the recursion relationships of Bessel functions to leave
# only zero and first order functions.simplify(ans2);
# EXAMPLE. Amplitude-dependence of impedance for the HER with
# recombination. Must have hit the HER mechanism part above. Goes much
# more slowly.
# Method is non-linear numerical solution. Case of omega=10, ks as
# noted, b varies.
 rateconsts:={k1=1.,km1=1.,k2=1.,km2=2.,k3=2.,km3=1,omega=10.};
 for k from 1 to 10 do
   print(k);
   pars:={op(rateconsts),omega=10.,b=evalf(k/10)};
   ans:=fnonlin(eqnssc,varssc,pars);
   Yreal:=subs(op(pars),op(ans),currentsin1/b);
   Yim:=subs(op(pars),op(ans),currentcos1/b);
   Y:=Yreal+j*Yim:
   Z:=1/Y:Zreal:=evalc(Re(Z)):Zim:=evalc(Im(Z)):
   Yr[k]:=Yreal;Yi[k]:=Yim;Zr[k]:=Zreal;Zi[k]:=Zim;
   YYreal:=subs(op(pars),op(ans),currentsin2/b^2);
   YYim:=subs(op(pars),op(ans),currentcos2/b^2);
   YY:=YYreal+j*YYim:
   ZZ:=1/YY:ZZreal:=evalc(Re(ZZ)):ZZim:=evalc(Im(ZZ)):
   YYr[k]:=YYreal;YYi[k]:=YYim;ZZr[k]:=ZZreal;ZZi[k]:=ZZim; 
 od:
# Amplitude dependence of real part of admittance
 pointplot([seq([k/10,Yr[k]],k=1..10)]);