Overview
Comment:Scripts by Dieter (Olli) Egger (November 30th 2019)
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | descendants | 2009.11.30 | trunk
Files: files | file ages | folders
SHA3-256: 5ce5316a33430f7c474932758c3e40015fc8d5ae521ac6b2a0abdc4de0128e9e
User & Date: jeff@gridfinity.com on 2021-03-01 07:23:02
Other Links: manifest | tags
Context
2021-03-01
07:27:08
Update README.md: Small update, markdownlint fixes check-in: bc03f92dc3 user: jeff@gridfinity.com tags: master, trunk
07:23:02
Scripts by Dieter (Olli) Egger (November 30th 2019) Leaf check-in: 5ce5316a33 user: jeff@gridfinity.com tags: 2009.11.30, trunk
Changes

Added .gitattributes version [f9d1f1840e].













>
>
>
>
>
>
1
2
3
4
5
6
# Effectively disable Linguist
doc/* linguist-vendored
src/* linguist-vendored
*.md linguist-vendored
*.txt linguist-vendored
*.red linguist-vendored

Added .github/issue-close-app.yml version [9a3b66da03].

















>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
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"

Added .github/mergeable.yml version [7fb133d986].























>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
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

Added .gitignore version [a7ffc6f8bf].

Added README.md version [a1a82c2af5].















































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# 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)

Added SECURITY.md version [4bf8ab023a].











>
>
>
>
>
1
2
3
4
5
# Security Information

## Reporting a Vulnerability

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

Added doc/MetricUniverse.pdf version [bdfb9288e6].

cannot compute difference between binary files

Added doc/Metrik.pdf version [2a57d66e37].

cannot compute difference between binary files

Added src/algebra.txt version [4a09a9e43d].



















































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
%%%%%%%%%%%%%%%%%%%%%
%  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;

Added src/analyze.red version [5edece8dce].











































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
% 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;

Added src/boolean.txt version [62850ba20a].





















































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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";

Added src/calculus.txt version [4984952ae5].









































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
%%%%%%%%%%%%%%%%%%%%%
%  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;

Added src/clear_constants.red version [366dfd4dfc].

































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
% 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;

Added src/constants.red version [3bd159fe32].

































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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;

Added src/example01.red version [26599d6c9f].

















>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
a:=sin (x);
b:=cos (x);
c:=a^2+b^2;
d:=a/b;
f:=tan (x)-d;
trigsimp c;
trigsimp f;
end;

Added src/example02.red version [db6f1a63d4].



























>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
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;

Added src/example03.red version [f5853dc345].









































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
% 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;

Added src/example04.red version [0addcc0385].



















>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
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;

Added src/example05.red version [f430f73384].















































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
(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;

Added src/heron.txt version [3e5894bb3b].















>
>
>
>
>
>
>
1
2
3
4
5
6
7
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;

Added src/hypermeijerg.txt version [1bd35ffc38].

























































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%%%%%%%%%%%%%%%%%%%%%
%  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;

Added src/integLimit.red version [849ec37768].





















>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
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;

Added src/integral.txt version [9ab6dc4c0e].























>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
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;

Added src/introReduce.txt version [938a1ca51e].









































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
% 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;
%

Added src/linalg.txt version [e9c1c226ce].





































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
%%%%%%%%%%%%%%%%%%%%%
%  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;

Added src/mat.red version [24a6c2b757].















>
>
>
>
>
>
>
1
2
3
4
5
6
7
% 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;

Added src/metric2calc.red version [027137da09].





























































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
% 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;

Added src/metric2d.red version [041b81ff2e].







































































































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
% 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;

Added src/metric3calc.red version [3b9fd98e09].

























































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
% 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;

Added src/metric3d.red version [87125fc9f6].

































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
% 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;

Added src/metric4calc.red version [12312f33f4].











































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
% 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;

Added src/metric4d.red version [105692f385].







































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
% 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;

Added src/mkid.txt version [6aa00a225d].











>
>
>
>
>
1
2
3
4
5
% 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;

Added src/partialFraction.txt version [f1283d1b2c].





















>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
% 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;

Added src/photonenergy.txt version [4d32c62a43].





















>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
% 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;

Added src/prefix.txt version [9a0ec2ab9b].























































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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;

Added src/programming.txt version [e35a9918fc].



















































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
%%%%%%%%%%%%%%%%%%%%%
%  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;

Added src/rules.red version [e266788fa7].















>
>
>
>
>
>
>
1
2
3
4
5
6
7
% 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;

Added src/scalprod.red version [085ebbfed8].















>
>
>
>
>
>
>
1
2
3
4
5
6
7
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;

Added src/speedoflight.txt version [8a84a36ba8].





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


olli-scripts
English Homepage | German Homepage | DL2MIE | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]