"****************************************"
==> copy and run this script in Eigenmath
"****************************************"
""
""
"*******************************************************************************"
"*                                 data types :                                *"
"*******************************************************************************"
" a = 450001                                  integer"
" b = 0.66666                                 real"
" c = 3 + 5*i                                 complex"
" d = 'bonjour'                               string"
do(a=450001, b=0.66666, c=3+5*i, d="bonjour")
a
b
c
d
" True = 1, False = 0                         boolean"
"                                             1==2"
1==2
"                                             2==2"
2==2
" v = (2, 6, 1)                               vector"
v=(2,6,1)
v
" A =  ((1, 5),(3,i))                         matrix"
A =((1,5),(3,2*i))
A
" pi, e, i                                    predefined"
pi
e
i
" complex number :"
" can be entered in rectangular or polar form :"
x = 1/2 + i/2
x
"                                            x=polar(x)"
x=polar(x)
x
"                                            x=rect(x)"
x=rect(x)
x
"                                            abs(x)"
abs(x)
""
"*******************************************************************************"
"*                         variables, expressions and functions :              *"
"*******************************************************************************"
"x"
x
" x = 126/3 + 17"
x = 126/3 + 17
x
" reinitialize x :  x = quote(x)"
x=quote(x)
x
" y = cos(theta) + i*sin(theta)"
y = cos(theta) + i*sin(theta)
y
" z = alpha+beta/3"
z = alpha+beta/3
z
" w = u^2 + v^2 + u/2 + v/2"
v=quote(v)
w = u^2 + v^2 + u/2 + v/2
w
" epsilon = 10^(-20)"
epsilon = 10^(-20)
epsilon
" Amunu = 3.0 "
Amunu=3.0
Amunu
" alphamunu = 3.0 "
alphamunu=3.0
alphamunu
" NA = 1"
NA=1
NA
" Na = 1"
Na=1
Na
" na = 1"
na=1
na

" S = 'this is a string' "
S = "this is a string"
S
" V = (1, 5, 6)    a vector"
V = (1, 5, 6)
V
" M = ((1, 4, 2),(5, 3, 9),(3, 8, 1))    a matrix"
M = ((1, 4, 2),(5, 3, 9),(3, 8, 1))
M
""
"f(x, y) = x + 3*y      a function"
f(x, y) = x + 3*y
"                                                   f(2,3)"
f(2,3)
""
"*******************************************************************************"
"*                     arithmetic and logical operations :                     *"
"*******************************************************************************"
" +   -   *   / "
" x^(y)"
" ! "
" <   <=    ==   >   >= "
" not(x)    and(x,y,...)    or(x,y,...) "
""
" examples :"
"                                          21^7"
21^7
"                                          float(21^7)"
float(21^7)
"                                          4^(1/2)"
4^(1/2)
"                                          10^(-3)"
10^(-3)
"                                          z^(-1/2)"
z^(-1/2)
"                                          100!"
17!
"                                          factor(last)"
factor(17!)
"                                          (a1 + i*a2)^2"
(a1 + i*a2)^2
"                                          and(3>=2, 5<15, 7==7)"
and(3>=2, 5<15)
"                                          and(x1>y1, x2>y2, not(y1==y2))"
""
"*******************************************************************************"
"*                                 control instructions :                      *"
"*******************************************************************************"
" do block = sequence of expressions :"
" do(expression1, expression2,...,expressionk)"
" do(x = sqrt(6), y = exp(x/2), z = x + y, float(z))"
do(x = sqrt(6), y = exp(x/2), z = x + y, float(z))
""
" test instructions :"
""
" if predicate is true execute block1"
"                 else execute block2"
" test (predicate, block1, block2)"
""
"      if predicate1 is true execute block1"
"           else if predicate is true execute block2, etc ..."
""
" test (predicate1,block1,predicate2,block2, ...)"
""
" examples :"
""
" f(a,b) = test(a>b, a 'is greater', b 'is greater')"
f(a,b) = test(a>b, a "is greater", b "is greater")
"                                                      f(3,4) "
f(3,4)

