| blob | 79e111b044e96650ea184a6888cb27afc1ae753b |
1 /*-*MACSYMA-*-*/
2 /*
3 * This program is free software; you can redistribute it and/or
4 * modify it under the terms of the GNU General Public License as
5 * published by the Free Software Foundation; either version 2 of
6 * the License, or (at your option) any later version.
7 *
8 * This program is distributed in the hope that it will be
9 * useful, but WITHOUT ANY WARRANTY; without even the implied
10 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 * PURPOSE. See the GNU General Public License for more details.
12 *
13 * Comments:
14 *
15 * The ctensor package was downcased, cleaned up, and moving frames
16 * functionality was added by Viktor Toth (http://www.vttoth.com/).
17 *
18 * As of November, 2004, the naming conventions in this package now
19 * correspond with the naming conventions in commercial MACSYMA.
20 */
22 if get('ctensor,'version)#false then error("CTENSOR already loaded!")$
24 /***********************************
25 PACKAGE INITIALIZATION
26 ***********************************/
27 tensorkill:true$
28 define_variable(dim,4,fixnum,"The dimension of the problem.")$
29 define_variable(diagmetric,false,boolean,
30 "True if the metric entered is a diagonal matrix.")$
31 define_variable(ctrgsimp,false,boolean,
32 "Apply trigonometric simplifications (TRIGSIMP).")$
33 define_variable(cframe_flag,false,boolean,
34 "Perform computations in a moving frame.")$
35 define_variable(ctorsion_flag,false,boolean,"Apply torsion")$
36 define_variable(cnonmet_flag,false,boolean,"Apply nonmetricity")$
37 /* define_variable(derivabbrev,true,boolean)$ */
38 define_variable(ctayswitch,false,boolean)$
39 define_variable(ratchristof,true,boolean)$
40 define_variable(rateinstein,true,boolean)$
41 define_variable(ratriemann,true,boolean)$
42 define_variable(ratweyl,true,boolean)$
44 init_ctensor():=
45 (
46 tensorkill:true,
47 dim:4,
48 diagmetric:false,
49 setflags(),
50 kill(ct_coords,lg,ug,lcs,mcs,geod,ric,uric,tracer,ein,lein,
51 riem,lriem,uriem,kinvariant,lfg,ufg,fr,fri,fb,tr,kt,nm,nmc),
52 lfg:diagmatrix(dim,1),
53 lfg[1,1]:-1,
54 done
55 )$
57 _ug():=if cframe_flag then 'ufg else 'ug;
58 _lg():=if cframe_flag then 'lfg else 'lg;
60 /***********************************
61 HELPER FUNCTIONS
62 ***********************************/
63 setflags():=(derivabbrev:true,ctayswitch:false,
64 ratchristof:true,rateinstein:true,ratriemann:true,ratweyl:true)$
66 readvalue(message,pred,[badboy]):=
67 block([value],
68 loop,
69 value:read(message),
70 if apply(pred,[value])=true then return(value),
71 if badboy#[] then apply('print,badboy),
72 go(loop))$
74 modedeclare(function(yesp),boolean)$
75 yesp([messages]):=block
76 (
77 [reply],
78 loop,
79 reply:apply('read,messages),
80 if member(reply,['yes,'ye,'y,'yes,'yes,'y]) then return(true),
81 if member(reply,['no,'n,'no,'no,'n]) then return(false),
82 print("Please reply YES or NO."),
83 go(loop)
84 )$
86 resimp(exp):=apply('ev,[exp,'noeval,'simp])$
88 kdelt(i,j):=(modedeclare([i,j],fixnum),if i = j then 1 else 0)$
90 /* Optional Taylor-series approximation */
91 ctaylor(x):=
92 if ctayswitch then ratdisrep(taylor(x,ctayvar,ctaypt,ctaypov-1)) else x$
94 /* This function transforms second-rank covariant tensors given the components
95 and the transformation law, which is in the form Yi=Yi(X1,...,Xdim) */
96 ctransform(qxyz,[equ]):=block
97 (
98 [ttemp,omtemp,vartemp,newtemp],
99 if length(equ)=0 then
100 (
101 for i thru dim do vartemp[i]:read("Transform #",i)
102 )
103 else
104 (
105 if length(equ)=1 then
106 (
107 if listp(equ[1])#true then print("argument 2 needs to be list")
108 else for i thru dim do vartemp[i]:equ[1][i]
109 )
110 else print("too many arguments")
111 ),
112 for i thru dim do
113 for j thru dim do
114 (
115 ttemp[i,j]:qxyz[i,j],
116 for k thru dim do ttemp[i,j]:subst(omtemp[k],ct_coords[k],ttemp[i,j])
117 ),
118 for i thru dim do
119 for j thru dim do
120 for k thru dim do ttemp[i,j]:subst(vartemp[k],omtemp[k],ttemp[i,j]),
121 for a thru dim do
122 for b thru dim do newtemp[a,b]:txyzsum(vartemp,ct_coords,a,b,ttemp),
123 if ctrgsimp=true then ratsimp(trigsimp(genmatrix(newtemp,dim,dim)))
124 else genmatrix(newtemp,dim,dim)
125 )$
127 txyzsum(vartemp,ct_coords,a,b,ttemp):=block(
128 [temp:0],
129 for i from 1 thru dim do
130 for j from 1 thru dim do
131 temp:temp+diff(vartemp[i],ct_coords[a])*diff(vartemp[j],ct_coords[b])*
132 ttemp[i,j],
133 temp
134 )$
136 /* Check diagonality of a matrix or (2D) array */
137 diagmatrixp(mat,nlim):=
138 (
139 modedeclare(nlim,fixnum),
140 block
141 (
142 [diagflag:true],
143 for i thru nlim do
144 for j thru nlim do
145 if i#j and mat[i,j]#0 then return(diagflag:false),
146 diagflag
147 )
148 )$
150 /* Check symmetry of a matrix or (2D) array */
151 symmetricp(m,n):=
152 (
153 modedeclare(n,fixnum),
154 block
155 (
156 [symflag:true],
157 if n#1 then for i thru n-1 do
158 for j from i+1 thru n do
159 if m[j,i]#m[i,j] then return(symflag:false),
160 symflag
161 )
162 )$
164 /* Count the number of terms in a (2D) array */
165 ntermst(f):=for i thru dim do for j thru dim do print([[i,j],nterms(f[i,j])])$
167 /* Display the contents of a tensor */
168 cdisplay(_t,[_k]):=block
169 (
170 [_len,_tmp,_i,_j,_n],
171 if member(_t,arrays) then
172 (
173 _tmp : funmake(arrayinfo, [_t]),
174 _len : ev(_tmp)[2] - length(_k),
175 if _len > 2 then
176 (
177 for _n thru dim do apply(cdisplay,append([_t],append(_k,[_n])))
178 )
179 else if _len = 2 then block
180 (
181 [_tmp:zeromatrix(dim,dim)],
182 for _i thru dim do for _j thru dim do
183 _tmp[_i,_j]:arrayapply(_t,append(_k,[_i,_j])),
184 disp((if length(_k) > 0 then arrayapply(nounify(_t),_k)
185 else nounify(_t))=_tmp)
186 )
187 else if _len = 1 then (true)
188 else print("Array is empty")
189 )
190 else disp(_t=ev(_t)),
191 done
192 )$
194 /* Delete an element */
195 deleten(l,n):=
196 (
197 modedeclare(n,fixnum),
198 block
199 (
200 [len],
201 modedeclare(len,fixnum),
202 if l = [] then l
203 else
204 (
205 len:length(l),
206 if n>len or n<1 then error("Second argument out of range")
207 else if n=1 then rest(l)
208 else if n=len then rest(l,-1)
209 else append(rest(l,n-len-1),rest(l,n))
210 )
211 )
212 )$
214 /***********************************
215 SETTING UP THE METRIC
216 ***********************************/
217 csetup():=block
218 (
219 [],
220 if not tensorkill then
221 error("Type KILL(ALL)\; and then TENSORKILL:TRUE\;
222 before you enter a new metric."),
223 tensorkill:false,
224 setflags(),
225 dim:mode_identity(fixnum,
226 readvalue("Enter the dimension of the coordinate system:",
227 lambda([v], if integerp(v) then block
228 (
229 [u:v],
230 modedeclare(u,fixnum),
231 if u>0 then true
232 )
233 ),
234 "Must be a positive integer!")
235 ),
236 ct_coords:if dim=2 then [x,y]
237 else if dim = 3 then [x,y,z]
238 else if dim = 4 then [x,y,z,t]
239 else read("Enter a list containing the names of the coordinates in order"),
240 if yesp("Do you wish to change the coordinate names?") then
241 ct_coords:
242 read("Enter a list containing the names of the coordinates in order"),
243 if length(ct_coords)#dim then
244 error("Length of the coordinate list is not equal to the dimension"),
246 readvalue("Do you want to
247 1. Enter a new metric?
248 2. Enter a metric from a file?
249 3. Approximate a metric with a Taylor series?",
250 lambda([opt], if opt = 1 then(newmet(),true)
251 else if opt = 2 then(filemet(),true)
252 else if opt = 3 then(sermet(),true)
253 else false
254 ),
255 "Invalid option, please enter 1, 2, or 3."
256 ),
257 done
258 )$
260 /* Enter a new metric by hand */
261 newmet():=block
262 (
263 [],
264 lg:entermatrix(dim,dim),
265 read("Enter functional dependencies with DEPENDS or 'N' if none"),
266 if yesp("Do you wish to see the metric?") then disp(lg),
267 cmetric()
268 )$
270 /* Read a metric from a file */
271 filemet():=block
272 (
273 [file,fpos],
274 file:read("Specify the file as [filename1,filename2,directory]?"),
275 fpos:(print("What is the ordinal position of the metric in this file?"),
276 read("(Note, the name LG must be assigned to your metric in the file)")),
277 apply(batch,[file,off,fpos,fpos]),cmetric()
278 )$
280 /* Approximate the metric with a Taylor series */
281 sermet():=block
282 (
283 [],
284 ctayswitch:true,
285 ctayvar:read("Enter expansion parameter"),
286 ctaypov:read("Enter minimum power to drop"),
287 ctaypt:read("Enter the point to expand the series around"),
288 if yesp("Is the metric in a file?") then filemet() else newmet()
289 )$
291 /* Computing the frame field from the inverse frame field (user-defined) */
292 tmetric(disp):=block(
293 [],
294 /* fr are the contravariant frame vector components. fri are covariant.
295 i.e.,
296 i
297 fr[a,i] = e , fri[a,i] = e
298 (a) (a)i
299 */
300 fr:lfg.(transpose(fri)^^(-1)),
301 ufg:lfg^^(-1),
302 lg:transpose(lfg.fri).fri,
304 if ctrgsimp then lg:trigsimp(lg)
305 else lg:ratsimp(lg),
306 if ratfac then lg:factor(lg),
307 if disp then ldisplay(lg)
308 );
310 /* Computing the contravariant metric or the frame field */
311 cmetric([dis]):=if cframe_flag then tmetric(if dis#[] then dis[1] else false)
312 else block
313 (
314 [],
315 if length(lg)#dim or length(transpose(lg))#dim then
316 error("The rank of the metric is not equal to the dimension of the space"),
317 /* else if not symmetricp(lg,dim) then
318 error("The metric tensor must be symmetric!"),*/
320 diagmetric:diagmatrixp(lg,dim),
321 gdet:factor(determinant(lg)),
322 ug:ident(length(first(lg))),
323 if diagmetric then for ii:1 thru length(first(lg)) do ug[ii,ii]:1/lg[ii,ii]
324 else ug:block([detout:true],transpose(lg^^(-1))),
325 if ratfac then ug:ratsimp(ug),
326 if dis#[] and dis[1]#false then
327 if yesp("Do you wish to see the metric inverse?") then ldisp(ug),
328 ug:resimp(ug),
329 done
330 )$
332 ct_coordsys(coordinate_system, [extra_args]):=
333 (
334 if coordinate_system='cartesian2d then
335 (
336 dim:2,
337 ct_coords:[x,y],
338 if cframe_flag then
339 (
340 fri:ident(2),
341 lfg:ident(2)
342 )
343 else lg:ident(2)
344 )
345 else if coordinate_system='polar then
346 (
347 dim:2,
348 ct_coords:[r,phi],
349 assume(r>=0),
350 if cframe_flag then
351 (
352 fri:matrix([cos(phi),-sin(phi)*r],
353 [sin(phi), cos(phi)*r]),
354 lfg:ident(2)
355 )
356 else lg:matrix([1,0],[0,r^2])
357 )
358 else if coordinate_system='elliptic then
359 (
360 dim:2,
361 ct_coords:[u,v],
362 declare(e,constant),
363 if cframe_flag then
364 (
365 fri:matrix([e*sinh(u)*cos(v),-e*cosh(u)*sin(v)],
366 [e*cosh(u)*sin(v), e*sinh(u)*cos(v)]),
367 lfg:ident(2)
368 )
369 else lg:matrix([e^2*(cosh(u)^2-cos(v)^2), 0],
370 [ 0, e^2*(cosh(u)^2-cos(v)^2)])
371 )
372 else if coordinate_system='confocalelliptic then
373 block(
374 [tmp:sqrt((u^2-1)*(1-v^2))],
375 dim:2,
376 ct_coords:[u,v],
377 declare(e,constant),
378 if cframe_flag then
379 (
380 fri:matrix([ e*v, e*u],
381 [e*u*(1-v^2)/tmp, e*v*(1-u^2)/tmp]),
382 lfg:ident(2)
383 )
384 else lg:matrix([e^2*(u^2-v^2)/(u^2-1), 0],
385 [0, e^2*(v^2-u^2)/(v^2-1)])
386 )
387 else if coordinate_system='bipolar then
388 block(
389 [tmp:(cosh(v)-cos(u))^2],
390 dim:2,
391 ct_coords:[u,v],
392 declare(e,constant),
393 if cframe_flag then
394 (
395 fri:matrix([ -e*sin(u)*sinh(v)/tmp, e*(1-cos(u)*cosh(v))/tmp],
396 [e*(cos(u)*cosh(v)-1)/tmp, -e*sin(u)*sinh(v)/tmp]),
397 lfg:ident(2)
398 )
399 else lg:matrix([e^2/tmp,0],[0,e^2/tmp])
400 )
401 else if coordinate_system='parabolic then
402 (
403 dim:2,
404 ct_coords:[u,v],
405 if cframe_flag then
406 (
407 fri:matrix([u,-v],[v,u]),
408 lfg:ident(2)
409 )
410 else lg:matrix([u^2+v^2, 0], [0, u^2+v^2])
411 )
412 else if coordinate_system='cartesian3d then
413 (
414 dim:3,
415 ct_coords:[x,y,z],
416 if cframe_flag then
417 (
418 fri:ident(3),
419 lfg:ident(3)
420 )
421 else lg:ident(3)
422 )
423 else if coordinate_system='polarcylindrical then
424 (
425 dim:3,
426 ct_coords:[r,theta,z],
427 if cframe_flag then
428 (
429 fri:matrix([cos(theta), -r*sin(theta), 0],
430 [sin(theta), r*cos(theta), 0],
431 [ 0, 0, 1]),
432 lfg:ident(3)
433 )
434 else lg:matrix([1, 0, 0], [0, r^2, 0], [0, 0, 1])
435 )
436 else if coordinate_system='ellipticcylindrical then
437 (
438 dim:3,
439 ct_coords:[u,v,z],
440 if cframe_flag then
441 (
442 fri:matrix([e*sinh(u)*cos(v), -e*cosh(u)*sin(v), 0],
443 [e*cosh(u)*sin(v), e*sinh(u)*cos(v), 0],
444 [ 0, 0, 1]),
445 lfg:ident(3)
446 )
447 else lg:matrix([e^2*(sin(v)^2+sinh(u)^2), 0, 0],
448 [ 0, e^2*(sin(v)^2+sinh(u)^2), 0],
449 [ 0, 0, 1])
450 )
451 else if coordinate_system='confocalellipsoidal then
452 (
453 dim:3,
454 ct_coords:[u,v,w],
455 declare([e,f,g],constant),
456 if cframe_flag then
457 (
458 fri:matrix([-sqrt((v-e^2)*(w-e^2)/(f^2-e^2)/(g^2-e^2)/(e^2-u))/2,
459 -sqrt((u-e^2)*(w-e^2)/(f^2-e^2)/(g^2-e^2)/(e^2-v))/2,
460 -sqrt((u-e^2)*(v-e^2)/(f^2-e^2)/(g^2-e^2)/(e^2-w))/2],
462 [-sqrt((v-f^2)*(w-f^2)/(f^2-e^2)/(g^2-f^2)/(u-f^2))/2,
463 -sqrt((u-f^2)*(w-f^2)/(f^2-e^2)/(g^2-f^2)/(v-f^2))/2,
464 sqrt((u-f^2)*(v-f^2)/(f^2-e^2)/(g^2-f^2)/(w-f^2))/2],
466 [-sqrt((v-g^2)*(w-g^2)/(g^2-e^2)/(g^2-f^2)/(g^2-u))/2,
467 -sqrt((u-g^2)*(w-g^2)/(g^2-e^2)/(g^2-f^2)/(g^2-v))/2,
468 -sqrt((u-g^2)*(v-g^2)/(g^2-e^2)/(g^2-f^2)/(g^2-w))/2]),
469 lfg:ident(3)
470 )
471 else lg:matrix([(v-u)*(w-u)/4/(e^2-u)/(u-f^2)/(u-g^2), 0, 0],
472 [0, (v-u)*(w-v)/4/(v-e^2)/(v-f^2)/(v-g^2), 0],
473 [0, 0, (w-u)*(w-v)/4/(e^2-w)/(w-f^2)/(w-g^2)])
474 )
475 else if coordinate_system='bipolarcylindrical then
476 block(
477 [tmp:(cosh(v)-cos(u))^2],
478 dim:3,
479 ct_coords:[u,v,z],
480 declare(e,constant),
481 if cframe_flag then
482 (
483 fri:matrix([e*sin(u)*sinh(v)/tmp, e*(1-cos(u)*cosh(v))/tmp, 0],
484 [e*(cos(u)*cosh(v)-1)/tmp, e*sin(u)*sinh(v)/tmp, 0],
485 [ 0, 0, 1]),
486 lfg:ident(3)
487 )
488 else lg:matrix([e^2/tmp, 0, 0],
489 [0, e^2/tmp, 0],
490 [0, 0, 1])
491 )
492 else if coordinate_system='paraboliccylindrical then
493 (
494 dim:3,
495 ct_coords:[u,v,z],
496 if cframe_flag then
497 (
498 fri:matrix([u,-v,0],[v,u,0],[0,0,1]),
499 lfg:ident(3)
500 )
501 else lg:matrix([u^2+v^2, 0, 0],[0, u^2+v^2, 0],[0, 0, 1])
502 )
503 else if coordinate_system='paraboloidal then
504 (
505 dim:3,
506 ct_coords:[phi,u,v],
507 if cframe_flag then
508 (
509 fri:matrix([ 0, u, -v],
510 [-sin(phi)*u*v, cos(phi)*v, cos(phi)*u],
511 [ cos(phi)*u*v, sin(phi)*v, sin(phi)*u]),
512 lfg:ident(3)
513 )
514 else lg:matrix([u^2*v^2, 0, 0],[0, u^2+v^2, 0], [0, 0, u^2+v^2])
515 )
516 else if coordinate_system='conical then
517 (
518 dim:3,
519 ct_coords:[u,v,w],
520 declare([e,f],constant),
521 if cframe_flag then
522 (
523 fri:matrix([v*w/e/f, u*w/e/f, u*v/e/f],
524 [-u*w/sqrt((f^2-e^2)*(u^2-e^2)/(e^2-v^2))/e,
525 v*w/sqrt((f^2-e^2)*(v^2-e^2)/(e^2-u^2))/e,
526 -sqrt((u^2-e^2)*(v^2-e^2)/(e^2-f^2))/e],
527 [ u*w/sqrt((f^2-e^2)*(u^2-f^2)/(v^2-f^2))/f,
528 v*w/sqrt((f^2-e^2)*(v^2-f^2)/(u^2-f^2))/f,
529 sqrt((u^2-f^2)*(v^2-f^2)/(f^2-e^2))/f]),
530 lfg:ident(3)
531 )
532 else lg:matrix([(v-u)*(v+u)*w^2/((u-e)*(u+e)*(u-f)*(u+f)), 0, 0],
533 [0, (u-v)*(u+v)*w^2/((v-e)*(v+e)*(v-f)*(v+f)), 0],
534 [0, 0, 1])
535 )
536 else if coordinate_system='toroidal then
537 (
538 dim:3,
539 ct_coords:[phi,u,v],
540 declare(e,constant),
541 if cframe_flag then
542 (
543 fri:matrix([-e*sin(phi)*sinh(v)/(cosh(v)-cos(u)),
544 -e*cos(phi)*sin(u)*sinh(v)/(cosh(v)-cos(u))^2,
545 -e*cos(phi)*(cos(u)*cosh(v)-1)/(cosh(v)-cos(u))^2],
546 [ e*cos(phi)*sinh(v)/(cosh(v)-cos(u)),
547 -e*sin(phi)*sin(u)*sinh(v)/(cosh(v)-cos(u))^2,
548 -e*sin(phi)*(cos(u)*cosh(v)-1)/(cosh(v)-cos(u))^2],
549 [ 0, e*(cos(u)*cosh(v)-1)/(cosh(v)-cos(u))^2,
550 -e*sin(u)*sinh(v)/(cosh(v)-cos(u))^2]),
551 lfg:ident(3)
552 )
553 else lg:matrix([e^2*sinh(v)^2/(cosh(v)-cos(u))^2, 0, 0],
554 [0, e^2/(cosh(v)-cos(u))^2, 0],
555 [0, 0, e^2/(cosh(v)-cos(u))^2])
556 )
557 else if coordinate_system='spherical then
558 (
559 dim:3,
560 ct_coords:[r,theta,phi],
561 assume(r>=0),
562 if cframe_flag then
563 (
564 fri:matrix([1,0,0],[0,r,0],[0,0,r*sin(theta)]),
565 lfg:ident(3)
566 )
567 else lg:matrix([1,0,0],[0,r^2,0],[0,0,r^2*sin(theta)^2])
568 )
569 else if coordinate_system='oblatespheroidal then
570 (
571 dim:3,
572 ct_coords:[u,v,phi],
573 declare(e,constant),
574 if cframe_flag then
575 (
576 fri:matrix([abs(e)*sqrt(sin(v)^2+sinh(u)^2), 0, 0],
577 [0, abs(e)*sqrt(sin(v)^2+sinh(u)^2), 0],
578 [0, 0, abs(e)*cosh(u)*abs(cos(v))]),
579 lfg:ident(3)
580 )
581 else lg:matrix([e^2*(sin(v)^2+sinh(u)^2), 0, 0],
582 [0, e^2*(sin(v)^2+sinh(u)^2), 0],
583 [0, 0, e^2*cosh(u)^2*cos(v)^2])
584 )
585 else if coordinate_system='oblatespheroidalsqrt then
586 (
587 dim:3,
588 ct_coords:[u,v,phi],
589 declare(e,constant),
590 if cframe_flag then
591 (
592 fri:matrix([abs(e)*sqrt((u^2-v^2)/(u^2-1)), 0, 0],
593 [0, abs(e)*sqrt((u^2-v^2)/(1-v^2)), 0],
594 [0, 0, abs(e*u*v)]),
595 lfg:ident(3)
596 )
597 else lg:matrix([e^2*(u^2-v^2)/(u^2-1), 0, 0],
598 [0, e^2*(u^2-v^2)/(1-v^2), 0],
599 [0, 0, e^2*u^2*v^2])
600 )
601 else if coordinate_system='prolatespheroidal then
602 (
603 dim:3,
604 ct_coords:[u,v,phi],
605 declare(e,constant),
606 if cframe_flag then
607 (
608 fri:matrix([abs(e)*sqrt(sin(v)^2+sinh(u)^2), 0, 0],
609 [0, abs(e)*sqrt(sin(v)^2+sinh(u)^2), 0],
610 [ 0, 0, abs(e*sinh(u)*sin(v))]),
611 lfg:ident(3)
612 )
613 else lg:matrix([e^2*(sin(v)^2+sinh(u)^2), 0, 0],
614 [0, e^2*(sin(v)^2+sinh(u)^2), 0],
615 [0, 0, e^2*sin(v)^2*sinh(u)^2])
616 )
617 else if coordinate_system='prolatespheroidalsqrt then
618 (
619 dim:3,
620 ct_coords:[u,v,phi],
621 declare(e,constant),
622 if cframe_flag then
623 (
624 fri:matrix([abs(e)*sqrt((v^2-u^2)/(1-u^2)), 0, 0],
625 [0, abs(e)*sqrt((v^2-u^2)/(v^2-1)), 0],
626 [0, 0, abs(e)*sqrt((1-u^2)*(v^2-1))]),
627 lfg:ident(3)
628 )
629 else lg:matrix([e^2*(v^2-u^2)/(1-u^2), 0, 0],
630 [0, e^2*(v^2-u^2)/(v^2-1), 0],
631 [0, 0, e^2*(1-u^2)*(v^2-1)])
632 )
633 else if coordinate_system='ellipsoidal then
634 (
635 dim:3,
636 ct_coords:[r,theta,phi],
637 declare([a,b,c],constant),
638 if cframe_flag then error("No frame fields yet")
639 else lg:matrix(
640 [(a^2*cos(phi)^2+b^2*sin(phi)^2)*sin(theta)^2+c^2*cos(theta)^2,
641 (a^2*cos(phi)^2+b^2*sin(phi)^2-c^2)*r*cos(theta)*sin(theta),
642 (b^2-a^2)*cos(phi)*sin(phi)*r*sin(theta)^2],
643 [(a^2*cos(phi)^2+b^2*sin(phi)^2-c^2)*r*cos(theta)*sin(theta),
644 r^2*((a^2*cos(phi)^2+b^2*sin(phi)^2)*cos(theta)^2+c^2*sin(theta)^2),
645 (b^2-a^2)*cos(phi)*sin(phi)*r^2*cos(theta)*sin(theta)],
646 [(b^2-a^2)*cos(phi)*sin(phi)*r*sin(theta)^2,
647 (b^2-a^2)*cos(phi)*sin(phi)*r^2*cos(theta)*sin(theta),
648 (a^2*sin(phi)^2+b^2*cos(phi)^2)*r^2*sin(theta)^2])
649 )
650 else if coordinate_system='cartesian4d then
651 (
652 dim:4,
653 ct_coords:[x,y,z,t],
654 if cframe_flag then fri:ident(4)
655 else lg:ident(4),
656 if cframe_flag then lfg:ident(4)
657 )
658 else if coordinate_system='spherical4d then
659 (
660 dim:4,
661 ct_coords:[r,theta,eta,phi],
662 /* declare(kurv,constant),*/
663 if cframe_flag then
664 (
665 fri:matrix([1 /* sqrt(1/(1-kurv*r^2)) */, 0, 0, 0],
666 [ 0, r, 0, 0],
667 [ 0, 0, r*sin(theta), 0],
668 [ 0, 0, 0, r*sin(eta)*sin(theta)]),
669 lfg:ident(4)
670 )
671 else lg:matrix([1 /* /(1-kurv*r^2) */, 0, 0, 0],
672 [0, r^2, 0, 0],
673 [0, 0, r^2*sin(theta)^2, 0],
674 [0, 0, 0, r^2*sin(eta)^2*sin(theta)^2])
675 )
676 else if coordinate_system='exteriorschwarzschild then
677 (
678 declare(m,constant),
679 assume(r>2*m),
680 dim:4,
681 ct_coords:[t,r,theta,phi],
682 if cframe_flag then
683 (
684 fri:matrix(
685 [sqrt((r-2*m)/r), 0,0, 0],
686 [ 0,sqrt(r/(r-2*m)),0, 0],
687 [ 0, 0,r, 0],
688 [ 0, 0,0,sin(theta)*r]),
689 lfg:ident(4),
690 lfg[1,1]:-1
691 )
692 else lg:matrix(
693 [(2*m-r)/r, 0, 0, 0],
694 [ 0,r/(r-2*m), 0, 0],
695 [ 0, 0,r^2, 0],
696 [ 0, 0, 0,sin(theta)^2*r^2])
697 )
698 else if coordinate_system='interiorschwarzschild then
699 (
700 dim:4,
701 ct_coords:[t,z,u,v],
702 declare([m],constant),
703 assume(t<2*m),
704 if cframe_flag then
705 (
706 fri:matrix(
707 [sqrt(t/(2*m-t)), 0, 0, 0],
708 [ 0, sqrt((2*m-t)/t), 0, 0],
709 [ 0, 0, t, 0],
710 [ 0, 0, 0, sin(u)*t]),
711 lfg:ident(4),
712 lfg[1,1]:-1
713 )
714 else lg:matrix([-t/(2*m-t), 0, 0, 0],
715 [ 0, (2*m-t)/t, 0, 0],
716 [ 0, 0, t^2, 0],
717 [ 0, 0, 0, sin(u)^2*t^2])
718 )
719 else if coordinate_system='kerr_newman then
720 (
721 dim:4,
722 declare([a,m,e],constant),
723 rho:sqrt(r^2+(a*cos(theta))^2),
724 delta:a^2-2*m*r+r^2+e^2,
725 ct_coords:[ct,r,theta,phi],
726 if cframe_flag then
727 (
728 fri:matrix(
729 [sqrt(delta)/rho, 0, 0, -a*sqrt(delta)*sin(theta)^2/rho],
730 [0, rho/sqrt(delta), 0, 0],
731 [0, 0, rho, 0],
732 [-a*sin(theta)/rho,0, 0, (r^2+a^2)*sin(theta)/rho]),
733 lfg:ident(4),
734 lfg[1,1]:-1
735 )
736 else lg:matrix(
737 [(a^2*sin(theta)^2-delta)/rho^2,0, 0,
738 a*sin(theta)^2*(delta-r^2-a^2)/rho^2],
739 [0, rho^2/delta, 0, 0],
740 [0, 0,rho^2, 0],
741 [a*sin(theta)^2*(delta-r^2-a^2)/rho^2,0,0,
742 ((r^2+a^2)^2-a^2*delta*sin(theta)^2)*sin(theta)^2/rho^2])
743 )
745 /* Now some undocumented metrics from metrics.mac and other sources */
747 else if coordinate_system='standard then
748 (
749 dim:4,
750 remove([a,d],constant),
751 depends([a,d],x),
752 ct_coords:[x,y,z,t],
753 if cframe_flag then error("No frames implementation")
754 else lg:matrix([a,0,0,0],[0,x^2,0,0],[0,0,x^2*sin(y)^2,0],[0,0,0,d])
755 )
756 else if coordinate_system='nondiagonal3d then
757 (
758 dim:3,
759 ct_coords:[x,y,z],
760 if cframe_flag then error("No frames implementation")
761 else lg:lg:matrix([x^2,y,0],[y,y^2,0],[0,0,z])
762 )
763 else if coordinate_system='alternateflat4d then
764 (
765 dim:4,
766 ct_coords:[x,y,z,t],
767 if cframe_flag then error("No frames implementation")
768 else lg:matrix([1-t^2,0,0,1-t*x],[0,1,0,0],[0,0,1,0],[1-t*x,0,0,1-x^2])
769 )
770 else if coordinate_system='eddington_finkelstein then
771 (
772 dim:4,
773 ct_coords:[r,theta,phi,v],
774 if cframe_flag then error("No frame fields yet")
775 else lg:matrix([0, 0, 0, 1],
776 [0, r^2, 0, 0],
777 [0, 0, sin(theta)^2*r^2, 0],
778 [1, 0, 0, (2*m-r)/r])
779 )
780 else if coordinate_system='isotropicschwarzschild then
781 (
782 dim:4,
783 ct_coords:[t,r,theta,phi],
784 if cframe_flag then
785 (
786 fri:matrix(
787 [(1-m/2/r)/(1+m/2/r), 0, 0, 0],
788 [0, (1+m/2/r)^2, 0, 0],
789 [0, 0, (1+m/2/r)^2*r, 0],
790 [0, 0, 0, (1+m/2/r)^2*r*sin(theta)]),
791 lfg:ident(4),
792 lfg[1,1]:-1
793 )
794 else lg:matrix([-((1-m/2/r)/(1+m/2/r))^2, 0, 0, 0],
795 [ 0, (1+m/2/r)^4, 0, 0],
796 [ 0, 0, (1+m/2/r)^4*r^2, 0],
797 [ 0, 0, 0, (1+m/2/r)^4*r^2*sin(theta)^2])
798 )
799 else if coordinate_system='vacuum then
800 (
801 dim:4,
802 ct_coords:[x,y,z,t],
803 if cframe_flag then error("No frames implementation")
804 else lg:matrix([0,1,a*(2*y-t),1],[1,0,a*(2*x-t),1],
805 [a*(2*y-t),a*(2*x-t),-6*a^2*(y+x)^2
806 +2*a^2*(2*x-t)*(2*y-t)-3/2,-a*(y+x+2*t)],
807 [1,1,-a*(y+x+2*t),1/2])
808 )
809 else if coordinate_system='bondi then
810 (
811 dim:4,
812 ct_coords:[x,y,z,t],
813 dependencies(a(x,y,t),b(x,y,t),c(x,y,t),d(x,y,t)),
814 if cframe_flag then error("No frames implementation")
815 else lg:matrix([0,0,0,%e^(2*b)],[0,-%e^(2*a)*x^2,0,%e^(2*a)*d*x^2],
816 [0,0,-%e^(-2*a)*x^2*sin(y)^2,0],
817 [%e^(2*b),%e^(2*a)*d*x^2,0,%e^(2*b)*c/x-%e^(2*a)*d^2*x^2])
818 )
819 else if coordinate_system='brans_dicke then
820 (
821 dim:4,
822 ct_coords:[x,y,z,t],
823 if cframe_flag then error("No frames implementation")
824 else lg:at(matrix([a,0,0,0],[0,a*x^2,0,0],
825 [0,0,a*x^2*sin(y)^2,0],[0,0,0,-d]),
826 [d=((1-b/x)/(1+b/x))^(2/l),
827 a=(1+b/x)^4*((1-b/x)/(1+b/x))^(2*(l-c-1)/l),
828 p=((1-b/x)/(1+b/x))^(c/l)])
829 /* ell:l=sqrt((c+1)^2-c*(1-w*c/2)) */
830 )
831 else if coordinate_system='lapides then
832 (
833 /* RICCI FLAT METRIC FROM QUANTUM GRAVITY-FOUND BY ALAN
834 LAPIDES- LOS ALAMOS APRIL 1982- PRESUME IT IS ONE OF
835 THE HARRISON METRICS */
836 dim:4,
837 ct_coords:[x,y,z,t],
838 if cframe_flag then error("No frames implementation")
839 else lg:at(matrix([1/c,0,0,0],[0,c*a,0,0],[0,0,c*b,0],[0,0,0,c*a]),
840 [c=%e^-(2*atan(sinh(t))),a=(cosh(t)^2-sin(y)^2)^2/cosh(t)^2,
841 b=cosh(t)^2*sin(y)^2])
842 )
843 else if coordinate_system='kantowski_sachs then
844 (
845 dim:4,
846 ct_coords:[r,theta,phi,t],
847 if cframe_flag then error("No frames implementation")
848 else
849 (
850 lg:matrix([X(t)^2,0,0,0],[0,Y(t)^2,0,0],
851 [0,0,Y(t)^2*sin(theta)^2,0],[0,0,0,-N^2]),
852 if member('misner,extra_args) then
853 lg:at(lg,[X(t)=exp(2*sqrt(3)*b(t)),Y(t)=exp(-2*sqrt(3)*o(t))]))
854 )
856 else if listp(coordinate_system) then
857 (
858 dim:length(coordinate_system)-1,
859 ct_coords:coordinate_system[dim+1],
860 if cframe_flag then
861 (
862 lfg:ident(dim),
863 fri:zeromatrix(dim,dim),
864 for i thru dim do for j thru dim do
865 fri[i,j]:diff(coordinate_system[i],ct_coords[j])
866 )
867 else
868 (
869 lg:zeromatrix(dim,dim),
870 for i thru dim do for j thru i do
871 (
872 lg[i,j]:sum(diff(coordinate_system[k%],ct_coords[i])*
873 diff(coordinate_system[k%],ct_coords[j]),k%,1,dim),
874 if i#j then lg[j,i]:lg[i,j]
875 )
876 )
877 )
878 else error("Unknown coordinate system type!"),
879 if member('cylindrical,extra_args) then
880 (
881 if member('z, ct_coords) then
882 error("Cannot add Z coordinate to this metric, as it already exists")
883 else
884 (
885 dim:dim+1,
886 ct_coords:append(ct_coords,['z]),
887 if cframe_flag then
888 (
889 fri:addrow(addcol(fri,zeromatrix(dim-1,1)),ident(dim)[dim]),
890 lfg:addrow(addcol(lfg,zeromatrix(dim-1,1)),ident(dim)[dim])
891 )
892 else lg:addrow(addcol(lg,zeromatrix(dim-1,1)),ident(dim)[dim])
893 )
894 )
895 else if member('minkowski,extra_args) then
896 (
897 if member('ct, ct_coords) then
898 error("Cannot add CT coordinate to this metric, as it already exists")
899 else
900 (
901 dim:dim+1,
902 /* ct_coords:append(ct_coords,['ct]),*/
903 ct_coords:cons('ct,ct_coords),
904 if cframe_flag then
905 (
906 /* fri:addrow(addcol(fri,zeromatrix(dim-1,1)),ident(dim)[dim]),
907 lfg:addrow(addcol(lfg,zeromatrix(dim-1,1)),-ident(dim)[dim])*/
908 fri:apply(matrix,cons(-ident(dim)[1],
909 makelist(transpose(apply(matrix,cons(zeromatrix(dim-1,dim-1)[1],
910 makelist(transpose(fri)[i],i,1,dim-1))))[j],j,1,dim-1))),
911 lfg:apply(matrix,cons(-ident(dim)[1],
912 makelist(transpose(apply(matrix,cons(zeromatrix(dim-1,dim-1)[1],
913 makelist(transpose(lfg)[i],i,1,dim-1))))[j],j,1,dim-1)))
914 )
915 /* else lg:addrow(addcol(lg,zeromatrix(dim-1,1)),-ident(dim)[dim])*/
916 else lg:apply(matrix,cons(-ident(dim)[1],
917 makelist(transpose(apply(matrix,cons(zeromatrix(dim-1,dim-1)[1],
918 makelist(transpose(lg)[i],i,1,dim-1))))[j],j,1,dim-1)))
919 )
920 ),
921 if member('all,extra_args) then
922 (
923 cmetric(),
924 christof(false)
925 ),
926 if verbose then
927 (
928 ldisplay(dim),
929 ldisplay(ct_coords),
930 if cframe_flag then
931 (
932 ldisplay(lfg),
933 ldisplay(fri)
934 )
935 else ldisplay(lg)
936 ),
937 done
938 )$
940 /***********************************
941 THE TENSORS OF GENERAL RELATIVITY
942 ***********************************/
944 /* Compute the Ricci rotation coefficients in a frame base */
945 trrc(disp):=block
946 (
947 [a,b,c,d],
948 array(_l,dim,dim,dim),
949 if not symmetricp(_lg(),dim) then
950 error("Frame base calculations require a symmetric metric."),
951 for a thru dim do for b thru dim do for c from b thru dim do
952 if b=c then _l[a,b,c]:0
953 else
954 (
955 /* Landau-Lifshitz II p98, pp379 (Hungarian edition) */
956 _l[a,b,c]:sum(sum((diff(fri[a][i%],ct_coords[k%])-
957 diff(fri[a][k%],ct_coords[i%]))*fr[b][i%]*fr[c][k%],
958 k%,1,dim),i%,1,dim),
959 if ctrgsimp then _l[a,b,c]:trigsimp(_l[a,b,c])
960 else if ratchristof then _l[a,b,c]:ratsimp(_l[a,b,c]),
961 _l[a,c,b]:-_l[a,b,c],
962 if disp#false and debugmode then ldisplay(_l[a,b,c])
963 ),
964 for c thru dim do for a thru dim do for b from a thru dim do
965 if a=b then lcs[c,a,b]:0
966 else
967 (
968 lcs[c,b,a]:(_l[a,b,c]+_l[b,c,a]-_l[c,a,b])/2,
969 if ctrgsimp then lcs[c,b,a]:trigsimp(lcs[c,b,a])
970 else if ratchristof then lcs[c,b,a]:ratsimp(lcs[c,b,a]),
971 if ratfac then lcs[c,b,a]:factor(lcs[c,b,a]),
972 lcs[c,a,b]:-lcs[c,b,a]
973 ),
974 remarray(_l),
975 if disp=all or disp=lcs then
976 for a thru dim do for b thru dim-1 do for c from b+1 thru dim do
977 if lcs[a,b,c]#0 then ldisplay(lcs[a,b,c])
978 );
980 /* Compute the Christoffel-symbols, or, if a frame base is used, the Ricci
981 rotation coefficients and place them in the arrays lcs[] and mcs[] */
982 christof(dis):=block
983 (
984 [symmmetric:symmetricp(_lg(),dim)],
985 if cframe_flag then
986 (
987 if not symmmetric then
988 error("Frame base calculations require a symmetric metric."),
989 trrc(dis)
990 ),
991 if (not symmmetric and dis = lcs) then
992 error("Christoffel symbols of the first kind require a symmetric metric"),
993 if (not symmmetric and ctorsion_flag) then
994 error("Torsion requires a symmetric metric"),
995 if (not symmmetric and cnonmet_flag) then
996 error("Nonmetricity requires a symmetric metric"),
997 if ctorsion_flag and not member('kt,arrays) then contortion(tr),
998 if cnonmet_flag and not member('nmc,arrays) then nonmetricity(nm),
999 if symmmetric then
1000 (
1001 for i thru dim do
1002 (
1003 for j from i thru dim do
1004 (
1005 if not cframe_flag then for k thru dim do
1006 (
1007 lcs[i,j,k]:(diff(lg[j,k],ct_coords[i])
1008 +diff(lg[k,i],ct_coords[j])-diff(lg[i,j],ct_coords[k]))/2,
1009 lcs[j,i,k]:lcs[i,j,k]
1010 ),
1011 for k thru dim do
1012 (
1013 mcs[i,j,k]:ctaylor
1014 (
1015 if diagmetric then ratsimp(_ug()[k,k]*lcs[i,j,k])
1016 else sum(_ug()[l%,k]*lcs[i,j,l%],l%,1,dim)
1017 ),
1018 if ratchristof then mcs[i,j,k]:ratsimp(mcs[i,j,k]),
1019 mcs[j,i,k]:mcs[i,j,k],
1020 if ctorsion_flag then
1021 (
1022 mcs[i,j,k]:mcs[i,j,k]-kt[i,j,k],
1023 if j#i then mcs[j,i,k]:mcs[j,i,k]-kt[j,i,k]
1024 ),
1025 if cnonmet_flag then
1026 (
1027 mcs[i,j,k]:mcs[i,j,k]-nmc[i,j,k],
1028 if j#i then mcs[j,i,k]:mcs[j,i,k]-nmc[j,i,k]
1029 )
1030 )
1031 )
1032 )
1033 )
1034 else block
1035 (
1036 /* This is a direct solution of the field equations,
1037 g_{ab,c}-g_{bd}Gamma^d_{ac}-g_{dc}Gamma^d_{ba}=0.
1038 See, e.g., Moffat & Boal, Phys Rev D V11 N6 p1375, Eq 2.5 */
1039 [_mcs:zeromatrix(dim^3,dim^3),_ucs,_abc,_gam,a,b,c,d],
1040 for a thru dim do
1041 for b thru dim do
1042 for c thru dim do
1043 for d thru dim do
1044 (
1045 _mcs[(a-1)*dim^2+(b-1)*dim+c,(a-1)*dim^2+(c-1)*dim+d]:_mcs[(a-1)*dim^2+(b-1)*dim+c,(a-1)*dim^2+(c-1)*dim+d]+lg[b,d],
1046 _mcs[(a-1)*dim^2+(b-1)*dim+c,(b-1)*dim^2+(a-1)*dim+d]:_mcs[(a-1)*dim^2+(b-1)*dim+c,(b-1)*dim^2+(a-1)*dim+d]+lg[d,c]
1047 ),
1048 _ucs:_mcs^^(-1),
1049 if ratchristof then _ucs:ratsimp(_ucs),
1050 _abc:zeromatrix(dim^3,1),
1051 for a thru dim do
1052 for b thru dim do
1053 for c thru dim do
1054 (
1055 _abc[(a-1)*dim^2+(b-1)*dim+c,1]:diff(lg[b,c],ct_coords[a]),
1056 if ratchristof then _abc[(a-1)*dim^2+(b-1)*dim+c,1]:ratsimp(_abc[(a-1)*dim^2+(b-1)*dim+c,1])
1057 ),
1058 _gam:_ucs._abc,
1059 if ratchristof then _gam:ratsimp(_gam),
1060 for a thru dim do
1061 for b thru dim do
1062 for c thru dim do
1063 mcs[a,b,c]:_gam[(a-1)*dim^2+(b-1)*dim+c][1]
1064 ),
1065 if dis = all or dis = lcs then for i thru dim do
1066 for j from i thru dim do
1067 for k thru dim do
1068 if lcs[i,j,k]#0 then ldisplay(lcs[i,j,k]),
1069 if dis = mcs or dis = all then for i thru dim do
1070 for j from (if symmmetric then i else 1) thru dim do
1071 for k thru dim do
1072 if mcs[i,j,k]#0 then ldisplay(mcs[i,j,k]),
1073 done
1074 )$
1076 /* In a moving frame, the Ricci-tensor is computed from the covariant
1077 curvature tensor */
1078 tricci(disp):=block
1079 (
1080 [df],
1081 if not symmetricp(_lg(),dim) then
1082 error("Frame base calculations require a symmetric metric."),
1083 if not member('lriem,arrays) then triemann(false),
1084 /* kill(ric), */
1085 df:true,
1086 for i thru dim do for j from i thru dim do
1087 (
1088 ric[i,j]:sum(sum(ufg[k%,l%]*lriem[l%,k%,i,j],k%,1,dim),l%,1,dim),
1089 if ctrgsimp then ric[i,j]:trigsimp(ric[i,j])
1090 else if ratfac then ric[i,j]:factor(ric[i,j])
1091 else ric[i,j]:ctaylor(ratsimp(ric[i,j])),
1092 if i#j then ric[j,i]:ric[i,j],
1093 if disp and ric[i,j]#0 then
1094 (
1095 ldisplay(ric[i,j]),
1096 df:false
1097 )
1098 ),
1099 tracer:sum(sum(ric[i%,j%]*ufg[i%,j%],i%,1,dim),j%,1,dim),
1100 if df and disp then print("THIS SPACETIME IS EMPTY AND/OR FLAT"),
1101 done
1102 );
1104 /* In case of a holonomic metric, the Christoffel-symbols of the
1105 second kind are used. */
1106 ricci(dis):=if cframe_flag then tricci(dis) else block
1107 (
1108 [suma,sumb,flat:true,symmmetric:symmetricp(_lg(),dim)],
1109 modedeclare(flat,boolean),
1110 if not member('mcs,arrays) then christof(false),
1111 if symmmetric then for i thru dim do
1112 (
1113 for j from i thru dim do
1114 (
1115 suma:0,sumb:0,
1116 for k thru dim do
1117 (
1118 if k#i then suma:suma+(diff(mcs[j,i,k],ct_coords[k])
1119 -diff(mcs[j,k,k],ct_coords[i])),
1120 sumb:sumb+sum(mcs[l%,k,k]*mcs[i,j,l%]-mcs[k,i,l%]*mcs[j,l%,k],l%,1,dim)
1121 ),
1122 ric[i,j]:ctaylor(suma+sumb),
1123 if ratfac then ric[i,j]:factor(ctaylor(ratsimp(ric[i,j]))),
1124 ric[j,i]:ric[i,j]
1125 )
1126 )
1127 else
1128 (
1129 for m thru dim do for n thru dim do
1130 ric[m,n]:ctaylor(sum(diff(mcs[m,n,a],ct_coords[a])-
1131 (diff(mcs[m,a,a]+mcs[a,m,a],ct_coords[n])+
1132 diff(mcs[n,a,a]+mcs[a,n,a],ct_coords[m]))/4+
1133 sum(-mcs[m,o,a]*mcs[a,n,o]+mcs[m,n,a]*mcs[a,o,o],o,1,dim)
1134 ,a,1,dim)),
1135 if ratfac then ric[i,j]:factor(ctaylor(ratsimp(ric[i,j])))
1136 ),
1137 if dis#false then
1138 (
1139 for i thru dim do for j from (if symmmetric then i else 1) thru dim do if ric[i,j]#0 then
1140 (
1141 flat:false,
1142 ldisplay(ric[i,j])
1143 ),
1144 if flat then print("THIS SPACETIME IS EMPTY AND/OR FLAT")
1145 ),
1146 tracer:if diagmetric
1147 then ctaylor(ratsimp(sum(ric[i%,i%]*ug[i%,i%],i%,1,dim)))
1148 else ctaylor(sum(sum(ric[i%,j%]*ug[i%,j%],i%,1,dim),j%,1,dim)),
1149 done
1150 )$
1152 /* Computing the mixed Ricci-tensor */
1153 uricci(dis):=block(
1154 [flat:true],
1155 modedeclare(flat,boolean),
1156 if not (member('ric,arrays) or matrixp(ric)) then ricci(false),
1157 for i thru dim do for j thru dim do uric[i,j]:0,
1158 for i thru dim do
1159 (
1160 for j thru dim do
1161 (
1162 uric[i,j]:if diagmetric then ratsimp(ctaylor(ric[i,j]*_ug()[j,j]))
1163 else ctaylor(sum(ric[i,k%]*_ug()[k%,j],k%,1,dim)),
1164 if ratfac then uric[i,j]:factor(ctaylor(ratsimp(uric[i,j])))
1165 )
1166 ),
1167 if dis#false then
1168 (
1169 for i thru dim do for j thru dim do if uric[i,j]#0 then
1170 (
1171 flat:false,
1172 ldisplay(uric[i,j])
1173 ),
1174 if flat then print("THIS SPACETIME IS EMPTY AND/OR FLAT")
1175 ),
1176 done
1177 )$
1179 /* The scalar curvature */
1181 scurvature():=
1182 (
1183 if not (member('uric,arrays) or matrixp(uric)) then uricci(false),
1184 if ratfac then factor(tracer) else tracer
1185 )$
1187 /* Computing the mixed Einstein-tensor */
1188 einstein(dis):=block
1189 (
1190 [flat:true],
1191 modedeclare(flat,boolean),
1192 if not (member('uric,arrays) or matrixp(uric)) then uricci(false),
1193 for i thru dim do for j thru dim do
1194 (
1195 if i=j then ein[i,j]:ctaylor(uric[i,j]-tracer/2)
1196 else ein[i,j]:ctaylor(uric[i,j]),
1197 if ratfac then ein[i,j]:factor(ratsimp(ein[i,j]))
1198 else if rateinstein then ein[i,j]:ratsimp(ein[i,j]),
1199 if dis#false and ein[i,j]#0 then
1200 (
1201 flat:false,
1202 ldisplay(ein[i,j])
1203 )
1204 ),
1205 if dis#false and flat then print("THIS SPACETIME IS EMPTY AND/OR FLAT"),
1206 done
1207 )$
1209 /* Computing the covariant Einstein-tensor */
1210 leinstein(dis):=block
1211 (
1212 if not (member('ein,arrays) or matrixp(ein)) then einstein(false),
1213 for i thru dim do for j thru i do
1214 (
1215 lein[i,j]:if diagmetric then ein[i,j]*_lg()[j,j]
1216 else sum(ein[i,k%]*_lg()[k%,j],k%,1,dim),
1217 if i#j then lein[j,i]:lein[i,j],
1218 if dis#false and lein[i,j]#0 then ldisplay(lein[i,j])
1219 )
1220 )$
1222 /* The Riemann-tensor is computed so that its indices are compatible with the
1223 index ordering in the itensor package:
1225 l _l _l _l _m _l _m
1226 R = | - | + | | - | |
1227 ijk ij,k ik,j mk ij mj ik
1229 This index ordering is somewhat unusual. In this form, the (3,1) Riemann
1230 tensor is antisymmetric in its second and third indices.
1232 The computation of the Weyl tensor follows the same index ordering.
1233 */
1235 /* This function computes the (1,3) Riemann tensor */
1236 riemann(dis):=if cframe_flag then block
1237 (
1238 [df],
1239 if not member('lriem,arrays) then triemann(false),
1240 df:true,
1241 for i thru dim do for j thru dim do for k thru dim do for l thru dim do
1242 (
1243 riem[i,j,k,l]:sum(lriem[i,j,k,m%]*ufg[l,m%],m%,1,dim),
1244 if dis and riem[i,j,k,l]#0 then
1245 (
1246 ldisplay(riem[i,j,k,l]),
1247 df:false
1248 )
1249 ),
1250 if df and dis then print("THIS SPACETIME IS FLAT"),
1251 done
1252 )
1253 else block
1254 (
1255 [flat : true],
1256 modedeclare(flat,boolean),
1257 if not member('mcs,arrays) then christof(false),
1258 for i thru dim do
1259 for j thru dim do
1260 for k thru dim do
1261 for l thru dim do riem[i,j,k,l]:0,
1262 for i thru dim do
1263 for j thru dim do
1264 for k thru j do
1265 for l thru dim do
1266 (
1267 riem[i,j,k,l]:ctaylor(diff(mcs[i,j,l], ct_coords[k]) -
1268 diff(mcs[i,k,l], ct_coords[j]) +
1269 sum(mcs[m%,k,l]*mcs[i,j,m%] -
1270 mcs[m%,j,l]*mcs[i,k,m%], m%, 1, dim)),
1271 if ratfac then riem[j,i,k,l] : factor(ctaylor(ratsimp(riem[j,i,k,l])))
1272 else if ratriemann then riem[j,i,k,l] : ctaylor(ratsimp(riem[j,i,k,l])),
1273 if k # j then riem[i,k,j,l] : -riem[i,j,k,l],
1274 if dis#false and riem[i,j,k,l] # 0 then
1275 (flat : false, ldisplay(riem[i,j,k,l]))
1276 ),
1277 if dis#false and flat then print("This spacetime is flat"),
1278 done
1279 )$
1281 /* Computing the covariant Riemann-tensor in a moving frame */
1282 triemann(disp):=block
1283 (
1284 [df],
1285 if not symmetricp(_lg(),dim) then
1286 error("Frame base calculations require a symmetric metric."),
1287 if not member('lcs,arrays) then trrc(false),
1288 /* kill(lriem), */
1289 df:true,
1290 for i thru dim do
1291 for j from i thru dim do
1292 for k thru dim do
1293 for l from k thru dim do
1294 (
1295 if i=j or k=l then lriem[l,j,i,k]:0
1296 /* Wald (1984), p51, eq. 3.4.21 */
1297 else lriem[l,j,i,k]:ctaylor(sum(fr[i][a%]*diff(lcs[j,k,l],ct_coords[a%])-
1298 fr[j][a%]*diff(lcs[i,k,l],ct_coords[a%]),a%,1,dim)-sum(sum(ufg[a%,b%]*
1299 (lcs[i,b%,k]*lcs[j,a%,l]-lcs[j,b%,k]*lcs[i,a%,l]+lcs[i,b%,j]*lcs[a%,k,l]-
1300 lcs[j,b%,i]*lcs[a%,k,l]),b%,1,dim),a%,1,dim)),
1301 if ctrgsimp then lriem[l,j,i,k]:trigsimp(ev(lriem[l,j,i,k]))
1302 else if ratriemann then lriem[l,j,i,k]:ratsimp(ev(lriem[l,j,i,k]))
1303 else lriem[l,j,i,k]:ev(lriem[l,j,i,k]),
1304 if ratfac then lriem[l,j,i,k]:factor(lriem[l,j,i,k]),
1305 lriem[l,i,j,k]:-lriem[l,j,i,k],
1306 lriem[k,j,i,l]:-lriem[l,j,i,k],
1307 lriem[k,i,j,l]:lriem[l,j,i,k],
1308 if disp and lriem[l,j,i,k]#0 then
1309 (
1310 ldisplay(lriem[l,j,i,k]),
1311 df:false
1312 )
1313 ),
1314 if df and disp then print("THIS SPACETIME IS FLAT"),
1315 done
1316 );
1317 /* This function computes the covariant Riemann tensor */
1318 lriemann(dis):=if cframe_flag then triemann(dis) else block
1319 (
1320 [symmmetric:symmetricp(_lg(),dim)],
1321 if not member('riem,arrays) then riemann(false),
1322 for i thru dim do
1323 for j thru dim do
1324 for k thru dim do
1325 for l thru dim do lriem[i,j,k,l]:0,
1326 for i thru dim do
1327 for j thru dim do
1328 for k thru (if symmmetric then j else dim) do
1329 for l thru (if symmmetric then i else dim) do
1330 (
1331 lriem[i,j,k,l] : if diagmetric then ratsimp(riem[i,j,k,l]*lg[l,l])
1332 else sum(riem[i,j,k,m%]*lg[m%,l], m%, 1, dim),
1333 if ratriemann then lriem[i,j,k,l] : ratsimp(lriem[i,j,k,l]),
1334 if symmmetric then
1335 (
1336 if j # k then lriem[i,k,j,l] : -lriem[i,j,k,l],
1337 if l # i then
1338 (
1339 lriem[l,j,k,i] : -lriem[i,j,k,l],
1340 if j # k then lriem[l,k,j,i] : lriem[i,j,k,l]
1341 )
1342 ),
1343 if dis#false and lriem[i,j,k,l] # 0 then ldisplay(lriem[i,j,k,l])
1344 ),
1345 done
1346 )$
1348 /* This function computes the contravariant Riemann tensor */
1349 uriemann(dis):=block
1350 (
1351 [symmmetric:symmetricp(_lg(),dim)],
1352 if not member('riem,arrays) then riemann(false),
1353 for i thru dim do
1354 for j thru dim do
1355 for k thru dim do
1356 for l thru dim do uriem[i,j,k,l]:0,
1357 for i thru dim do
1358 for j thru dim do
1359 for k thru (if symmmetric then j else dim) do
1360 for l thru (if symmmetric then i else dim) do
1361 (
1362 uriem[i,j,k,l] : if diagmetric then
1363 ratsimp(riem[i,j,k,l]*_ug()[i,i]*_ug()[j,j]*_ug()[k,k])
1364 else sum(sum(sum(
1365 riem[m%,n%,o%,l]*_ug()[i,m%]*_ug()[j,n%]*_ug()[k,o%],
1366 m%, 1, dim), n%, 1, dim), o%, 1, dim),
1367 if ratriemann then uriem[i,j,k,l] : ratsimp(uriem[i,j,k,l]),
1368 if symmmetric then
1369 (
1370 if j # k then uriem[i,k,j,l] : -uriem[i,j,k,l],
1371 if l # i then
1372 (
1373 uriem[l,j,k,i] : -uriem[i,j,k,l],
1374 if j # k then uriem[l,k,j,i] : uriem[i,j,k,l]
1375 )
1376 ),
1377 if dis#false and uriem[i,j,k,l] # 0 then ldisplay(uriem[i,j,k,l])
1378 ),
1379 done
1380 )$
1382 /* The Kretchmann invariant is R^2=lriem[i,j,k,l]*uriem[i,j,k,l] */
1383 rinvariant():=kinvariant:sum(sum(sum(sum(lriem[i%,j%,k%,l%]*uriem[i%,j%,k%,l%],
1384 i%,1,dim),j%,1,dim),k%,1,dim),l%,1,dim)$
1386 /* The covariant Weyl tensor. See the note on the Riemann tensor above
1387 for index ordering. */
1388 weyl(dis):=block
1389 (
1390 [flat:true,symmmetric:symmetricp(_lg(),dim)],
1391 modedeclare(flat,boolean),
1392 if dim=2 then
1393 (
1394 print("ALL 2 DIMENSIONAL SPACETIMES ARE CONFORMALLY FLAT"),
1395 return(done)
1396 ),
1397 if not ((member('ric,arrays) or matrixp(ric)) and member('lriem,arrays)) then
1398 (
1399 ricci(false),
1400 riemann(false),
1401 lriemann(false)
1402 ),
1403 for i thru dim do
1404 for j thru dim do
1405 for k thru dim do
1406 for l thru dim do weyl[i,j,k,l]:0,
1407 for i thru dim do
1408 for j from (if symmmetric then i+1 else 1) thru dim do
1409 for k from (if symmmetric then i else 1) thru dim do
1410 for l from (if symmmetric then k+1 else 1) thru (if (symmmetric and k=i) then j else dim) do
1411 (
1412 weyl[l,j,i,k]:ctaylor(lriem[l,j,i,k]+
1413 1/(dim-2)*(_lg()[i,l]*ric[j,k]-_lg()[i,k]*ric[l,j]
1414 +_lg()[j,k]*ric[l,i]-_lg()[j,l]*ric[k,i])-
1415 tracer/((dim-1)*(dim-2))*
1416 (_lg()[i,l]*_lg()[k,j]-_lg()[i,k]*_lg()[l,j])
1417 ),
1418 if ratfac then weyl[l,j,i,k]:factor(ratsimp(weyl[l,j,i,k]))
1419 else if ratweyl then weyl[l,j,i,k]:ratsimp(weyl[l,j,i,k]),
1420 if symmmetric then
1421 (
1422 weyl[k,i,j,l]:weyl[l,j,i,k],
1423 weyl[k,j,i,l]:weyl[l,i,j,k]:-weyl[l,j,i,k],
1424 if i#k or j>l then
1425 (
1426 weyl[j,l,k,i]:weyl[i,k,l,j]:weyl[l,j,i,k],
1427 weyl[j,k,l,i]:weyl[i,l,k,j]:-weyl[l,j,i,k]
1428 )
1429 ),
1430 if dis#false and weyl[l,j,i,k]#0 then
1431 (
1432 flat:false,
1433 ldisplay(weyl[l,j,i,k])
1434 )
1435 ),
1436 if dis#false and flat then print("THIS SPACETIME IS CONFORMALLY FLAT"),
1437 done
1438 )$
1440 /***********************************
1441 MISCELLANEOUS CALCULATIONS
1442 ***********************************/
1444 /* The geodesic equation of motion */
1445 cgeodesic(dis):=block
1446 (
1447 [s,lgeod],
1448 map(lambda([u],apply('remove,[u,constant])),ct_coords),
1449 depends(ct_coords,s),
1450 for i thru dim do
1451 lgeod[i]:if diagmetric then
1452 ratsimp(lg[i,i]*diff(ct_coords[i],s,2)+
1453 sum(diff(lg[i,i],ct_coords[a%])*diff(ct_coords[i],s)*
1454 diff(ct_coords[a%],s),a%,1,dim)-
1455 1/2*sum(diff(lg[a%,a%],ct_coords[i])*
1456 diff(ct_coords[a%],s)^2,a%,1,dim))
1457 else ratsimp(sum(lg[i,a%]*diff(ct_coords[a%],s,2),a%,1,dim)+
1458 sum(sum(diff(lg[i,b%],ct_coords[a%])*diff(ct_coords[b%],s)*
1459 diff(ct_coords[a%],s)-
1460 1/2*diff(lg[a%,b%],ct_coords[i])*diff(ct_coords[a%],s)*
1461 diff(ct_coords[b%],s),a%,1,dim),b%,1,dim)),
1462 for i thru dim do
1463 (
1464 geod[i]:factor(ratexpand(sum(ratsimp(_ug()[i,a%]*lgeod[a%]),a%,1,dim))),
1465 if dis#false then ldisplay(geod[i])
1466 ),
1467 done
1468 )$
1470 /* The 4-dimensional Brans-Dicke vacuum field equations are represented by the
1471 array bd and the scalar equation is generated by setting the d'Alembertian
1472 of the scalar field to zero. That is, one calls the function dscalar on the
1473 scalar field. */
1474 bdvac(zz):=block
1475 (
1476 [addd,boxq:0],
1477 if dim#4 then error("This program is restricted to 4 dimensions"),
1478 for i:1 thru 4 do for j:1 thru 4 do
1479 addd[i,j]:w/zz^2*(diff(zz,ct_coords[i])*diff(zz,ct_coords[j])-lg[i,j]*
1480 sum(diff(zz,ct_coords[k%])*diff(zz,ct_coords[k%])*_ug()[k%,k%],k%,1,4)/2)+
1481 (diff(diff(zz,ct_coords[i]),ct_coords[j])-sum(mcs[i,j,k%]*
1482 diff(zz,ct_coords[k%]),k%,1,4)-lg[i,j]*boxq)/zz,
1483 for i:1 thru 4 do for j:1 thru 4 do
1484 bd[i,j]:ratsimp(ric[i,j]-r*lg[i,j]/2-0*t[i,j]-addd[i,j]),
1485 remarray(addd),
1486 done
1487 )$
1489 /* Compute the Euler-Lagrange equations for the density of the invariant R^2.
1490 The form is (D is the Kronecker-delta):
1492 b 2 b b bc
1493 (1/2)*D R - 2 R R + 2 D []R -2 g R
1494 a a a ;ac
1496 The equations are sometimes less complex if tracer is given a parametric
1497 name with dependencies corresponding to those of the scalar curvature. */
1498 invariant1():=block
1499 (
1500 for aa thru dim do for bb thru dim do
1501 (
1502 inv1[aa,bb]:0
1503 ),
1504 if diagmetric then for aa thru dim do for bb thru dim do
1505 (
1506 inv1[aa,bb]:ratsimp(kdelt(aa,bb)*tracer^2/2-2*tracer*uric[aa,bb]-
1507 2*kdelt(aa,bb)*dscalar(tracer)+2*_ug()[aa,aa]*
1508 (diff(diff(tracer,ct_coords[bb]),ct_coords[aa])-
1509 sum(mcs[bb,aa,k%]*diff(tracer,ct_coords[k%]),k%,1,dim)))
1510 )
1511 else for aa thru dim do for bb thru dim do
1512 (
1513 inv1[aa,bb]:ratsimp(kdelt(aa,bb)*tracer^2/2-2*tracer*uric[aa,bb]-
1514 2*kdelt(aa,bb)*dscalar(tracer)+2*sum(_ug()[bb,c%]*
1515 (diff(diff(tracer,ct_coords[aa]),ct_coords[c%])-sum(
1516 mcs[aa,c%,k%]*diff(tracer,ct_coords[k%]),k%,1,dim)),c%,1,dim))
1517 )
1518 )$
1520 invariant2():="NOT YET IMPLEMENTED";
1521 bimetric():="NOT YET IMPLEMENTED";
1523 /* Compute a Newman-Penrose null tetrad */
1524 nptetrad(disp):=block(
1525 [nu],
1526 if cframe_flag#true then error("Frame base required"),
1527 if dim#4 then error("This calculation requires dim=4"),
1528 nu:ident(4),
1529 nu[1,1]:-1,
1530 if lfg#nu then error("Frame metric signature of (-,+,+,+) required"),
1531 np:ident(4),
1532 np[1]:(fri[1]+fri[2])/sqrt(2),
1533 np[2]:(fri[1]-fri[2])/sqrt(2),
1534 np[3]:(fri[3]+%i*fri[4])/sqrt(2),
1535 np[4]:(fri[3]-%i*fri[4])/sqrt(2),
1536 if ctrgsimp then np:trigsimp(np)
1537 else np:ratsimp(np),
1538 if ratfac then np:factor(np),
1539 if disp then ldisplay(np),
1540 npi:np.ug,
1541 if ctrgsimp then npi:trigsimp(npi)
1542 else npi:ratsimp(npi),
1543 if ratfac then npi:factor(npi),
1544 if disp then ldisplay(npi),
1545 done
1546 );
1548 contract4(w,l1,l2,l3,l4)::=
1549 (
1550 if dim#4 then error("This calculation requires dim=4"),
1551 modedeclare([l1,l2,l3,l4],fixnum),
1552 block
1553 (
1554 [ans:0,i,j,k,l],
1555 for i thru 4 do for j thru 4 do for k thru 4 do for l thru 4 do
1556 ans:ans+w[l,j,i,k]*npi[l1,i]*npi[l2,j]*npi[l3,k]*npi[l4,l],
1557 if ctrgsimp then ans:trigsimp(ans)
1558 else ratsimp(ans)
1559 )
1560 );
1562 /* Compute the Newman-Penrose coefficients */
1563 psi(disp):=block
1564 (
1565 /* To compute these values, we must transform the Weyl-tensor to
1566 coordinate representation if we've been working in a frame field. */
1567 [w],
1568 if cframe_flag then
1569 (
1570 for i thru dim do for j thru dim do for k thru dim do for l thru dim do
1571 w[i,j,k,l]:0,
1572 for i thru dim do
1573 for j from i+1 thru dim do
1574 for k from i thru dim do
1575 for l from k+1 thru (if k=i then j else dim) do
1576 (
1577 w[l,j,i,k]:sum(sum(sum(sum(weyl[l%,j%,i%,k%]*fri[i%,i]*fri[j%,j]*
1578 fri[k%,k]*fri[l%,l],i%,1,dim),j%,1,dim),k%,1,dim),l%,1,dim),
1579 if ratfac then w[l,j,i,k]:factor(ratsimp(w[l,j,i,k]))
1580 else if ratweyl then w[l,j,i,k]:ratsimp(w[l,j,i,k]),
1581 w[k,i,j,l]:w[l,j,i,k],
1582 w[k,j,i,l]:w[l,i,j,k]:-w[l,j,i,k],
1583 if i#k or j>l then
1584 (
1585 w[j,l,k,i]:w[i,k,l,j]:w[l,j,i,k],
1586 w[j,k,l,i]:w[i,l,k,j]:-w[l,j,i,k]
1587 )
1588 )
1589 )
1590 else w:'weyl,
1591 psi[0]:ratsimp(-contract4(w,1,3,1,3)),
1592 if disp then ldisplay(psi[0]),
1593 psi[1]:ratsimp(-contract4(w,1,2,1,3)),
1594 if disp then ldisplay(psi[1]),
1595 psi[2]:ratsimp(gfactor(-1/2*(contract4(w,1,2,1,2)+contract4(w,1,2,3,4)))),
1596 if disp then ldisplay(psi[2]),
1597 psi[3]:ratsimp(contract4(w,1,2,2,4)),
1598 if disp then ldisplay(psi[3]),
1599 psi[4]:ratsimp(-contract4(w,2,4,2,4)),
1600 if disp then ldisplay(psi[4]),
1601 done
1602 );
1604 /* Compute the Petrov-classification of the metric; algorithm based on
1605 Classifying geometries in general relativity: III Classification in practice
1606 by Pollney, Skea, d'Inverno, Class. Quant. Grav. 17 2885-2992 (2000)
1608 The algorithm has been optimized to compute only quantities that are
1609 required for a specific (sub)case. As of 12/19/2004, it is essentially
1610 untested, may even contain typos in some expressions.
1611 */
1612 petrov():=block
1613 (
1614 [
1615 T:[ 0,'N,'II,'III, 'D,'II,'II, 7,
1616 'II,'I, 'I, 11,'II, 13, 14,15,
1617 'N,'I, 'I, 19,'II, 21, 13,23,
1618 'III,19, 11, 27, 7, 23, 15,31],
1619 P:1
1620 ],
1622 if psi[4]#0 then P:P+1,
1623 if psi[3]#0 then P:P+2,
1624 if psi[2]#0 then P:P+4,
1625 if psi[1]#0 then P:P+8,
1626 if psi[0]#0 then P:P+16,
1628 if numberp(T[P]) and T[P]>0 then
1629 (
1630 if T[P]=7 then
1631 if ratsimp(psi[3]^2-3*psi[2]*psi[4])=0 then 'D
1632 else 'II
1633 else if T[P]=11 then
1634 if ratsimp(27*psi[4]^2*psi[1]+64*psi[3]^3)=0 then 'II
1635 else 'I
1636 else if T[P]=13 then
1637 if ratsimp(psi[1]^2*psi[4]+2*psi[2]^3)=0 then 'II
1638 else 'I
1639 else if T[P]=14 then
1640 if ratsimp(9*psi[2]^2-16*psi[1]*psi[3])=0 then 'II
1641 else 'I
1642 else if T[P]=15 then
1643 if ratsimp(3*psi[2]^2-4*psi[1]*psi[3])=0 and
1644 ratsimp(psi[2]*psi[3]-3*psi[1]*psi[4])=0 then 'II
1645 else 'I
1646 else if T[P]=19 then
1647 if ratsimp(psi[0]*psi[4]^3-27*psi[3]^4)=0 then 'II
1648 else 'I
1649 else if T[P]=21 then
1650 if ratsimp(9*psi[2]^2-psi[4]^2)=0 then 'D
1651 else 'I
1652 else if T[P]=23 then block
1653 (
1654 [I:ratsimp(psi[0]*psi[4]+3*psi[2]^2)],
1655 if I=0 and ratsimp(4*psi[2]*psi[4]-3*psi[3]^2)=0 then 'III
1656 else block
1657 (
1658 [J:ratsimp(4*psi[2]*psi[4]-3*psi[3]^2)],
1659 if ratsimp(psi[4]*I^2-3*J*(psi[0]*J-2*psi[2]*I))=0 then 'II
1660 else 'I
1661 )
1662 )
1663 else if T[P]=27 then
1664 if ratsimp(psi[0]*psi[3]^2-psi[1]^2*psi[4])=0 then
1665 if ratsimp(psi[0]*psi[4]+2*psi[1]*psi[3])=0 then 'D
1666 else if ratsimp(psi[0]*psi[4]-16*psi[1]*psi[3])=0 then 'II
1667 else 'I
1668 else block
1669 (
1670 [I:ratsimp(psi[0]*psi[4]+2*psi[1]*psi[3])],
1671 if I=0 then block
1672 (
1673 [J:ratsimp(-psi[0]*psi[3]^2-psi[1]^2*psi[4])],
1674 if J=0 then 'III
1675 else if ratsimp(I^3-27*J^2)=0 then 'II
1676 else 'I
1677 )
1678 else 'I
1679 )
1680 else block
1681 (
1682 [H:ratsimp(psi[0]*psi[2]-psi[1]^2)],
1683 if H=0 then
1684 if ratsimp(psi[0]*psi[3]-psi[1]*psi[2])=0 then
1685 if ratsimp(psi[0]*psi[4]-psi[2]^2)=0 then 'N
1686 else 'I
1687 else block
1688 (
1689 [E:ratsimp(psi[0]*psi[4]-psi[2]^2)],
1690 if E=0 then
1691 if ratsimp(37*psi[2]^2+27*psi[1]*psi[3])=0 then 'II
1692 else 'I
1693 else block
1694 (
1695 [
1696 A:ratsimp(psi[1]*psi[3]+psi[2]^2),
1697 I:ratsimp(E-4*A)
1698 ],
1699 if I#0 and ratsimp(I^3-27*(psi[4]*H-psi[3]^2*psi[0]+
1700 psi[1]*psi[2]*psi[3]+psi[2]*A)^2)=0 then 'II
1701 else 'I
1702 )
1703 )
1704 else block
1705 (
1706 [I:ratsimp(psi[0]*psi[4]-psi[2]^2-4*(psi[1]*psi[3]+psi[2]^2))],
1707 if I=0 then
1708 if ratsimp(psi[4]*H-psi[3]^2*psi[0]+psi[1]*psi[2]*psi[3]+
1709 psi[2]*(psi[1]*psi[3]+psi[2]^2))=0 then 'III
1710 else 'I
1711 else
1712 if ratsimp(psi[0]^2*psi[3]-psi[0]*psi[1]*psi[2]-2*psi[1]*H)=0 then
1713 if ratsimp(psi[0]^2*I-12*H^2)=0 then 'D
1714 else if ratsimp(psi[0]^2*I-3*H^2)=0 then 'II
1715 else 'I
1716 else block
1717 (
1718 [J:ratsimp(psi[4]*H-psi[3]^2*psi[0]+psi[1]*psi[2]*psi[3]+
1719 psi[2]*(psi[1]*psi[3]+psi[2]^2))],
1720 if J#0 and ratsimp(I^3-27*J^2)=0 then 'II
1721 else 'I
1722 )
1723 )
1724 )
1725 )
1726 else T[P]
1727 );
1729 /* Old algorithm by anonymous, intended for diagonal metrics
1730 petrov():=block
1731 (
1732 [ii,jj,gg,hh],
1733 ii:ratsimp(psi[0]*psi[4]-4*psi[1]*psi[3]+3*psi[2]^2),
1734 jj:ratsimp(determinant(matrix([psi[0],psi[1],psi[2]],
1735 [psi[1],psi[2],psi[3]],
1736 [psi[2],psi[3],psi[4]]))),
1737 if ratsimp(ii^3-27*jj^2)#0 then "Type is I"
1738 else
1739 (
1740 gg:ratsimp(psi[0]^2*psi[3]-3*psi[0]*psi[1]*psi[2]+2*psi[1]^3),
1741 if ii=0 and jj=0 then
1742 (
1743 if gg=0 then
1744 (
1745 hh:ratsimp(psi[0]*psi[2]-psi[1]^2),
1746 if hh=0 then
1747 (
1748 if psi[0] = 0 and psi[1] = 0 and psi[2] = 0 and psi[3] = 0 and
1749 psi[4] = 0 then "Type is FLAT(0)"
1750 else "Type is N"
1751 )
1752 else "Type is III"
1753 )
1754 else "Type is III"
1755 )
1756 else
1757 (
1758 hh:ratsimp(psi[0]*psi[2]-psi[1]^2),
1759 if ratsimp(gg-(psi[0]^2*ii-12*hh^2))#0 then "Type is II"
1760 else "Type is D"
1761 )
1762 )
1763 )$*/
1765 frame_bracket(fr,fri,diagframe):=block
1766 (
1767 for a thru dim do for b thru dim do for c thru dim do
1768 fb[a,b,c]:sum(sum(fr[a,e%]*(diff(fri[e%,c],ct_coords[f%])-
1769 diff(fri[f%,c],ct_coords[e%]))*fr[b,f%],e%,1,dim),f%,1,dim)
1770 );
1772 /* The findde function finds a list of unique differential equations in the
1773 list or 2-dimensional array a. The variable deindex is a global list
1774 of the indices corresponding to these unique differential equations. */
1775 findde(a,n):=
1776 (
1777 modedeclare(n,fixnum),
1778 if n=1 then findde1(a)
1779 else if n=2 then findde2(a)
1780 else if n=3 then findde3(a)
1781 else error("Invalid dimension:",n)
1782 )$
1784 findde1(list):=block
1785 (
1786 [inflag:true,deriv:nounify('derivative),l:[],l1,q],
1787 deindex:[],
1788 for i while list#[] do
1789 (
1790 l1:factor(num(first(list))),
1791 list:rest(list),
1792 if not freeof(deriv,l1) then
1793 (
1794 deindex:cons(i,deindex),
1795 if inpart(l1,0)#"*" then l:cons(l1,l)
1796 else
1797 (
1798 q:1,
1799 for j thru length(l1) do
1800 if not freeof(deriv,inpart(l1,j)) then q:q*inpart(l1,j),
1801 l:cons(q,l)
1802 )
1803 )
1804 ),
1805 cleanup(l)
1806 )$
1808 findde2(a):=block
1809 (
1810 [inflag:true,deriv:nounify('derivative),l:[],t,q],
1811 deindex:[],
1812 for i thru dim do for j thru dim do
1813 (
1814 t:factor(num(a[i,j])),
1815 if not freeof(deriv,t) then
1816 (
1817 deindex:cons([i,j],deindex),
1818 if inpart(t,0)#"*" then l:cons(t,l)
1819 else
1820 (
1821 q:1,
1822 for n thru length(t) do
1823 if not freeof(deriv,inpart(t,n)) then q:q*inpart(t,n),
1824 l:cons(q,l)
1825 )
1826 )
1827 ),
1828 cleanup(l)
1829 )$
1831 findde3(a):=block
1832 (
1833 [inflag:true,deriv:nounify('derivative),l:[],t,q],
1834 deindex:[],
1835 for i thru dim do for j thru dim do for k thru dim do
1836 (
1837 t:factor(num(a[i,j,k])),
1838 if not freeof(deriv,t) then
1839 (
1840 deindex:cons([i,j,k],deindex),
1841 if inpart(t,0)#"*" then l:cons(t,l)
1842 else
1843 (
1844 q:1,
1845 for n thru length(t) do
1846 if not freeof(deriv,inpart(t,n)) then q:q*inpart(t,n),
1847 l:cons(q,l)
1848 )
1849 )
1850 ),
1851 cleanup(l)
1852 )$
1854 cleanup(ll):=block
1855 (
1856 [a,l:[],index:[]],
1857 while ll#[] do
1858 (
1859 a:first(ll),
1860 ll:rest(ll),
1861 if not member(a,ll) then
1862 (
1863 l:cons(a,l),
1864 index:cons(first(deindex),index)
1865 ),
1866 deindex:rest(deindex)
1867 ),
1868 deindex:index,
1869 l
1870 )$
1872 /***********************************
1873 AUXILIARY FUNCTIONS (not used in ctensor)
1874 ***********************************/
1876 /* Covariant and contravariant gradients */
1877 cograd(f,xyz):=block
1878 (
1879 [],
1880 for mm thru dim do
1881 (
1882 arraysetapply(xyz,[mm],ratsimp(diff(f,ct_coords[mm])))
1883 )
1884 )$
1886 contragrad(f,xyz):=block
1887 (
1888 if diagmetric then for mm thru dim do
1889 arraysetapply(xyz,[mm],ratsimp(ratsimp(_ug()[mm,mm]*diff(f,ct_coords[mm]))))
1890 else for mm thru dim do
1891 arraysetapply(xyz,[mm],sum(_ug()[mm,n%]*diff(f,ct_coords[n%]),n%,1,dim))
1892 )$
1894 /* The d'Alembertian of a scalar */
1895 dscalar(phi):=block
1896 (
1897 if diagmetric then ratsimp(1/sqrt(-gdet)*sum(diff(_ug()[i%,i%]*
1898 sqrt(-gdet)*diff(phi,ct_coords[i%]),ct_coords[i%]),i%,1,dim))
1899 else ratsimp(1/sqrt(-gdet)*sum(sum(diff(_ug()[i%,j%]*sqrt(-gdet)*
1900 diff(phi,ct_coords[j%]),ct_coords[i%]),i%,1,dim),j%,1,dim))
1901 )$
1903 /* Covariant divergence. The object gxyz must be a mixed 2nd rank tensor
1904 whose first index is covariant. */
1905 checkdiv(gxyz):=block
1906 (
1907 modedeclare([i,j%],fixnum),
1908 if diagmatrixp(gxyz,dim) then for i thru dim do
1909 print(div[i]:ratsimp(ctaylor(1/sqrt(-gdet)*diff(sqrt(-gdet)*gxyz[i,i],
1910 ct_coords[i])-sum(mcs[i,j%,j%]*gxyz[j%,j%],j%,1,dim))))
1911 else for i thru dim do
1912 print(div[i]:ratsimp(ctaylor(1/sqrt(-gdet)*sum(diff(sqrt(-gdet)*gxyz[i,j%],
1913 ct_coords[j%]),j%,1,dim)-sum(sum(mcs[i,j%,a%]*gxyz[a%,j%],a%,1,dim),j%,1,dim))))
1914 )$
1916 /* Contortion. Mostly generated with ic_convert(). */
1917 contortion(tr):=block
1918 (
1919 for i thru dim do
1920 for j thru dim do
1921 for k thru dim do kt[i,j,k]:sum(_ug()[k,l%]*
1922 (-sum(tr[l%,i,m%]*lg[j,m%],m%,1,dim)-sum(lg[l%,m%]*tr[i,j,m%],m%,1,dim)-
1923 sum(tr[l%,j,m%]*lg[i,m%],m%,1,dim)), l%,1,dim)/2
1924 )$
1926 /* Nonmetricity. Also generated mainly with ic_convert(). */
1927 nonmetricity(nm):=block
1928 (
1929 for i thru dim do
1930 for j thru dim do
1931 for k thru dim do nmc[i,j,k]:-nm[i]*kdelt(j,k)-kdelt(i,k)*nm[j]+
1932 sum(_ug()[k,l%]*nm[l%],l%,1,dim)*lg[i,j]
1933 )$
1936 /********************************************************************
1937 * Untested code from miscellaneous files...
1939 /* The second covariant derivative of the covariant Ricci Tensor */
1941 riccicov(ii,jj,kk,ll):=block(
1942 modedeclare([ii,jj,kk,ll,a%,b%],fixnum),
1943 if diagmetric then h[ii,jj,kk,ll]:h[jj,ii,kk,ll]:
1944 ratsimp(sum(
1945 mcs[a%,ii,jj]*ric[jj,jj]*mcs[kk,ll,a%]+
1946 mcs[a%,jj,ii]*ric[ii,ii]*mcs[kk,ll,a%]+
1947 ric[a%,a%]*mcs[ii,kk,a%]*mcs[jj,ll,a%]+
1948 mcs[a%,kk,ii]*ric[ii,ii]*mcs[jj,ll,a%]+
1949 ric[a%,a%]*mcs[ii,ll,a%]*mcs[jj,kk,a%]+
1950 mcs[a%,kk,jj]*ric[jj,jj]*mcs[ii,ll,a%]-
1951 diff(ric[ii,jj],ct_coords[a%],1)*mcs[kk,ll,a%],a%,1,dim)
1953 -(diff(ric[ii,ii],ct_coords[kk],1)*mcs[jj,ll,ii]+
1954 ric[ii,ii]*diff(mcs[jj,kk,ii],ct_coords[ll],1)+
1955 diff(ric[ii,ii],ct_coords[ll],1)*mcs[jj,kk,ii]+
1956 diff(ric[jj,jj],ct_coords[kk],1)*mcs[ii,ll,jj]+
1957 ric[jj,jj]*diff(mcs[ii,kk,jj],ct_coords[ll],1)+
1958 diff(ric[jj,jj],ct_coords[ll],1)*mcs[ii,kk,jj])
1960 +diff(diff(ric[ii,jj],ct_coords[kk],1),ct_coords[ll],1))
1962 else h[ii,jj,kk,ll]:h[jj,ii,kk,ll]:
1963 ratsimp(sum(sum(
1964 mcs[a%,ii,b%]*ric[b%,jj]*mcs[kk,ll,a%]+
1965 mcs[a%,jj,b%]*ric[b%,ii]*mcs[kk,ll,a%]+
1966 ric[a%,b%]*mcs[ii,kk,a%]*mcs[jj,ll,b%]+
1967 mcs[a%,kk,b%]*ric[b%,ii]*mcs[jj,ll,a%]+
1968 ric[a%,b%]*mcs[ii,ll,a%]*mcs[jj,kk,b%]+
1969 mcs[a%,kk,b%]*ric[b%,jj]*mcs[ii,ll,a%],a%,1,dim),b%,1,dim)
1971 - sum(
1972 diff(ric[ii,jj],ct_coords[a%],1)*mcs[kk,ll,a%]+
1973 diff(ric[a%,ii],ct_coords[kk],1)*mcs[jj,ll,a%]+
1974 ric[a%,ii]*diff(mcs[jj,kk,a%],ct_coords[ll],1)+
1975 diff(ric[a%,ii],ct_coords[ll],1)*mcs[jj,kk,a%]+
1976 diff(ric[a%,jj],ct_coords[kk],1)*mcs[ii,ll,a%]+
1977 ric[a%,jj]*diff(mcs[ii,kk,a%],ct_coords[ll],1)+
1978 diff(ric[a%,jj],ct_coords[ll],1)*mcs[ii,kk,a%],a%,1,dim)
1980 +diff(diff(ric[ii,jj],ct_coords[kk],1),ct_coords[ll],1)))$
1982 *******************************************************************/
1986 put('ctensor,'v20211011,'version)$