?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Use %%PRETTY-FNAME in more quadpack error messages
[maxima.git] / share / tensor / gener.lisp
blob1e99d0f1f05435db8046501b4923dd9ffbd159c9
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 "~
62 IC_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 \f
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 \f
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 \f
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 \f
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 \f
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
Maxima, a Computer Algebra System