" f(a) = test(a<0, 0, a==0, 1, a>0, 0)"
f(a) = test(a<0, 0, a==0, 1, a>0, 0)
"                                                      f(0)"
f(0)
""
" range(u,a,b) = test(and(u>=a,u<=b),1,and(u>b,u<a),0)"
range(u,a,b) = test(and(u>=a,u<=b),1,and(u>b,u<a),0)
"                                                      range(12,1,99)"
range(12,11,20)
""
"f(a,b)=        test(a>=b,"
"                  print(a),"
"               b>5,"
"                  print(b),"
"                  print(a + b))"
""
f(a,b)=        test(a>=b, print(a), b>5, print(b), print(a + b))
"                                                       f(1,2)"
f(1,2)
""
" check(predicate)   stop if predicate False, else continue"
" check(a>b)"
check(1==1)
""
" number(x)        True if x a number"
"                                                       number(333)"
number(3.14116)
" isinteger(x)     True if x an integer"
"                                                       isinteger(123)"
isinteger(123)
""
"*******************************************************************************"
"*                                loop instruction :                           *"
"*******************************************************************************"
" for(i, min, max, do (...))"
" for(i, 1, 10,  print(i^2))"
for(i, 1, 10,  print(i^2))
""
" A=zero(3,5)"
" for(i, 1, 3,"
"   for(j, 1, 5,"
"      A[i,j] = i + 2*j ))"
A=zero(3,5)
for(i, 1, 3, for(j, 1, 5, A[i,j] = i+2*j ))
A
""
" product(k, imin ,imax, f)"
" Pi = 2 * product(k, 1 ,10 ,(2*k)^2/((2*k)^2-1))"
Pi=2*product(k, 1 ,200 ,(2*k)^2/((2*k)^2-1))
Pi
Pi=float(Pi)
Pi
""
" sum(k, imin, imax ,f)"
""
" S = sum (k, 1, 15, 1 + 1/k)"
S=sum(k,1,15, 1 + 1/k)
S
" S1 = sum(k, 1, 10, 1/4*(k^2*(k+1)^2))"
S1=sum(k, 1, 10, 1/4*(k^2*(k+1)^2))
S1
""
"*******************************************************************************"
"*                                 recurrences :                               *"
"*******************************************************************************"
" a(n) = test(n==0, 0, sqrt(2 + a(n-1)))"
a(n) = test(n==1, 1, sqrt(2 + a(n-1)) )
"                                                 a(9)"
a(9)
float(a(9))

