######################################################################
##ROKHLENKO.txt: Save this file as  ROKHLENKO.txt                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read ROKHLENKO.txt<Enter>                          #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################

#accompanying Fowler-Nordheim emission produced by laser radiation
#A.Rokhlenko and J.L.Lebowitz
 
#Created: Dec. 3, 2017
 
print(`Created: Dec. 3, 2017`):
print(` This is  ROKHLENKO.txt `):
print(`It is a short Maple package that is a computational appendix to the article `):
print(`Fowler-Nordheim emission produced by laser radiation`):
print(`by A.Rokhlenko and J.L.Lebowitz`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Sequ `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: S, Story `):
 print(` `):



elif nops([args])=1 and op(1,[args])=S then
print(`S(N,Di): The Rokhlenko equations (10)-(12) in terms of symbolic h,V, using Digits Di. `):
print(`It inputs a positive integer and solves, numerically, the equations  (10)-12 in the article`):
print(` "Fowler-Nordheim emission produced by laser radiation".  `):
print(`by A.Rokhlenko and J.L.Lebowitz`):
print(`It outputs the vector [c[0], ..., c[N]]`):

elif nops([args])=1 and op(1,[args])=Story then
print(`Story(Ntop,Di): tells the story of [c[0], ..., c[N]], satisfying the Rokhlenko equations (10)-(12) as given in the article`):
print(` "Fowler-Nordheim emission produced by laser radiation".  `):
print(`by A.Rokhlenko and J.L.Lebowitz`):
print(`For N from 1 to Ntop`):
print(`using Di Digits. Try:`):
print(`Story(7,20);`):

elif nops([args])=1 and op(1,[args])=Sequ then
 print(`Sequ(N,a,c,qnM,qnP): The Rokhlenko equations (10)-12 in terms of symbolic h,V0,hw,m,eF,r. N`):
 print(`Only the equations, w/o solving it`):

else
print(`There is no ezra for`,args):
fi:
 
end:



delt:=proc(n)
if n=0 then
  RETURN(1):
 else
  RETURN(0):
fi:
end:

#Sequ(N,a,c,qnM,qnP): The Rokhlenko equations (10)-12 in terms of symbolic h,V0,hw,m,eF,r. N is a positive integer
Sequ:=proc(N, a,c,qnM,qnP) local eq,var,gamma,k,n,eq1,eq2,eq3,eq4,i1,W,eF,r:
var:={seq(a[n],n=0..N), seq(c[n],n=0..N), seq(qnP[n],n=-1..N+1), seq(qnM[n],n=-1..N+1)}:

#1. W=V/2= 6.88*10^{-12},
#2. W_{n}=(6.88+2.48*n)*10^{-12},
#3. W_{-n}=(6.88-2.48*n)*10^{-12},
#4. k_n=4.04*sqrt{W_n}*10^{13},
#5. gamma_n= 4.04*sqrt{W_{-n}}*10^{13},
#6. eF=6.4*10^{-5},
#7. r=1.5*10^{-7}.


for n from -N to N do
 W[n]:=(6.88+2.48*n)*10^(-12):
od:

for n from -N to N do
k[n]:=4.04*sqrt(W[n])*10^(13):
gamma[n]:=4.04*sqrt(W[-n])*10^(13):
od:

eF:=6.4*10^(-5):
r:=1.5*10^(-7):


#(10a)
eq1:=[seq(delt(n,0)+a[n]=qnP[n]+qnM[n],n=0..N)]:

#(10b)
eq2:= [seq(I*k[n]*(delt(n)-a[n])=gamma[n]*(qnP[n]-qnM[n])-eF/2/I/W[n]*(qnP[n+1]-qnM[n+1]-qnP[n-1]+qnM[n-1]),n=0..N)]:

#(11a) (typo corrected)
eq3:=[seq(
qnP[n]*exp(gamma[n]*r)+qnM[n]*exp(-gamma[n]*r)
-eF*r/2/I/W[n]*
     (qnP[n+1]*exp(gamma[n+1]*r)+qnM[n+1]*exp(-gamma[n+1]*r)-
      qnP[n-1]*exp(gamma[n-1]*r)-qnM[n-1]*exp(-gamma[n-1]*r))=
     c[n]*exp(I*k[n]*r)+c[n-1]*exp(I*k[n-1]*r)+ c[n+1]*exp(I*k[n+1]*r),n=0..N)]:

#(11b)
eq4:=[seq(
gamma[n]*(qnP[n]*exp(gamma[n]*r) - qnM[n]*exp(-gamma[n]*r))
-eF/2/I/W[n]*
     ( qnP[n+1]*exp(gamma[n+1]*r)+qnM[n+1]*exp(-gamma[n+1]*r)-
      qnP[n-1]*exp(gamma[n-1]*r)-qnM[n-1]*exp(-gamma[n-1]*r))=
     c[n]*exp(I*k[n]*r)+c[n-1]*exp(I*k[n-1]*r)+ c[n+1]*exp(I*k[n+1]*r),
n=0..N)]:



