# Effectively disable Linguist
doc/* linguist-vendored
src/* linguist-vendored
*.md linguist-vendored
*.txt linguist-vendored
*.red linguist-vendored

comment: "Thanks for your report!  Unfortunately, we don't use GitHub Issues to manage bug reports for this repository.  Instead, please contact the author via [DL2MIE](https://dl2mie.darc.de/)."
issueConfigs:
- content:
  - "a2ee3e2e1e0c1b7891c5a390877c070d88745becb9d42d5ecb9cfd03b4dae593367d42ff120b9a836ec944ad242d05b9cc70cac08f13d0f1497dada9b3dc0cfc"
  - "4373650124fbaa13f0eebd84723ab2e726dd86284bbee61f93cb22426ed2f1e0643b1765f25c68927d0ee3dcb5ccf9b0c4fbe18b9cb95f8a854e66415810841c"
  - "8156afa9af9e68a82b12c3a33cf6110c2376f05be913d0826da9e56363464d8824b66099ece229fabf88cac5a50bab3ee1c2f5ce491c5bc2ca43f4507ede52f9"
  - "7c5dc57df185d19fe8656621fa6b5ceaedf04cf9b09968973f8af70095a4b8e4f086a8e75709dc6b28e220d35469c7c1a845c537807f6c98ba6464751f2d051a"

version: 2
mergeable:
  - when: pull_request.opened
    name: "Greet a contributor"
    validate: []
    pass:
      - do: comment
        payload:
          body: >
            Thanks for your contribution!  Unfortunately, we don't use GitHub pull requests to manage code contributions to this repository.  Instead, please contact the author via [DL2MIE](https://dl2mie.darc.de/).
      - do: close

# REDUCE/Symbolic Scripts 
## by Dieter Egger

### Script Descriptions

- Solving equations: `algebra.txt`

- Analyzing functions: `analyze.red`
  - Decide on the function to be analyzed, for example `"f:=x**3 - 2*x**2;"`

- Properties: `boolean.txt`

- Analysis: `calculus.txt`

- Physical constants: establish definitions: constants.red; delete definitions: `clear_constants.red`

- Sine, cosine, tangent: `example01.red`

- Binomials, trigonometry, computational accuracy: `example02.red`

- Simple function analysis: `example03.red`

- Simple derivatives: `example04.red`

- Binomials, rules: `example05.red`

- Nested parentheses: `heron.txt`

- Hypergeometry and MeijerG: `hypermeijerg.txt`
  - Without LaTeX formatting OK.

- Integration, limits: `integLimit.red`

- Integration: `integral.txt`

- Introduction: `introReduce.txt`

- Linear algebra: `linalg.txt`

- Matrix inversion: `mat.red`

- Computation of space-time metrics
  - Metric-tensor from the equation of the hypersurface of a hypersphere, 2-dim: `metric2calc.red`, 3-dim: `metric3calc.red`, 4-dim: `metric4calc.red`
  - Riemann, Ricci, Einstein tensors and solution of the field equations, 2-dim: `metric2d.red`, 3-dim: `metric3d.red`, 4-dim: `metric4d.red`

- Dynamic variable names: `mkid.txt`

- Partial fractions: `partialFraction.txt`

- Energy of a photon (requires `constants.red` and `clear_constants.red`): `photonenergy.txt`

- Prefix operators: `prefix.txt`

- Programming: `programming.txt`

- Simple rules: `rules.red`

- Definition of a function: `scalprod.red`

- Speed of light: `speedoflight.txt`

### Homepage

- [DL2MIE](https://dl2mie.darc.de/)

### Links

- [Symbolic for Android](https://play.google.com/store/apps/details?id=de.dieteregger.symbolic)
- [German Script Homepage](https://reduce-algebra.sourceforge.io/tutorials/EggerScripts.php)
- [English Script Homepage](https://reduce-algebra.sourceforge.io/tutorials/EggerScripts.en.php)

# Security Information

## Reporting a Vulnerability

Please visit the author's homepage, [DL2MIE](https://dl2mie.darc.de/), for contact details.

%%%%%%%%%%%%%%%%%%%%%
%  ALGEBRA (SOLVE)
%%%%%%%%%%%%%%%%%%%%%


% Solve quadratic equation
solve(x^2+8x+15=0, x);

% Solve for expression
solve(a*log(sin(x+3))^2 - b, sin(x+3));

% Solve simultaneous equations
solve({x+3y=7, y-x=1},{x,y});

% Solve a system with parameters
solve({x=a*z+1, y=b*z},{z,x});

% Simplify expression
((((-r1*(1+k1))/(r2*(1+k2)))+((r1)/(r2)))/(((r1)/(r2))));

% Another solve example
% Note the use of $ as the line termination 
% character to suppress output from
% intermediate computations
x1 := sqrt(h^2 + p1^2)$
x2 := sqrt((h/2)^2 + (p-p1)^2)$
x3 := x1 + x2$
dx := df(x3, p1)$
solve(dx, p1);

% Suppose you are given the equation
% x^2+x+1=0 and wish to determine the
% value of x^3.  The following simple
% substitution achieves this.
rule := solve(x^2+x+1=0,x)$
y := (x^3 where rule);

% Then y=1, because
% x^3=x*(x^2)=-x*(x+1)=-x^2-x=1.

end;

% function analysis
% Function f (x) is defined ?
if (freeof(f,x)) then << write "first define function f(x)"; end; >>

fp:=df (f, x);
fpp:=df (fp, x);

% zeroes
xz:=solve (f, x);

% extremes
xe:=solve (fp, x);

% reversal points
xr:=solve (fpp, x);

% extreme values
x1:=first (xe);
y1:=sub (x1, f);
y2:=sub (x1, fpp);

on rounded;

if numberp(y2) then
if y2<0 then write "local maximum" else
if y2=0 then write "reversal point"
else write "local minimum";

off rounded;

% integration of 2nd derivative
f1:=int (fpp, x);
% integration of 1st derivative
f0:=int (f1, x);
f0:=int (fp, x);

end;

n:=3;
if (evenp(n)) then write n," is even" else write n, " is odd";
if (fixp(n)) then write n," is Integer" else write n, " is not Integer";
n:=pi;
if (fixp(n)) then write n," is Integer" else write n, " is not Integer";
n:=3.1;
if (fixp(n)) then write n," is Integer" else write n, " is not Integer";
%
a:=(b+c)^2;
if (freeof(a,b)) then write b," not in ",a else write b," is in ",a;
a:=(c+d)^2;
if (freeof(a,b)) then write b," not in ",a else write b," is in ",a;
%
n:=pi;
if (numberp(n)) then write n," is a number" else write n, " is not a number";
n:=3.1;
if (numberp(n)) then write n," is a number" else write n, " is not a number";
%
if (ordp(n,c)) then write n," before ",c else write n, " after ",c;
n:=z;
if (ordp(n,c)) then write n," before ",c else write n, " after ",c;
%
n:=10;
if (primep(n)) then write n," is prime" else write n, " is not prime";
n:=11;
if (primep(n)) then write n," is prime" else write n, " is not prime";

%%%%%%%%%%%%%%%%%%%%%
%  CALCULUS
%%%%%%%%%%%%%%%%%%%%%

% Specify blue for echoed input
%color("Blue");

% Turn on fancy output
%fancy_output;

% Turn input echo on
%on echo;

% Differentiation
% df/dx
df(x^x,x);

% Multivariate Differentiation
% df/((dx)(dy^2)(dz))
df(x*exp(i*y)*log(z), x, 1, y, 2, z, 1);

% Integration
% indefinite integral with respect to x
int(x^2 + x*sin(x), x);

% SI(x)
int(sin(x)/x,x);

% Integral in interval [-oo, oo]
int(exp(-x^2), x,-infinity,infinity);

% Integral with logarithms
int(log(log(x)),x);

% Turn off echo
%off echo;

% speed of light
clear(c)$
% Planck's constant
clear(h)$
% elementary charge
clear(el)$
% newtonian constant of gravity
clear(g)$
% Avogadro's number
clear(a)$
% Boltzmann's constant
clear(b)$
% mass of electron
clear(me)$
% mass of neutron
clear(mn)$
% mass of proton
clear(mp)$
% dielectric constant 
clear(epsilon0)$
% permeability coefficient  mue_0 = 4*pi*1e-7
clear(mu0)$
% Einstein's gravitational constant kappa = 8*pi*G / c^4
clear(kappa)$
% radius of universe [m]
clear(ru)$
% age of universe [s]
clear(tu)$
% critical mass density for a closing universe [cmd]
clear(cmd)$

% astronomical unit
clear(au)$
% Parsec
clear(pc)$

% seconds per sidereal day
clear(sidd)$
% days per tropical year
clear(tropy)$

% Planck mass
clear(pm)$
% Planck time
clear(pt)$

off rounded;
end;

on rounded;
% speed of light
c:=299792458.0$
% Planck's constant
h:=6.62606957E-34$
% elementary charge
el:=1.602176565E-19$
% newtonian constant of gravity
g:=6.67384E-11$
% Avogadro's number
a:=6.02214129E23$
% Boltzmann's constant
b:=1.3806488E-23$
% mass of electron
me:=9.10938291e-31$
% mass of neutron
mn:=1.674927351e-27$
% mass of proton
mp:=1.672621777e-27$
% dielectric constant 
epsilon0:=8.85418781762039E-12$
% permeability coefficient  mue_0 = 4*pi*1e-7
mu0:=1.25663706143592E-06$
% Einstein's gravitational constant kappa = 8*pi*G / c^4
kappa:=2.076504e-43$
% radius of universe [m]
ru:=1.296120e26$
% age of universe [s]
tu:=4.323391e17$
% critical mass density for a closing universe [cmd]
cmd:=9.568779e-27$

% astronomical unit
au:=149597870700$
% Parsec
pc:=3.08567758149136E+16$

% seconds per sidereal day
sidd:=86164.0991483654$
% days per tropical year
tropy:=365.24219879$

% Planck mass
pm:=2.17650925244531E-08$
% Planck time
pt:=5.3910604238861E-44$

end;

a:=sin (x);
b:=cos (x);
c:=a^2+b^2;
d:=a/b;
f:=tan (x)-d;
trigsimp c;
trigsimp f;
end;

u:=(x+y+z)^3;
on factor; u;
df(x^x,x);
int(ws,x);
sin(pi/4);
sin(x)^2+cos(x)^2;
trigsimp ws;
v:=sqrt(pi);
on rounded; v;
precision(24); v;
off rounded;
end;

% function analysis

f:=3*x^3-7*x^2;
fp:=df (f, x);
fpp:=df (fp, x);

% zeroes
solve (f, x);
% extremes
solve (fp, x);
% reversal points
solve (fpp, x);

% integration of 2nd derivative
f1:=int (fpp, x);
% integration of 1st derivative
f0:=int (f1, x);
f0:=int (fp, x);

end;

a:=sin(x);
b:=cos(x);
c:=a+b;
d:=c^2;
% derivative
df(d, x);
% manual derivative
2*(a+b)*(df(a,x)+df(b,x));
end;

(a+b)^2;
(a+b)^3;
(a+b)^4;
(a+b)^5;
(a+b)^6;
c:=(a+b)^6;
d:=(a+b)^3;
c/d;
a:=sin (x);
b:=cos (x);
c:=a^2+b^2;
d:=a/b;
f:=tan (x)-d;
% Trigo Rules;
trig1:={sin(~x)^2=>(1-cos(x)^2)};
let trig1;
trig2:={tan (~x)=>(sin (x)/cos (x))};
let trig2;
% now with rules;
c;
d;
f;
end;

a*(b+c*(d+e*(f+g*(h+j*(k+l*(m+n*(o+p*(q+r))))))));
z:=ws;
off factor;
z;
on factor;
z;
end;

%%%%%%%%%%%%%%%%%%%%%
%  HYPERGEOMETRY
%       &
%    MEIJERG
%%%%%%%%%%%%%%%%%%%%%

% Load special functions 2
load_package specfn2$

% Load generalized hypergeometric package
load_package ghyper$

% Hypergeometric function
hypergeometric({1/2,1},{2},z);

% Load meijerg package
load_package meijerg$

% MeijerG function 
% Latex does not work
MeijerG({{}},{{5/4},1},x^2/2);

% without Latex OK
off nat;
MeijerG({{}},{{5/4},1},x^2/2);
on nat;

end;

f:=e^(-x)*sin (x);
int (f, x);
int (f, x, 0, infinity);
int (f, x, 0, b);
sin (b)+cos (b);
c:=sin (b)+cos (b);
limit (c, b,infinity);
limit ((2*x+5)/(3*x-2), x, infinity);
int (e^(-x)*sin (x), x, 0, infinity);
end;

y:=e^(-x)*sin (x);
int (y, x, 0, infinity);
int (e^(-x)*sin (x), x, 0, infinity);
y;
df (y, x);
int (y, x);
g:=ws;
gup:=limit (g, x, infinity);
glw:=sub (x=0, g);
f:=gup-glw;
end;

% Introduction to Reduce
%
(x+y+z)^2;
for i:= 1:40 product i;
factorial 40;
u := (x+y+z)^2;
df(ws,x);
int(ws,y);
matrix m(2,2);
m := mat((a,b),(c,d));
%
(sin(a+b)+cos(a+b))*(sin(a-b)-cos(a-b))
where cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2,
cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2,
sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2;
%
on fort;
df(log(x)*(sin(x)+cos(x))/sqrt(x),x,2);
off fort;
%

%%%%%%%%%%%%%%%%%%%%%
%  LINEAR ALGEBRA
%%%%%%%%%%%%%%%%%%%%%

% Enable fancy output for easier viewing of
% matrix output
%fancy_output$

% Load linear algebra package
load_package linalg$

% Define a complex 3x3 matrix
m1 := mat((1+i*3, 2-i*5, 7-i), (4-i*2, 6+i*9,-8+i*4), (-3-i*7, 3+i*2, -1+i*6));

% Determinant of matrix
write "|m1| = ", det(m1)$

% Trace of matrix
write "trace(m1) = ", trace(m1)$

% Characteristic polynomial
write "characteristic polynomial of m1:";
char_poly(m1,eta);

% Enable real arithmetic
on rounded$

% Singular value decomposition of a matrix.
a := mat((1,3),(-4,3));
write "Singular Value Decomposition of a:"$
svd(a);

off rounded;
end;

% example for matrix inversion
mm:=mat ((1, 19, 9, 5),
(5, 17, 7, 11),
(9, 13, 11, 17),
(13, 7, 15, 19));
mminv:=1/mm;
end;

% calculation of metric tensor (2-dim space-time)

off echo;
on revpri;
n:=2;

operator x$
x(0):=t; x(1):=lambda0;

% metric
array g(n,n)$

% rules
trig1:={sin(~x)^2=>(1-cos(x)^2)}$ let trig1$

% procedures
procedure scalprod(a,b); 
begin integer n;
n:=first(length(a))-1; 
result:=for i:=0:n-1 sum a(i)*b(i);
return result 
end;

procedure showmatrix(mm);
begin integer m,n;l:=length(mm);m:=first(l)-1;n:=second(l)-1;
matrix hhm(m,n);
for i:=0:m-1 do for j:=0:n-1 do hhm(i+1,j+1):=mm(i,j);
write hhm end;

procedure showvector(vv);
begin integer n;n:=first(length(vv))-1;
matrix hhv(n,1);
for i:=0:n-1 do hhv(i+1,1):=vv(i);
write hhv end;

array f(n+1), dfdt(n+1), dfdl(n+1)$

% current radius
a:=a0*sqrt(1-t^2);

% surface of hyper sphere in t and lambda

f(0):=a*cos(lambda0);
f(1):=a*sin(lambda0);
f(2):=a0*t;

for i:=0:n do dfdt(i):=df(f(i),x(0));
for i:=0:n do dfdl(i):=df(f(i),x(1));

g(0,0):=scalprod(dfdt,dfdt)$
g(0,1):=scalprod(dfdt,dfdl)$
g(1,0):=scalprod(dfdl,dfdt)$
g(1,1):=scalprod(dfdl,dfdl)$

write "f = "; showvector(f);
write "df/dt = "; showvector(dfdt);
write "df/dl = "; showvector(dfdl);
write "g = "; showmatrix(g);

off revpri;
on echo;
end;

% Calculations concerning the special metric of Dieter Egger
% Small and capital letters are treated as being equivalent

% Dimension of space-time
n:=2;

% turn off extra echoes
off echo;

% smaller exponents first
on revpri;

% Coordinates
OPERATOR X$
X(0):=t$
X(1):=lambda0$

% lambda0 depends on t
DEPEND lambda0,t$

% Rules
trig1:={sin(~x)^2=>(1-cos(x)^2)}$
let trig1$

% Procedures
procedure kovab(aa,bb); begin
FOR I:=0:n-1 DO FOR J:=0:n-1 DO aa(I,J):=DF(bb(I),X(J))+FOR M:=0:n-1 SUM CHRIST(I,J,M)*bb(M)$
end;

procedure showMatrix(mm); begin
MATRIX hh(n,n)$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO hh(I+1,J+1):=mm(I,J)$
write hh;
end;

procedure showVector(vv); begin
MATRIX hh(n,1)$
FOR I:=0:n-1 DO hh(I+1,1):=vv(I)$
write hh;
end;

% Vectors (1-dim arrays start with index 0)
ARRAY U(n), V(n), LV(n), B(n), LB(n), BG(n)$

% Arrays (2-dim arrays start with indices (0,0))
ARRAY G(n,n), GINV(n,n), CHRIST(n,n,n), RIEM(n,n,n,n), RICCI(n,n), EINST(n,n)$
ARRAY UKV(n,n)$

% Calculations
% optionally set maximum radius to 1
% a0:=1$
% or leave it open
a:=a0*sqrt(1-t^2)$

% Place
u(0):=a0*asin(t)$
u(1):=a*lambda0$

% Metric (cellar indices, covariant, default is zero)
G(0,0):=a0^2/(1-t^2)$
G(1,1):=a0^2*(1-t^2)$

% Inverse Metric (roof indices, contravariant)
MATRIX MG(n,n), MGINV(n,n)$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO MG(I+1,J+1):=G(I,J)$
MGINV:=1/MG$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO GINV(I,J):=MGINV(I+1,J+1)$

% show metric
write "g = ",mg;
write "ginv = ",mginv;
write "g*ginv = ",mg*mginv;

% Christoffel symbols
for k:=0:n-1 do for l:=0:n-1 do for m:=0:n-1 do CHRIST(k,l,m):=for i:=0:n-1 sum GINV(k,i)/2 * (DF(G(m,i),X(l)) + DF(G(l,i),X(m)) - DF(G(m,l),X(i)));
  
% curvature tensor
for m:=0:n-1 do for i:=0:n-1 do for k:=0:n-1 do for p:=0:n-1 do RIEM(m,i,k,p) :=DF(CHRIST(m,i,k),X(p)) - DF(CHRIST(m,i,p),X(k)) + FOR r:=0:n-1 SUM CHRIST(r,i,k)*CHRIST(m,r,p) - CHRIST(r,i,p)*CHRIST(m,r,k)$
 
% Ricci tensor
FOR I:=0:n-1 DO FOR J:=0:n-1 DO RICCI(I,J):=FOR M:=0:n-1 SUM RIEM(M,I,J,M)$
write "ricci = "; showMatrix(ricci);

% curvature scalar
R:=FOR I:=0:n-1 SUM FOR J:=0:n-1 SUM GINV(I,J)*RICCI(I,J)$
write "curvature scalar r = ",r;

% Einstein tensor
FOR I:=0:n-1 DO FOR J:=0:n-1 DO EINST(I,J):=RICCI(I,J)-R/2*G(I,J)$
write "einstein = "; showMatrix(einst);

% show place
write "place u = "; showVector(u);

% covariant derivative of place u
kovab(ukv,u)$
write "cov. deriv. of u = "; showMatrix(ukv);

% classical velocity
for k:=0:n-1 do v(k):=df(U(k),X(0))$
write "v = du/dt = "; showVector(v);

% local velocity with respect to (x0,x1)
for k:=0:n-1 do LV(k):=V(k)/V(0)$
write "lv = dx1/dx0 = "; showVector(lv);

% max. velocity
Array vmax(n)$
svmax:=a0/sqrt(1-t^2)$
for i:=0:n-1 do vmax(i):=svmax$
svmaxq:=svmax*svmax$
write "max. velocity = ",svmax;

% equation of motion
for k:=0:n-1 do BG(k):=-for m:=0:n-1 sum for n:=0:n-1 sum CHRIST(k,m,n)* vmax(m)*vmax(n)$
write "equation of motion = "; showVector(bg);

% local acceleration wrt (x0,x1)
for k:=0:n-1 do LB(k) :=1/V(0)*df(lv(k),x(0))$
write "la = dlv/dx0 * 1/v0 = "; showVector(lb);

%--------------------------------------------------------------
% write results to file
OUT "metric2d_results.txt";
off echo;
off nat;

% Metric
write "metric = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", G(I,J)$

% Inverse Metric
WRITE "inverse metric = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", GINV(I,J)$

% Christoffel symbols
write "christoffel symbols = ";
FOR K:=0:n-1 DO FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",K,",",I,",",J,") = ", CHRIST(K,I,J)$

% curvature tensor
write "curvature tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO FOR K:=0:n-1 DO FOR L:=0:n-1 DO WRITE "(",I,",",J,",",K,",",L,") = ", RIEM(I,J,K,L)$
  
% Ricci tensor
write "ricci tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", RICCI(I,J)$

% curvature scalar
write "curvature scalar = ",R$

% Einstein tensor
write "einstein tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",EINST(I,J)$

% place U
write "place u = ";
FOR I:=0:n-1 DO WRITE "(",I,") = ", U(I)$

% covariant derivative of U
write "covariant derivative of u = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",Ukv(I,J)$

% velocity V
write "velocity v = ";
FOR I:=0:n-1 DO WRITE "(",I,") = ", V(I)$

% local velocity wrt (x0,x1)
write "local velocity wrt (x0,x1) = ";
FOR I:=0:n-1 DO WRITE "(",I,") = ", LV(I)$

% acceleration
write "acceleration = ";
FOR I:=0:n-1 DO WRITE "(",I,") = ", B(I)$

% local acceleration wrt (x0,x1)
write "local acceleration wrt (x0,x1) = ";
FOR I:=0:n-1 DO WRITE "(",I,") = ", LB(I)$

% equation of motion
write "equation of motion = ";
FOR I:=0:n-1 DO WRITE "(",I,") = ", BG(I)$

% equation of motion 
on factor;
write "equation of motion = ";
FOR I:=0:n-1 DO WRITE "(",I,")  =", BG(I)$
off factor;

SHUT "metric2d_results.txt";

off revpri;
on nat;

END;

% calculation of metric tensor (3-dim space-time)

off echo;
on revpri;

n:=3;

operator x$
x(0):=t; x(1):=lambda0; x(2):=lambda1;

% rules
trig1:={sin(~x)^2=>(1-cos(x)^2)}$ let trig1$

% procedures
procedure scalprod(a,b); begin
integer n; n:=first(length(a))-1; 
result:=for i:=0:n-1 sum a(i)*b(i);
return result end;

procedure showmatrix(mm);begin integer m,n;l:=length(mm);m:=first(l)-1;n:=second(l)-1;matrix hm(m,n);for i:=0:m-1 do for j:=0:n-1 do hm(i+1,j+1):=mm(i,j); write hm end;

procedure showvector(vv);begin scalar n;n:=first(length(vv))-1;matrix hv(n,1);for i:=0:n-1 do hv(i+1,1):=vv(i);write hv end;

array f(n+1), dfdt(n+1), dfdl0(n+1), dfdl1(n+1)$

% current radius
a:=a0*sqrt(1-t^2);

% surface of hyper sphere in t and lambda

f(0):=a*cos(lambda0)*cos(lambda1);
f(1):=a*cos(lambda0)*sin(lambda1);
f(2):=a*sin(lambda0);
f(3):=a0*t;

for i:=0:n do dfdt(i):=df(f(i),x(0))$
for i:=0:n do dfdl0(i):=df(f(i),x(1))$
for i:=0:n do dfdl1(i):=df(f(i),x(2))$

array g(n,n)$

g(0,0):=scalprod(dfdt,dfdt)$
g(0,1):=scalprod(dfdt,dfdl0)$
g(0,2):=scalprod(dfdt,dfdl1)$

g(1,0):=scalprod(dfdl0,dfdt)$
g(1,1):=scalprod(dfdl0,dfdl0)$
g(1,2):=scalprod(dfdl0,dfdl1)$

g(2,0):=scalprod(dfdl1,dfdt)$
g(2,1):=scalprod(dfdl1,dfdl0)$
g(2,2):=scalprod(dfdl1,dfdl1)$

write "f = "; showvector(f);
write "df/dt = "; showvector(dfdt);
write "g = "; showmatrix(g);

off revpri;
on echo;
end;

% Calculations concerning the special metric of Dieter Egger
% Small and capital letters are treated as being equivalent

% Dimension of space-time
n:=3;

% turn off extra echoes
off echo;

% smaller exponents first
on revpri;

% Coordinates
OPERATOR X$
X(0):=t$
X(1):=lambda0$
X(2):=lambda1$

% Vectors (1-dim arrays start with index 0)
ARRAY U(n), V(n)$

% place (fixed to origin)
U(0):=a0*asin(t)$
U(1):=0$
U(2):=0$

% Rule
trig1:={sin(~x)^2=>(1-cos(x)^2)}$
let trig1$

% Procedure
procedure showMatrix(mm); begin
MATRIX hh(n,n)$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO hh(I+1,J+1):=mm(I,J)$
write hh; end;

% Arrays (2-dim arrays start with indices (0,0))
ARRAY G(n,n), GINV(n,n), CHRIST(n,n,n), RIEM(n,n,n,n), RICCI(n,n), EINST(n,n)$
ARRAY EIT(n,n), ENI(n,n)$

% Metric (cellar indices)
G(0,0):=a0^2/(1-t^2)$
G(1,1):=a0^2*(1-t^2)$
G(2,2):=a0^2*(1-t^2)*cos(lambda0)^2$

% Inverse Metric (roof indices)
MATRIX MG(n,n), MGINV(n,n)$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO MG(I+1,J+1):=G(I,J)$
MGINV:=1/MG$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO GINV(I,J):=MGINV(I+1,J+1)$

write "g = ",mg;
write "ginv = ",mginv;
write "g*ginv = ",mg*mginv;

% velocity
for i:=0:n-1 do v(i):=df(u(i),t)$

% max. velocity
Array vmax(n)$
svmax:=a0/sqrt(1-t^2)$
for i:=0:n-1 do vmax(i):=svmax$
svmaxq:=svmax*svmax$
write "max. velocity = ",svmax;

% energy impulse tensor (eit, roof indices)
for i:=0:n-1 do for j:=0:n-1 do eit(i,j):=v(i)*v(j)*(p/svmaxq + rho) - p * ginv(i,j)$  
write "eit roof = "; showMatrix(eit);

% energy impulse tensor (eni, cellar indices, including kappa)
for i:=0:n-1 do for j:=0:n-1 do eni(i,j) := - kappa * for k:=0:n-1 sum g(i,k)* for l:=0:n-1 sum g(j,l)*eit(k,l)$
write "eni = -kappa*(eit cellar) = "; showMatrix(eni);

% Christoffel symbols (Fliessbach)
for k:=0:n-1 do for l:=0:n-1 do for m:=0:n-1 do CHRIST(k,l,m):= for n:=0:n-1 sum GINV(k,n)/2 * (DF(G(m,n),X(l)) + DF(G(l,n),X(m)) - DF(G(m,l),X(n)))$
  
% curvature tensor (Fliessbach)
for m:=0:n-1 do for i:=0:n-1 do for k:=0:n-1 do for p:=0:n-1 do RIEM(m,i,k,p) :=  DF(CHRIST(m,i,k),X(p)) - DF(CHRIST(m,i,p),X(k)) + FOR r:=0:n-1 SUM CHRIST(r,i,k)*CHRIST(m,r,p) - CHRIST(r,i,p)*CHRIST(m,r,k)$
 
% Ricci tensor (Fliessbach)
FOR I:=0:n-1 DO FOR J:=0:n-1 DO RICCI(I,J):= FOR M:=0:n-1 SUM RIEM(M,I,M,J)$
write "ricci = "; showMatrix(ricci);

% curvature scalar
R:= FOR I:=0:n-1 SUM FOR J:=0:n-1 SUM GINV(I,J)*RICCI(I,J)$
write "curvature scalar r = ",r;

% Einstein tensor
FOR I:=0:n-1 DO FOR J:=0:n-1 DO EINST(I,J):=RICCI(I,J)-R/2*G(I,J)$
write "einstein = "; showMatrix(einst);

% solving field equations
write "solving field equations ...";
on factor;
erho:=solve(eni(0,0)=einst(0,0),rho)$
write "mass density = ", erho;

ep:=solve(eni(1,1)=einst(1,1),p)$
write "pressure = ", ep;
off factor;

ferho:=sub(ep,erho)$
write "final mass density = ", ferho;

%--------------------------------------------------------------
% write results to file
OUT "metric3d_results.txt";
off echo;
off nat;

% Metric
write "metric = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", G(I,J)$

% Inverse Metric
WRITE "inverse metric = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", GINV(I,J)$

% Christoffel symbols
write "christoffel symbols = ";
FOR K:=0:n-1 DO FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",K,",",I,",",J,") = ", CHRIST(K,I,J)$

% curvature tensor
write "curvature tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO FOR K:=0:n-1 DO FOR L:=0:n-1 DO WRITE "(",I,",",J,",",K,",",L,") = ", RIEM(I,J,K,L)$
  
% Ricci tensor
write "ricci tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", RICCI(I,J)$

% curvature scalar
write "curvature scalar = ",R$

% Einstein tensor
write "einstein tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",EINST(I,J)$

% energy impulse tensor
write "energy impulse tensor eit (roof) = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",eit(I,J);

% energy impulse tensor
write "energy impulse tensor eni (with kappa, cellar) = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",eni(I,J);

% solving field equations
write "solving field equations ...";
on factor;
write "mass density = ", erho;
write "pressure = ", ep;
off factor;
write "final mass density = ", ferho;

SHUT "metric3d_results.txt";

off revpri;
on nat;

END;

% calculation of metric tensor (4-dim space-time)

off echo;
on revpri;
n:=4;

array f(n+1), dfdt(n+1), dfdl0(n+1), dfdl1(n+1), dfdl2(n+1)$

array g(n,n)$

operator x$
x(0):=t; x(1):=lambda0; x(2):=lambda1; x(3):=lambda2;

% rules
trig1:={sin(~x)^2=>(1-cos(x)^2)}$ let trig1$

% procedures
procedure scalprod(a,b); begin
integer n; n:=first(length(a))-1; 
result:=for i:=0:n-1 sum a(i)*b(i);
return result end;

procedure showmatrix(mm);begin integer m,n;l:=length(mm);m:=first(l)-1;n:=second(l)-1;matrix hm(m,n);for i:=0:m-1 do for j:=0:n-1 do hm(i+1,j+1):=mm(i,j); write hm end;

procedure showvector(vv);begin integer n;n:=first(length(vv))-1; matrix hv(n,1); for i:=0:n-1 do hv(i+1,1):=vv(i); write hv end;

% current radius
a:=a0*sqrt(1-t^2);

% surface of hyper sphere in t and lambda

f(0):=a*cos(lambda0)*cos(lambda1)*cos(lambda2);
f(1):=a*cos(lambda0)*cos(lambda1)*sin(lambda2);
f(2):=a*cos(lambda0)*sin(lambda1);
f(3):=a*sin(lambda0);
f(4):=a0*t;

for i:=0:n do dfdt(i):=df(f(i),x(0))$
for i:=0:n do dfdl0(i):=df(f(i),x(1))$
for i:=0:n do dfdl1(i):=df(f(i),x(2))$
for i:=0:n do dfdl2(i):=df(f(i),x(3))$

g(0,0):=scalprod(dfdt,dfdt)$
g(0,1):=scalprod(dfdt,dfdl0)$
g(0,2):=scalprod(dfdt,dfdl1)$
g(0,3):=scalprod(dfdt,dfdl2)$

g(1,0):=scalprod(dfdl0,dfdt)$
g(1,1):=scalprod(dfdl0,dfdl0)$
g(1,2):=scalprod(dfdl0,dfdl1)$
g(1,3):=scalprod(dfdl0,dfdl2)$

g(2,0):=scalprod(dfdl1,dfdt)$
g(2,1):=scalprod(dfdl1,dfdl0)$
g(2,2):=scalprod(dfdl1,dfdl1)$
g(1,3):=scalprod(dfdl1,dfdl2)$

g(3,0):=scalprod(dfdl2,dfdt)$
g(3,1):=scalprod(dfdl2,dfdl0)$
g(3,2):=scalprod(dfdl2,dfdl1)$
g(3,3):=scalprod(dfdl2,dfdl2)$

write "f = "; showvector(f);
write "df/dt = "; showvector(dfdt);
write "g = "; showmatrix(g);

off revpri;
on echo;
end;

% Calculations concerning the special metric of Dieter Egger
% Small and capital letters are treated as being equivalent

% Dimension of space-time
n:=4;

% turn off extra echoes
off echo;

% smaller exponents first
on revpri;

% Coordinates
OPERATOR X$
X(0):=t$
X(1):=lambda0$
X(2):=lambda1$
X(3):=lambda2$

% Vectors (1-dim arrays start with index 0)
ARRAY U(n), V(n)$

% place (fixed to origin)
U(0):=a0*asin(t)$
U(1):=0$
U(2):=0$
U(3):=0$

% Rule
trig1:={sin(~x)^2=>(1-cos(x)^2)}$
let trig1$

% Procedure
procedure showMatrix(mm); begin
MATRIX hh(n,n)$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO hh(I+1,J+1):=mm(I,J)$
write hh; end;

% Arrays (2-dim arrays start with indices (0,0))
ARRAY G(n,n), GINV(n,n), CHRIST(n,n,n), RIEM(n,n,n,n), RICCI(n,n), EINST(n,n)$
ARRAY EIT(n,n), ENI(n,n)$

% Metric (cellar indices)
G(0,0):=a0^2/(1-t^2)$
G(1,1):=a0^2*(1-t^2)$
G(2,2):=a0^2*(1-t^2)*cos(lambda0)^2$
G(3,3):=a0^2*(1-t^2)*cos(lambda0)^2*cos(lambda1)^2$

% Inverse Metric (roof indices)
MATRIX MG(n,n), MGINV(n,n)$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO MG(I+1,J+1):=G(I,J)$
MGINV:=1/MG$
FOR I:=0:n-1 DO FOR J:=0:n-1 DO GINV(I,J):=MGINV(I+1,J+1)$

write "g = ",mg$
write "ginv = ",mginv$
write "g*ginv = ",mg*mginv$

% velocity
for i:=0:n-1 do v(i):=df(u(i),t)$

% max. velocity
Array vmax(n)$
svmax:=a0/sqrt(1-t^2)$
for i:=0:n-1 do vmax(i):=svmax$
svmaxq:=svmax*svmax$
write "max. velocity = ",svmax$

% energy impulse tensor (eit, roof indices)
for i:=0:n-1 do for j:=0:n-1 do eit(i,j):=v(i)*v(j)*(p/svmaxq + rho) - p * ginv(i,j)$  
write "eit roof = "$ showMatrix(eit)$

% energy impulse tensor (eni, cellar indices, including kappa)
for i:=0:n-1 do for j:=0:n-1 do eni(i,j) := - kappa * for k:=0:n-1 sum g(i,k)* for l:=0:n-1 sum g(j,l)*eit(k,l)$
write "eni = -kappa*(eit cellar) = "$ showMatrix(eni)$

% Christoffel symbols (Fliessbach)
for k:=0:n-1 do for l:=0:n-1 do for m:=0:n-1 do CHRIST(k,l,m):= for n:=0:n-1 sum GINV(k,n)/2 * (DF(G(m,n),X(l)) + DF(G(l,n),X(m)) - DF(G(m,l),X(n)))$
  
% curvature tensor (Fliessbach)
for m:=0:n-1 do for i:=0:n-1 do for k:=0:n-1 do for p:=0:n-1 do RIEM(m,i,k,p) :=  DF(CHRIST(m,i,k),X(p)) - DF(CHRIST(m,i,p),X(k)) + FOR r:=0:n-1 SUM CHRIST(r,i,k)*CHRIST(m,r,p) - CHRIST(r,i,p)*CHRIST(m,r,k)$
 
% Ricci tensor (Fliessbach)
FOR I:=0:n-1 DO FOR J:=0:n-1 DO RICCI(I,J):= FOR M:=0:n-1 SUM RIEM(M,I,M,J)$
write "ricci = "$ showMatrix(ricci)$

% curvature scalar
R:= FOR I:=0:n-1 SUM FOR J:=0:n-1 SUM GINV(I,J)*RICCI(I,J)$
write "curvature scalar r = ",r;

% Einstein tensor
FOR I:=0:n-1 DO FOR J:=0:n-1 DO EINST(I,J):=RICCI(I,J)-R/2*G(I,J)$
write "einstein = "$ showMatrix(einst)$

% solving field equations
write "solving field equations ...";
on factor;
erho:=solve(eni(0,0)=einst(0,0),rho)$
write "mass density = ", erho$

ep:=solve(eni(1,1)=einst(1,1),p)$
write "pressure = ", ep$
off factor;

ferho:=sub(ep,erho)$
write "final mass density = ", ferho$

%--------------------------------------------------------------
% write results to file
OUT "metric4d_results.txt";
off echo;
off nat;

% Metric
write "metric = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", G(I,J)$

% Inverse Metric
WRITE "inverse metric = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", GINV(I,J)$

% Christoffel symbols
write "christoffel symbols = ";
FOR K:=0:n-1 DO FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",K,",",I,",",J,") = ", CHRIST(K,I,J)$

% curvature tensor
write "curvature tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO FOR K:=0:n-1 DO FOR L:=0:n-1 DO WRITE "(",I,",",J,",",K,",",L,") = ", RIEM(I,J,K,L)$
  
% Ricci tensor
write "ricci tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ", RICCI(I,J)$

% curvature scalar
write "curvature scalar = ",R$

% Einstein tensor
write "einstein tensor = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",EINST(I,J)$

% energy impulse tensor
write "energy impulse tensor eit (roof) = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",eit(I,J);

% energy impulse tensor
write "energy impulse tensor eni (with kappa, cellar) = ";
FOR I:=0:n-1 DO FOR J:=0:n-1 DO WRITE "(",I,",",J,") = ",eni(I,J);

% solving field equations
write "solving field equations ...";
on factor;
write "mass density = ", erho;
write "pressure = ", ep;
off factor;
write "final mass density = ", ferho;

SHUT "metric4d_results.txt";

off revpri;
on nat;

END;

% dynamical identifiers
for i:=1:5 do set (mkid (z, i), i);
for i:=1:5 do write mkid (z, i);
for i:=1:5 sum mkid (z, i);
end;

% partial fractions
b:=(x+1)^2*(x+2);
a:=2/b;
pf (a, x);
c:=ws;
first (c)+second (c)+third (c);
off exp;
ws;
first (c)+second (c)+third (c);
end;

% define full path to desired Reduce input file
in "/storage/extSdCard/Reduce/constants.red";
lambda0:=5*10^(-6);
nu:=c/lambda0;
lambda0:=5*10^(-7);
nu:=c/lambda0;
h*nu/c;
h*nu;
h*nu/el;
end;

abs (-3/4);
abs (2a);
abs (i);
abs (-x);
ceiling (-5/4);
ceiling (-3/4);
ceiling (3/4);
ceiling (5/4);
ceiling (-a);
conj (1+i);
conj (a+i*b);
factorial 5;
factorial a;
fix (-5/4);
fix (-5/4);
fix (-3/4);
fix (3/4);
fix (5/4);
floor (-5/4);
floor (-3/4);
floor (3/4);
floor (5/4);
floor (a);
impart (1+i);
impart (1+2i);
impart (a+i*b);
max (2,-3, 4, 5);
min (2,-3, 4, 5);
min (2,-3, 4, 5, a);
max (2,-3, 4, 5, a);
nextprime 1111;
nextprime 0;
nextprime 11;
nextprime 111;
random 49;
random 49;
random 49;
random 49;
random 49;
random 49;
random_new_seed (1000);
random 49;
random 49;
random_new_seed (1000);
random 49;
random 49;
repart (4+2i);
repart (4a+2i);
repart (a+i*b);
round (-5/4);
round (-3/4);
round (3/4);
round (5/4);
round (7/4);
sign (-5);
sign (-5/4);
sign (5/4);
sign (0);
end;

%%%%%%%%%%%%%%%%%%%%%
%  PROGRAMMING
%%%%%%%%%%%%%%%%%%%%%

% Define number to factorize
x:=42;

% Factorize x and write out each individual 
% factor
factors:=factorize(fix(x))$
x:=0$
for i:=1:length(factors) do begin
   q:=part(factors,i);
   for j:=1:part(q,2) do begin
      x:=x+1;
      write "factor ", x, ": ", part(q,1);
   end;
end;

% Procedure to calculate Legendre polynomial
% using recursion 
procedure p(n,x);
   if n<0 then rederr "Invalid argument to p(n,x)"
   else if n=0 then 1
   else if n=1 then x
   else ((2*n-1)*x*p(n-1,x)-(n-1)*p(n-2,x))/n$

% Enable fancy output
%fancy_output$

% Calculate p(2,w)
write "P(2,w) = ", p(2,w);

% Incidentally, p(n,x) can be calculated more
% efficiently as follows
procedure p(n,x);
   sub(y=0,df(1/(y^2-2*x*y+1)^(1/2),y,n))/(for i:=1:n product i)$

write "P(3,w) = ", p(3,w);

end;

% Rules
%trig1:={sin(~x)^2+cos(~x)^2=>1, cos(~x)^2+sin(~x)^2=>1};
trig1:={sin(~x)^2=>(1-cos(x)^2)};
let trig1;
trig2:={tan (~x)=>(sin (x)/cos (x))};
let trig2;
end;

procedure scalprod(a,b); begin;
scalar n; 
n:=first (length(a))-1;
result:=for i:=0:n sum a(i)*b(i);
return result;
end;
end;

c:=1/sqrt (epsilon0*mu0);
end;