""
" defining functions :"
""
" f(arg1, arg2, ..., argn) = expression"
""
" f(arg1, arg2, ..., argn) = do (expression1, expression2,..., expressionk)"
" return expressionk "
" f(x,y) = sin(x)^2  + 2*cos(y)"
f(x,y) = sin(x)^2  + 2*cos(y)
"                                                 float(f(3,1))"
float(f(3,1))
""
" g(x,y) = do( check(x > 2), check(y > 0), x^3 + 3/y)"
g(x,y) = do( check(x < 0), check(y > 0), x^3 + 3/y)
"                                                 g(-5,2)"
g(-5,2)
""
" abs(x) = test(x < 0, -x, x)"
abs(x) = test(x < 0, -x, x)
"                                                  abs(4-3*5)"
abs(4-3*5)
""
" power(x, n) = do( check(n > 0), p=1, for (i, 1, n, p = p * x), p)"
power(x, n) = do( check(n > 0), p=1, for (i, 1, n, p = p * x), p)
"                                                  power(2,5)"
power(2,5)
""
" fib(n) =  do( check(n > 0), test(n==1, 1, test(n==2, 2, fib(n-1) + fib(n-2) )))"
fib(n)= do( check(n > 0), test(n==1, 1, test(n==2, 2, fib(n-1) + fib(n-2) )))
"                                                  fib(15)"
fib(15)
""
" fac(n) = do( check (n > 0),"
"            test(n == 1, 1, n * fac(n-1) ))"
fac(n) =   do( check (n > 0),
test(n == 1, 1, n * fac(n-1) ))
"                                                   fac(15)"
fac(15)
" mutually recurcive functions :"
""
" even(n) = test(n==0, 1, odd(n-1))"
" odd(n) = test(n==0, 0, even(n-1))"
even(n) = test(n==0, 1, odd(n-1))
odd(n) = test(n==0, 0, even(n-1))
"                                                   even(3), odd(3)"
even(3)
odd(3)
""
" x(n) = test(n==0, 1, x(n-1) + y(n-1))"
" y(n) = test(n==0, 2, y(n-1) * x(n-1))"
x(n) = test(n==0, 1, x(n-1) + y(n-1))
y(n) = test(n==0, 2, y(n-1) * x(n-1))
"                                                   x(3), y(3)"
x(3)
y(3)
""
"*******************************************************************************"
"*                                    graphics :                               *"
"*******************************************************************************"
x=quote(x)
" draw(x^2-4)"
draw(x^2-4)
""
" xrange = (-20,20)"
" yrange = (-1,1)"
" draw(sin(x)/x)"
xrange = (-20,20)
yrange = (-1,1)
draw(sin(x)/x)
""
" draw(x*sin(x))"
" yrange = (-20,20)"
yrange = (-20,20)
draw(x*sin(x))
""
" xrange = (-10,10)"
" yrange = (-50,50)"
" f=exp(-x^2/2)*hermite(x,5)"
" draw(f)"
xrange = (-10,10)
yrange = (-50,50)
f=exp(-x^2/2)*hermite(x,5)
draw(f)
""
" trange = (-pi,pi)"
" xrange = (-1,1)"
" yrange = (-1,1)"
" f=(cos(t),sin(t))"
" draw(f)"
trange = (-pi,pi)
xrange = (-1,1)
yrange = (-1,1)
f=(cos(t),sin(t))
draw(f)
""
" trange = (-2*pi,2*pi)"
" xrange = (-5,5)"
" yrange = (-5,5)"
" x = cos(t)/(1+sin(t)^2)"
" y = sin(t)*cos(t)/(1+sin(t)^2)"
" draw(5*(x,y))"
trange = (-2*pi,2*pi)
xrange = (-5,5)
yrange = (-5,5)
x = cos(t)/(1+sin(t)^2)
y = sin(t)*cos(t)/(1+sin(t)^2)
draw(5*(x,y))
""
" trange = (-29,29)"
" xrange = (-15,15)"
" yrange = (-15,15)"
" x = 5.5*cos(t)+8*cos(-5.5*t/4.5)"
" y = 5.5*sin(t)+8*sin(-5.5*t/4.5)"
trange = (-29,29)
xrange = (-15,15)
yrange = (-15,15)
x = 5.5*cos(t)+8*cos(-5.5*t/4.5)
y = 5.5*sin(t)+8*sin(-5.5*t/4.5)
draw((x,y))

""
" trange = (pi,pi)"
" xrange = (-1,1)"
" yrange = (-1,1)"
" r=(1+cos(t))/2"
" u=(cos(t),sin(t))"
" draw(r*u)"
trange = (-2*pi,2*pi)
xrange = (-1,1)
yrange = (-1,1)
r=(1+cos(t))/2
u=(cos(t),sin(t))
draw(r*u)
""