eq:=[op(eq1),op(eq2),op(eq3),op(eq4)]:
eq:=subs({qnP[-1]=0,qnM[-1]=0,qnP[N+1]=0,qnM[N+1]=0,c[-1]=0,c[N+1]=0},eq):
var:=[seq(a[i1],i1=0..N), seq(qnM[i1],i1=0..N), seq(qnP[i1],i1=0..N),seq(c[i1],i1=0..N)]:
eq,var:

end:




#S(N,Di): The Rokhlenko equations (10)-(12) in terms of symbolic h,V, using Digits Di. 
#It inputs a positive integer and solves, numerically, the equations  (10)-12 in the article
#"Fowler-Nordheim emission produced by laser radiation". 
#by A.Rokhlenko and J.L.Lebowitz`):
#It outputs the vector [c[0], ..., c[N]]

S:=proc(N,Di) local eq,var,gamma,k,n,eq1,eq2,eq3,eq4,W,eF,r, a,c,qnM,qnP,var1,i:
var:={seq(a[n],n=0..N), seq(c[n],n=0..N), seq(qnP[n],n=-1..N+1), seq(qnM[n],n=-1..N+1)}:
Digits:=Di:
#1. W=V/2= 6.88*10^{-12},
#2. W_{n}=(6.88+2.48*n)*10^{-12},
#3. W_{-n}=(6.88-2.48*n)*10^{-12},
#4. k_n=4.04*sqrt{W_n}*10^{13},
#5. gamma_n= 4.04*sqrt{W_{-n}}*10^{13},
#6. eF=6.4*10^{-5},
#7. r=1.5*10^{-7}.


for n from -N to N do
 W[n]:=(6.88+2.48*n)*10^(-12):
od:

for n from -N to N do
k[n]:=4.04*sqrt(W[n])*10^(13):
gamma[n]:=4.04*sqrt(W[-n])*10^(13):
od:

eF:=6.4*10^(-5):
r:=1.5*10^(-7):


#(10a)
eq1:=[seq(delt(n,0)+a[n]=qnP[n]+qnM[n],n=0..N)]:

#(10b)
eq2:= [seq(I*k[n]*(delt(n)-a[n])=gamma[n]*(qnP[n]-qnM[n])-eF/2/I/W[n]*(qnP[n+1]-qnM[n+1]-qnP[n-1]+qnM[n-1]),n=0..N)]:

#(11a) (typo corrected)
eq3:=[seq(
qnP[n]*exp(gamma[n]*r)+qnM[n]*exp(-gamma[n]*r)
-eF*r/2/I/W[n]*
     (qnP[n+1]*exp(gamma[n+1]*r)+qnM[n+1]*exp(-gamma[n+1]*r)-
      qnP[n-1]*exp(gamma[n-1]*r)-qnM[n-1]*exp(-gamma[n-1]*r))=
     c[n]*exp(I*k[n]*r)+c[n-1]*exp(I*k[n-1]*r)+ c[n+1]*exp(I*k[n+1]*r),n=0..N)]:

#(11b)
eq4:=[seq(
gamma[n]*(qnP[n]*exp(gamma[n]*r) - qnM[n]*exp(-gamma[n]*r))
-eF/2/I/W[n]*
     ( qnP[n+1]*exp(gamma[n+1]*r)+qnM[n+1]*exp(-gamma[n+1]*r)-
      qnP[n-1]*exp(gamma[n-1]*r)-qnM[n-1]*exp(-gamma[n-1]*r))=
     c[n]*exp(I*k[n]*r)+c[n-1]*exp(I*k[n-1]*r)+ c[n+1]*exp(I*k[n+1]*r),
n=0..N)]:



eq:=[op(eq1),op(eq2),op(eq3),op(eq4)]:
eq:=subs({qnP[-1]=0,qnM[-1]=0,qnP[N+1]=0,qnM[N+1]=0,c[-1]=0,c[N+1]=0},eq):

var1:=solve(eq,var):

if var1=NULL then
 RETURN(FAIL):
fi:

evalf(subs(var1,[seq(c[i],i=0..N)]),3):



end:


#Story(Ntop,Di): tells the story of [c[0], ..., c[N]], satisfying the Rokhlenko equations for N from 1 to Ntop
#using Di Digits. Try:
#Story(7,20);
Story:=proc(Ntop,Di)  local N,c,i,gu,t0:
t0:=time():
print(`These are the solutions of the equations (10)-(12) in the article `):
print(` "Fowler-Nordheim emission produced by laser radiation".  `):
print(`by A.Rokhlenko and J.L.Lebowitz`):

print(`only giving [c[0], ..., c[N]] for all N for N from 1 to `, Ntop ):
print(``):

for N from 1 to Ntop do
 gu:=S(N,Di):
  if gu<>FAIL then
 print(`When, N=`, N, `we have`):
 print(``):
 print([seq(c[i],i=0..N)]=gu):
 print(``):
 fi:
od:

print(`This ends this file, that took`, time()-t0, `seconds. `):

end: