1;;; -*- Mode:LISP; Package: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: Code to generate ctensor programs from itensor expressions 14;; 15 16; ** (c) Copyright 1979 Massachusetts Institute of Technology ** 17 18(in-package :maxima) 19 20 21(declare-top (special $imetric $metricconvert indlist empty)) 22 23;$METRICCONVERT if non-NIL will allow $IC_CONVERT to rename the metric tensor 24; ($IMETRIC must be bound) with 2 covariant indices to LG and with 2 25; contravariant indices to UG. 26 27(defun $ic_convert (ee) 28 (prog (e free lhs rhs) 29 (setq e ($expand ee)) 30 (cond ((or (atom e) (not (eq (caar e) 'mequal))) 31 (merror "IC_CONVERT requires an equation as an argument")) 32 ((equal (setq free ($indices e)) empty) 33 (return (cons '(msetq) (cdr e)))) 34 ((or (symbolp (cadr e)) ;If a symbol or 35 (and (rpobj (cadr e)) ;an indexed object with no dummy 36 (null (cdaddr ($indices2 (cadr e)))))) ;indices 37 (setq lhs (cadr e) rhs (caddr e))) 38 ((or (symbolp (caddr e)) 39 (and (rpobj (caddr e)) 40 (null (cdaddr ($indices2 (caddr e)))))) 41 (setq lhs (caddr e) rhs (cadr e))) 42 (t (merror "At least one side of the equation must be a~ 43 ~%symbol or a single indexed object"))) 44 (cond ((and (not (symbolp lhs)) 45 (not (null (cdddr lhs)))) 46 (merror "Cannot assign to indexed objects with derivative ~ 47 indices:~%~M" 48 (ishow lhs)))) 49 (setq free (nreverse (itensor-sort (cdadr free))) ;Set FREE to just the 50 indlist nil) ;free indices 51 (and $metricconvert (boundp '$imetric) 52 (setq lhs (changename $imetric t 0 2 '$ug 53 (changename $imetric t 2 0 '$lg lhs)) 54 rhs (changename $imetric t 0 2 '$ug 55 (changename $imetric t 2 0 '$lg rhs)))) 56 (tabulate rhs) 57 (setq indlist (unique indlist)) 58 (do ((q (mapcar 'car indlist) (cdr q))) 59 ((null q)) 60 (cond ((member (car q) (cdr q) :test #'eq) 61 (merror "~ 62IC_CONVERT cannot currently handle indexed objects of the same name~ 63~%with different numbers of covariant and//or contravariant indices:~%~M" 64 (car q))))) 65 (cond ((not (symbolp lhs)) 66 (do ((test) (flag) (name)) 67 (flag) 68 (setq test (list (caar lhs) (length (cdadr lhs)) 69 (length (cdaddr lhs)))) 70 (cond ((or (member test indlist :test #'equal) 71 (not (member (car test) 72 (mapcar 'car indlist) :test #'eq))) 73 (setq flag t)) 74 (t 75 (mtell "Assignment is to be made to ~M~ 76~%This name with a different number of covariant and//or contravariant~ 77~%indices appears on the other side of the equation. To avoid array name~ 78~%conflicts, choose a new name for this object:~%" 79 (ishow lhs)) 80 (cond ((not (symbolp (setq name (retrieve nil nil)))) 81 (merror "Name not an atom"))) 82 (setq lhs (cons (ncons name) (cdr lhs)))))))) 83 (return (do ((free free (cdr free)) 84 (equ (cons '(msetq) (list (changeform lhs) 85 (t-convert 86 (summer1 rhs)))))) 87 ((null free) equ) 88 (setq equ (append '((mdo)) (ncons (car free)) 89 '(1 1 nil $dim nil) 90 (ncons equ))))))) 91 92(defun tabulate (e) ;For each indexed object in E, appends a list of the 93 (cond ((atom e)) ;name of that object and the number of covariant and 94 ((rpobj e) ;contravariant indices to the global list INDLIST 95 (setq indlist (cons (list (caar e) (length (cdadr e)) 96 (length (cdaddr e))) 97 indlist))) 98 ((or (eq (caar e) 'mplus) (eq (caar e) 'mtimes)) 99 (mapcar 'tabulate (cdr e))))) 100 101(defun unique (l) ;Returns a list of the unique elements of L 102 (do ((a l (cdr a)) (b)) 103 ((null a) b) 104 (cond ((not (member (car a) b :test #'equal)) 105 (setq b (cons (car a) b)))))) 106 107(defun summer1 (e) ;Applies SUMMER to the products and indexed objects in E 108 (cond ((atom e) e) 109 ((eq (caar e) 'mplus) 110 (cons (car e) (mapcar 'summer1 (cdr e)))) 111 ((or (eq (caar e) 'mtimes) (rpobj e)) 112 (summer e (cdaddr ($indices e)))) 113 (t e))) 114 115(defun summer (p dummy) ;Makes implicit sums explicit in the product or indexed 116 ;object P where DUMMY is the list of dummy indices of P 117 (prog (dummy2 scalars indexed s dummy3) ;at this level 118 (setq dummy2 (intersect (all ($indices2 p)) dummy)) 119 (do ((p (cond ((eq (caar p) 'mtimes) (cdr p)) 120 (t (ncons p))) (cdr p)) 121 (obj)) 122 ((null p)) 123 (setq obj (car p)) 124 (cond ((atom obj) 125 (setq scalars (cons obj scalars))) 126 ((rpobj obj) 127 (cond ((null (intersect dummy2 (all ($indices2 obj)))) 128 (setq scalars (cons obj scalars))) 129 (t (setq indexed (cons obj indexed))))) 130 ((eq (caar obj) 'mplus) 131 (setq s t) 132 (cond ((null (intersect dummy (all ($indices obj)))) 133 (setq scalars 134 (cons (summer1 obj) scalars))) 135 (t (setq indexed 136 (cons (summer1 obj) indexed))))) 137 (t (setq scalars (cons obj scalars))))) 138 (cond ((and s 139 (not (samelists dummy2 140 (setq s 141 (cdaddr 142 ($indices 143 (append '((mtimes)) 144 scalars indexed))))))) 145 (setq dummy3 s 146 s scalars 147 scalars nil) 148 (do ((p s (cdr p)) (obj)) 149 ((null p)) 150 (setq obj (car p)) 151 (cond ((null (intersect dummy3 (all ($indices obj)))) 152 (setq scalars (cons obj scalars))) 153 (t (setq indexed (cons obj indexed))))))) 154 (return 155 (simptimes 156 (nconc (ncons '(mtimes)) 157 scalars 158 (cond ((not (null indexed)) 159 (do ((indxd (simptimes (cons '(mtimes) indexed) 160 1 nil)) 161 (dummy (itensor-sort dummy2) (cdr dummy))) 162 ((null dummy) (ncons indxd)) 163 (setq indxd (nconc (ncons '($sum)) 164 (ncons indxd) 165 (ncons (car dummy)) 166 '(1 $dim))))) 167 (t nil))) 168 1 nil)))) 169 170(defun all (l) ;Converts [[A, B], [C, D]] into (A B C D) 171 (append (cdadr l) (cdaddr l))) 172 173(defun t-convert (e) ;Applies CHANGEFORM to each individual object in an 174 (cond ((atom e) e) ;expression 175 ((or (eq (caar e) 'mplus) (eq (caar e) 'mtimes)) 176 (cons (car e) (mapcar 't-convert (cdr e)))) 177 ((eq (caar e) '$sum) 178 (append (ncons (car e)) (ncons (t-convert (cadr e))) (cddr e))) 179 (t (changeform e)))) 180 181(defun changeform (e) ;Converts a single object from ITENSOR format to 182 (cond ((atom e) e) ;ETENSR format 183 ((rpobj e) 184 (do ((deriv (cdddr e) (cdr deriv)) 185; (new (cond ((and (null (cdadr e)) (null (cdaddr e))) 186 (new (cond ((and (null (covi e)) (null (conti e))) 187 (caar e)) ;If no covariant and contravariant 188 ;indices then make into an atom 189 (t (cons (cons (equiv-table (caar e)) '(array)) 190; (append (cdadr e) (cdaddr e))))))) 191 (append (covi e) (conti e))))))) 192 ((null deriv) new) 193 (setq new (append '(($diff)) (ncons new) 194 (ncons (cons '($ct_coords array) 195 (ncons (car deriv)))))))) 196 (t e))) 197 198(defun equiv-table (a) ;Makes appropriate name changes converting 199 (cond ((member a '($ichr1 %ichr1) :test #'eq) '$lcs) ;from ITENSOR to ETENSR 200 ((member a '($ichr2 %ichr2) :test #'eq) '$mcs) 201 (t a))) 202 203(declare-top (unspecial indlist)) 204 205(declare-top (special smlist $funcs)) 206(setq $funcs '((mlist))) 207 208(defun $makebox (e name) 209 (cond ((atom e) e) 210 ((mtimesp e) (makebox e name)) 211 ((mplusp e) 212 (mysubst0 (simplifya (cons '(mplus) 213 (mapcar 214 (function 215 (lambda (q) ($makebox q name))) 216 (cdr e))) 217 nil) 218 e)) 219 ((mexptp e) (list (car e) ($makebox (cadr e) name) (caddr e))) 220 (t e))) 221 222(defun makebox (e name) 223 (prog (l1 l2 x l3 l) 224 (setq l (cdr e)) 225 again(setq x (car l)) 226 (cond ((rpobj x) 227 (cond ((and (eq (caar x) name) (null (cdddr x)) 228 (null (cdadr x)) (= (length (cdaddr x)) 2)) 229 (setq l1 (cons x l1))) 230 ((cdddr x) (setq l2 (cons x l2))) 231 (t (setq l3 (cons x l3))))) 232 (t (setq l3 (cons x l3)))) 233 (and (setq l (cdr l)) (go again)) 234 (cond ((or (null l1) (null l2)) (return e))) 235 (do ((l2 l2 (cdr l2))) 236 ((null l2) ) 237 (setq l l1) 238; (do l2 l2 (cdr l2) 239; (null l2) 240; (setq l l1)..) 241 242 (tagbody 243 loop 244 (setq x (car l)) 245 (cond 246 ((and (member (car (cdaddr x)) (cdddar l2) :test #'eq) 247 (member (cadr (cdaddr x))(cdddar l2) :test #'eq)) 248 (setq 249 l3 250 (cons (nconc 251 (list 252 (ncons 253 (implode (append '([ ]) 254 (cdr (explodec (caaar l2)))))) 255 (cadar l2) 256 (caddar l2)) (setdiff (cdddar l2)(cdaddr x))) 257 l3)) 258 (setq l1 (delete x l1 :count 1 :test #'equal))) 259 ((setq l (cdr l)) (go loop)) 260 (t (setq l3 (cons (car l2) l3)))))) 261 (return (simptimes (cons '(mtimes) (nconc l1 l3)) 262 1. 263 nil)))) 264 265(declare-top (special tensr)) 266 267(defmfun $average n ((lambda (tensr) (simplifya (average (arg 1)) nil)) 268 (and (= n 2) (arg 2)))) 269 270(defun average (e) 271 (cond ((atom e ) e) 272 ((rpobj e) (cond ((or (not tensr) (eq (caar e) tensr)) 273 (average1 e)) 274 (t e))) 275 (t (cons (ncons (caar e)) (mapcar (function average) (cdr e)))))) 276 277(defun average1 (e) 278 (cond ((= (length (cdadr e)) 2) 279 (setq e (list '(mtimes) '((rat simp) 1 2) 280 (list '(mplus) 281 (cons (car e) 282 (cons (arev (cadr e)) (cddr e))) e)))) 283 ((= (length (cdaddr e)) 2) 284 (setq e (list '(mtimes) '((rat smp) 1 2) 285 (list '(mplus) 286 (cons (car e) 287 (cons (cadr e) 288 (cons (arev (caddr e)) 289 (cdddr e)))) e))))) 290 e) 291 292(defun arev (l) (list (car l) (caddr l) (cadr l))) 293 294(declare-top (unspecial tensr)) 295(add2lnc '(($average) $tensor) $funcs) 296 297(defun $conmetderiv (e g) 298 (cond ((not (symbolp g)) 299 (merror "Invalid metric name: ~M" g)) 300 (t (conmetderiv e g ((lambda (l) (append (cdadr l) (cdaddr l))) 301 ($indices e)))))) 302 303(defun conmetderiv (e g indexl) 304 (cond ((atom e) e) 305 ((rpobj e) 306 (cond ((and (eq (caar e) g) (null (cdadr e)) 307 (equal (length (cdaddr e)) 2) 308 (not (null (cdddr e)))) 309 (do ((e (cmdexpand (car e) (car (cdaddr e)) 310 (cadr (cdaddr e)) (cadddr e) indexl)) 311 (deriv (cddddr e) (cdr deriv))) 312 ((null deriv) e) 313 (setq e (conmetderiv ($idiff e (car deriv)) 314 g indexl)))) 315 (t e))) 316 (t (mysubst0 (cons (car e) 317 (mapcar 318 (function (lambda (q) 319 (conmetderiv q g indexl))) 320 (cdr e))) e)))) 321 322(defun cmdexpand (g i j k indexl) 323 (do ((dummy) (flag)) 324 (flag (list '(mplus simp) 325 (list '(mtimes simp) -1 326 (list g (ncons smlist) (list smlist dummy i)) 327 (list '($ichr2 simp) (list smlist dummy k) 328 (list smlist j))) 329 (list '(mtimes simp) -1 330 (list g (ncons smlist) (list smlist dummy j)) 331 (list '($ichr2 simp) (list smlist dummy k) 332 (list smlist i))))) 333 (setq dummy ($idummy)) 334 (and (not (member dummy indexl :test #'eq)) (setq flag t)))) 335 336(add2lnc '(($conmetderiv) $exp $name) $funcs) 337 338(defun $flush1deriv (e g) 339 (cond ((not (symbolp g)) 340 (merror "Invalid metric name: ~M" g)) 341 (t (flush1deriv e g)))) 342 343(defun flush1deriv (e g) 344 (cond ((atom e) e) 345 ((rpobj e) 346 (cond ((and (eq (caar e) g) (equal (length (cdddr e)) 1) 347 (or (and (equal (length (cdadr e)) 2) 348 (null (cdaddr e))) 349 (and (equal (length (cdaddr e)) 2) 350 (null (cdadr e))))) 351 0) 352 (t e))) 353 (t (subst0 (cons (ncons (caar e)) 354 (mapcar 355 (function (lambda (q) (flush1deriv q g))) 356 (cdr e))) e)))) 357 358(add2lnc '(($flush1deriv) $exp $name) $funcs) 359 360(defun $igeodesic_coords (exp g) 361 ($flush1deriv ($flush exp '$ichr2 '%ichr2) g)) 362 363(add2lnc '(($igeodesic_coords) $exp $name) $funcs) 364 365 366