"*******************************************************************************"
"*                               linear algebra :                              *"
"*******************************************************************************"
do(a=quote(a),b=quote(b),c=quote(c),d=quote(d),x=quote(x),y=quote(y))
"                                            A=((a,b),(c,d)) a 2x2 matrix"
A=((a,b),(c,d))
A
"                                            A[1,2]    element row 1 column 2"
A[1,2]
"                                            inv(A)"
inv(A)
"                                            det(A)"
det(A)
"                                            transpose(A)"
transpose(A)
"                                            v=(x,y)      a 2x1 vector"
v=(x,y)
v
"                                            abs(v)"
abs(v)
"                                            dot(v,v)     scalar product v.v"
dot(v,v)
"                                            dot(A,v)     = v1"
dot(A,v)
"                                            dot(v,A,v)   scalar product v.v1 "
dot(v,A,v)
""
B=((0,1),(1,0))
B
"                                            eigenvec(B)"
eigenvec(B)
"                                            eigenval(B)"
eigenval(B)
""
"                                            A=((0,1),(-1,0))"
"                                            I=unit(2)"
"                                            A^2+3*A-I"
A=((0,1),(-1,0))
I=unit(2)
A^2+3*A-I
"*******************************************************************************"
"*                                3D  vectors and tensors :                    *"
"*******************************************************************************"
" variables used to describe physical quantities are of a number of types,"
" including scalars, vectors, and tensors"
" scalars are used to represent physical quantities with no directional qualities"
" such as temperature, volume, and time. Vectors are used for quantities which"
" have a single directional quality such as velocity and force. Tensors (we"
" consider only second-order tensors) are associated with quantities which have"
" two directional characteristics "
""
" vectors and Vector Operations :"
""
" given a coordinate system in three dimensions, a vector may thus be represented"
" by an ordered set of three components which represent its projections"
" v1; v2; v3 on the coordinate axes 1; 2; 3 "
" v = (v1, v2, v3)"
v = (v1, v2, v3)
v
" the three most commonly used coordinate systems are rectangular, cylindrical,"
" and spherical. Alternatively, a vector may be represented by the sum of the"
" magnitudes of its projections on three mutually perpendicular axes :"
v = v1*delta1 + v2*delta2 + v3*delta3
""
" the unit vectors"
delta1  delta2  delta3
" are in the rectangular coordinate system :"
deltax deltay deltaz
" in cylindrical coordinates :"
deltar  deltatheta deltaz
" in spherical :"

deltar  deltatheta  deltaphi
""
" each of these unit vectors points in the direction of the indicated spatial"
" coordinate"
" the magnitude of a vector is given by :"
" v = sqrt(v1^2 + v2^2 + v3^2)     =  abs(v)"
sqrt(v1^2 + v2^2 + v3^2)
" addition and subtraction of vectors is easily executed :"
"v + w = (v1 + w1)delta1 + (v2 + w2)delta2 + (v3 + w3)delta3"
(v1 + w1)delta1 + (v2 + w2)delta2 + (v3 + w3)delta3
" as well as multiplication by a scalar :"
" c*v = (c*v1)delta1 + (c*v2)delta2 + (c*v3)delta2"
(c*v1)delta1 + (c*v2)delta2 + (c*v3)delta2
""
" the dot product of two vectors results in a scalar :"
" v.w = v1 w1 + v2 w2 + v3 w3    =  dot(v,w)"
do(v=(v1,v2,v3), w=(w1,w2,w3))
dot(v,w)
""
" tensors and tensors operations :"
""

" a tensor is represented by an ordered array of nine components:"
" T = ((T11,T12,T13),(T21,T22,T23),(T31,T32,T33))"
T = ((T11,T12,T13),(T21,T22,T23),(T31,T32,T33))
T
" the diagonal elements of a tensor are those which have two identical subscripts"
" while the other elements are termed nondiagonal"
" the transpose of a tensor is obtained by interchanging the subscripts on"
" each element:"
TT=transpose(T)
TT
" a tensor is described as symmetric when T = transpose(T)"
" one special tensor is the unit tensor:"
delta=unit(3)
delta
" the dyadic product of two vectors results in a tensor, as follows:"
"                             v=(v1,v2,v3)"
"                             w=(w1,w2,w3)"
"                             D=outer(v,w)"
v=(v1,v2,v3)
w=(w1,w2,w3)
D=outer(v,w)
D
" this leads to the definition of the unit dyads, of which there are nine:"
do(e1=(1,0,0),e2=(0,1,0),e3=(0,0,1))
delta11=outer(e1,e1)
delta11
delta12=outer(e1,e2)
delta12
delta13=outer(e1,e3)
delta13
delta21=outer(e2,e1)
delta21
" in a similar manner to vectors, tensors are easily added  e.g.  A = T + U:"
U = ((U11,U12,U13),(U21,U22,U23),(U31,U32,U33))
T+U

" or multiplied by scalars e.g   c * T:"
c*T
""
" the double dot product of two tensors results in a scalar e. g.   T:U    :"
T[1,1]*U[1,1]+T[1,2]*U[1,2]+T[1,3]*U[1,3]+T[2,1]*U[2,1]+T[2,2]*U[2,2]+T[2,3]*U[2,3]+T[3,1]*U[3,1]+T[3,2]*U[3,2]+T[3,3]*U[3,3]
""
" the dot product of a tensor with a vector is  e.g   T.v   :"
dot(T,v)
""
" in general, (T.v) not= (v.T), however, they are equal if T is symmetric"
" the magnitude of a tensor is defined as :"
" |T| = sqrt(T:transpose(T))"
""
" contraction  :"
" contract(T,1,2) (1 for i, 2 for j), is the scalar  sum(A[i,j],with i=j) : "
" contract(T,1,2)"
contract(T,1,2)
" sum(i,1,3,(sum(j,i,i,T[i,j]))"
sum(i,1,3,(sum(j,i,i,T[i,j])))
" rank  :"
" return the number of indices that the tensor has  e.g.  rank(T) :"
rank(T)

"*******************************************************************************"
"*                                  calculus :                                 *"
"*******************************************************************************"
" derivatives :"
""
" d(f,x)                       -->   df(x)/dx"
" d(f,x,x)                   -->   d²f(x)/dx²"
" d(f,x,y)                   -->   d²f(x,y)/dxdy"
" d(f,x,m,y,n,...)"
" examples :"
"                                                            d(x^3)"
d(x^3)
"                                                            d(x^3,x,x)"
d(x^3,x,x)
"                             f=3*x^(2)*y-4*x^3              d(f,x,y)"
f=3*x^(2)*y-4*x^3
d(f,x,y)
"                             f=-4*x^(2)*y+4*x^3             d(f,y,x,x)"
f=-4*x^(2)*y+4*x^3
d(f,y,x,x)
"                             w=x*y^2*exp(x*y)               d(w,y)"
w=x*y^2*exp(x*y)
d(w,y)
""
" the gradient is obtained by using a vector for x in d(f,x)"
"                             r=sqrt(x^2+y^2)                d(r,(x,y))"

r=sqrt(x^2+y^2)
d(r,(x,y))
""
" the f in d(f,x) can be a tensor function. Gradient raise the rank by one"
"                             F=(x+2y,3x+4y) X=(x,y)         d(F,X)"

F=(x+2y,3x+4y)
X=(x,y)
d(F,X)
-----------------------
" integral  :"
" integral(f,x) returns the integral of f with respect to x"
" a multi-integral can be obtained by extending the argument list"
integral(x^2)
integral(x*y,x,y)

" defint(f,x,a,b,...) computes the definite integral of f with respect to x"
" evaluated from a to b. The argument list can be extended for multiple integrals"
" the following example computes the integral of f = x^2 over the domain of a"
" semicircle. For each x along the abscissa, y ranges from 0 to sqrt(1-x^2)"
""
defint(x^2,y,0,sqrt(1-x^2),x,-1,1)
""
" as an alternative, the eval function can be used to compute a definite"
" integral step by step :"
""
" I=integral(x^2,y)"
" I=eval(I,y,sqrt(1-x^2))-eval(I,y,0)"
" I=integral(I,x)"
" eval(I,x,1)-eval(I,x,-1)"
I=integral(x^2,y)
I=eval(I,y,sqrt(1-x^2))-eval(I,y,0)
I=integral(I,x)
eval(I,x,1)-eval(I,x,-1)
""
" here is a useful trick. Difficult integrals involving sine and cosine can"
" often be solved by using exponentials. Trigonometric simplifications involving"
" powers and multiple angles turn into simple algebra in the exponential domain"
""
" f=sin(t)^4-2*cos(t/2)^3*sin(t)"
" f=circexp(f)"
" defint(f,t,0,2*pi)"
f=sin(t)^4-2*cos(t/2)^3*sin(t)
f=circexp(f)
defint(f,t,0,2*pi)
" here is a check "
" g=integral(f,t)"
" f-d(g,t)"
g=integral(f,t)
f-d(g,t)

"*******************************************************************************"
"*                               build in functions :                          *"
"*******************************************************************************",
" usual :"
" sin(x), cos(x), tan(x), asin(x), acos(x), atan(x),"
" sinh(x), cosh(x), tanh(x), asinh(x), acosh(x), atanh(x)"
" circexp(x), expcos(x), expsin(x)"
" exp(x), sqrt(x), log(x)"
""
" logical :"
" and(a,b,c,...), or(a,b,c,...), not(a)"
""
" arithmetic :"
" mod(a,b)"
" ceiling(x), floor(x),"
" gcd(a,b,...), lcm(a,b,...), choose(n,k),"
" factor(n), isprime(n), prime(n),"
" factorial(n), n!,"
" divisors(x)"
""
" complex number :"
" mag(z), arg(z), conj(z),"
" imag(z), real(z),"
" polar(z), rect(z), clock(z)"
""
" polynomials :"
" coeff(p,x,n), deg(p,x), leading(p,x),"
" expand(r,x),  factor(p,x,y,...), quotient(p,q,x), rationalize(x),"
" roots(p,x), nroots(p,x),"
" hermite(x,n), laguerre(x,n,a), legendre(x,n,m)"
""
" vector, matrix, tensor :"
" rank(a), dim(a,n),"
" unit(n), zero(i,j,...),"
" det(m), inv(m), transpose(m), cofactor(m,i,j),"
" contract(a,i,j), transpose(a,i,j), dot(a,b,...), inner(a,b,...), outer(a,b,...),"
" abs(u), cross(u,v), dot(u,v), inner(u,v), outer(u,v),"
" hilbert(n)"
""
"*******************************************************************************"
"*                                    calculus :                               *"
"*******************************************************************************"
" d(f,x) :"
" d(f,x,x), d(f,x,n),"
" d(f,x,y), d(f,x,m,y,n,...)"
" d(r,(x,y))"
" div(u)"
" curl(u)"
""
" integral(f,x), integral(f,x,y,...)"
" defint(f,x,a,b,...)"
""
" Taylor series :"
" taylor(f,x,n,a)"
""
"*******************************************************************************"
"*                                     misclaneous :                           *"
"*******************************************************************************"
""
" abs(z)
" print(a,b,...)"
" simplify(x), condense(x), numerator(x), denominator(x), rationalize(x)"
" draw(f,x), xrange, yrange"
" sgn(x)"
" erf(x), erfc(x)"
" subs(a,b,c), eval(f,x,n), filter(f,a,b,...)"
" float(x)"
" quote(x)"
" clear"
" stop"
" tty = 0 o full display on"
" tty = 1 o full display off"
" trace = 0/1  trace off/on"
" autoexpand = 0/1  off/on"
" last        use last expression as argument"
" bake = 0/1"
" selftest"
""
" special functions :"
" hermite(x,n), laguerre(x,n,a), legendre(x,n,m)"
" dirac(x)"
""
"*******************************************************************************"
"*                                 EVA script functions :                      *"
"*******************************************************************************"
" isMvector(u)"
" isVector(u)"
" isScalar(u)"
" isInteger(u)"
" isNumber(u)"
""  