/* Copyright (C) 2004 Viktor T. Toth <http://www.vttoth.com/>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of
 * the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be
 * useful, but WITHOUT ANY WARRANTY; without even the implied
 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 * PURPOSE.  See the GNU General Public License for more details.
 *
 * Supplement to itensor.lisp: implementation of frames and torsion
 *
 */

inonmet_flag:false;
iframe_bracket_form:true;
defcon(ifr,ifri,ifg);
defcon(ifg,ifg,kdelta);
SYM;		/* So that 5.9.1 knows about this as a case-insensitive symbol */

/* Helper function to get the metric tensor or return an error */
_g([l]):=if iframe_flag then apply(nounify(ifg), l)
         else if (?boundp)('imetric) then
           apply(nounify(if true then imetric),l)
         else error("Name of metric must be specified");

/* Helper functions to conditionally apply the nonmetricity and
   torsion tensors only if itorsion_flag:true */
_inm([l]):=if inonmet_flag then apply('inm,l) else 0;
_itr([l]):=if itorsion_flag then apply('itr,l) else 0;

/* Coefficient used internally when computing the rotation coefficients */
/*%icc1(l):=block([i:idummy()],'ifr([l[1]],[i])*_g([l[2],l[3]],[],i)+
 *            _inm([l[1]],[])*_g([l[2],l[3]],[])-_itr([l[1],l[2],l[3]])-
 *            'ifb([l[1],l[2]],[i])*_g([i,l[3]],[]))/2; */

/* The frame bracket */
ifb(l,[ld]):=if length(ld)>0 and rest(ld)#[] then
    apply('idiff,cons(ifb(l),rest((?putinones)(rest(ld)))))
  else if length(ld)>0 and length(ld[1])>0 then
    block([e:idummy()],
      _g([],[e,ld[1][1]])*funmake(ifb,[append(l,[e]),rest(ld[1])])
    )
  else block([e:idummy(),f:idummy()],
     if iframe_bracket_form or itorsion_flag then
       'ifr([l[2]],[e])*'ifr([l[3]],[f])*
       ('ifri([l[1],e],[],f)-'ifri([l[1],f],[],e)-
        _itr([e,f],[m])*ifri([l[1],m],[])
       )
     else 'ifri([l[1],e],[])*('ifr([l[2]],[f])*'ifr([l[3]],[e],f)-
                         'ifr([l[2]],[e],f)*'ifr([l[3]],[f]))
);

/* The connection coefficients */
icc1(l,[ld]):=if length(ld)>0 and rest(ld)#[] then
  apply('idiff,cons(icc1(l),rest((?putinones)(rest(ld)))))
  else
    (if iframe_flag then 'ifc1(l,[])
     else 'ichr1(l,if length(ld)>0 then ld[1] else []))+
    (if itorsion_flag and not iframe_flag then -'ikt1(l,[]) else 0)+
    (if inonmet_flag then -'inmc1(l,[]) else 0);
icc2(l1,l2,[ld]):=
  if ld#[] then apply('idiff,cons(icc2(l1,l2),rest((?putinones)(ld))))
/*else block([d:idummy()],_g([],[l2[1],d])*(%icc1([l1[1],d,l1[2]])-
                           %icc1([d,l1[2],l1[1]])+%icc1([l1[2],l1[1],d]))/2);*/
  else
    (if iframe_flag then 'ifc2(l1,l2) else 'ichr2(l1,l2))+
    (if itorsion_flag and not iframe_flag then -'ikt2(l1,l2) else 0)+
    (if inonmet_flag then -'inmc2(l1,l2) else 0);

/* The frame coefficients */
ifc1(l,[ld]):=if length(ld)>0 and rest(ld)#[] then
  apply('idiff,cons(ifc1(l),rest((?putinones)(rest(ld)))))
  else ('ifb(l)+'ifb([l[2],l[3],l[1]])-'ifb([l[3],l[1],l[2]]))/2;
ifc2(l1,l2,[ld]):=if length(ld)>0 then
    apply('idiff,cons(ifc2(l1,l2),rest((?putinones)(ld))))
  else block([d:idummy()],_g([],[l2[1],d])*'ifc1([l1[1],l1[2],d]));


/* The nonmetricity coefficients */
inmc1(l,[ld]):=if not inonmet_flag then 0
  else if length(ld)>0 and rest(ld)#[] then
  apply('idiff,cons(inmc1(l),rest((?putinones)(rest(ld)))))
  else (-_inm([l[1]])*_g([l[2],l[3]])-_inm([l[2]])*_g([l[1],l[3]])+
        _inm([l[3]])*_g([l[1],l[2]]))/2;
inmc2(l1,l2,[ld]):=if not inonmet_flag then 0
  else if ld#[] then apply('idiff,cons(inmc2(l1,l2),rest((?putinones)(ld))))
  else block([m:idummy()],(-_inm([l1[1]])*'kdelta([l1[2]],[l2[1]])-
                           _inm([l1[2]])*'kdelta([l1[1]],[l2[1]])+
                           _g([],[l2[1],m])*_inm([m])*_g([l1[1],l1[2]]))/2);

/* Contortion */
ikt1(l,[ld]):=if not itorsion_flag then 0
  else if length(ld)>0 and rest(ld)#[] then
  apply('idiff,cons(ikt1(l),rest((?putinones)(rest(ld)))))
  else block([d:idummy()],(-_g([l[3],d])*_itr([l[1],l[2]],[d])-_g([l[2],d])*
                _itr([l[3],l[1]],[d])-_g([l[1],d])*_itr([l[3],l[2]],[d]))/2);
ikt2(l1,l2,[ld]):=if not itorsion_flag then 0
  else if ld#[] then apply('idiff,cons(ikt2(l1,l2),rest((?putinones)(ld))))
  else block([e:idummy()],_g([],[l2[1],e])*'ikt1([l1[1],l1[2],e]));

/* Simplify expressions containing the metric tensor's derivatives */
/* v1
simpmetderiv(exp):=
(
    if atom(exp) then exp
    else if op(exp)="-" then -simpmetderiv(-exp)
    else if op(exp)="+" then funmake("+", map(simpmetderiv, args(exp)))
    else if op(exp)="/" then
        simpmetderiv(part(exp,1))/simpmetderiv(part(exp,2))
    else if op(exp)="*" then
    block([sign:1,args:args(exp)],
        for i thru length(args) do
            for j thru length(args) do
        (
            if i#j and ?rpobj(args[i]) and ?rpobj(args[j]) and
                   op(args[i])=imetric and op(args[j])=imetric then
            block(
                [a:if length(covi(args[i]))>0 then args[i] else args[j],
                 b:if length(covi(args[i]))>0 then args[j] else args[i]],
                if length(covi(a)) = 2 and length(conti(a)) = 0 and
                   length(covi(b)) = 0 and length(conti(b)) = 2 and
                   length(?intersect(covi(a),conti(b))) = 1 then
                (
                    if (flipflag and length(deri(a)) = 1 and
                                     length(deri(b)) = 0) or
                       (not flipflag and length(deri(a)) = 0 and
                                         length(deri(b)) = 1) then
                    block(
                        [tmp:deri(a)],
                        args[i]:funmake(op(a),
                                        append([covi(a),conti(a)],deri(b))),
                        args[j]:funmake(op(b),append([covi(b),conti(b)],tmp)),
                        sign:-sign
                    )
                )
            )
        ),
        sign*funmake("*",args)
    )
    else exp
); */

simpmetderiv(exp,[stop]):=
(
    if atom(exp) then exp
    else if op(exp)="-" then -apply(simpmetderiv,cons(-exp,stop))
    else if op(exp)="+" then
      funmake("+", map(lambda([x],apply(simpmetderiv,cons(x,stop))), args(exp)))
    else if op(exp)="/" then apply(simpmetderiv,cons(part(exp,1),stop))/
                              apply(simpmetderiv,cons(part(exp,2),stop))
    else if op(exp)="*" then
    block([sign:1,args:args(exp)],
        for i thru length(args) do
            for j thru length(args) do
        (
            if i#j and ?rpobj(args[i]) and ?rpobj(args[j]) and
                   op(args[i])=imetric and op(args[j])=imetric then
            block(
                [a:if length(covi(args[i]))>0 then args[i] else args[j],
                 b:if length(covi(args[i]))>0 then args[j] else args[i]],
                if length(covi(a)) = 2 and length(conti(a)) = 0 and
                   length(covi(b)) = 0 and length(conti(b)) = 2 and
                    (
                        (
                            sort(covi(a)) = sort(conti(b)) and
                            length(deri(a)) = 1 and length(deri(b)) = 1 and
                            (
                                (flipflag and
                                 ordergreatp(deri(a)[1], deri(b)[1])) or
                                (not flipflag and
                                 ordergreatp(deri(b)[1], deri(a)[1]))
                            )
                        ) or
                        (
                            length(covi(a)) = 2 and length(conti(a)) = 0 and
                            length(covi(b)) = 0 and length(conti(b)) = 2 and
                            length(?intersect(covi(a),conti(b))) >= 1 and
                            (
                                (flipflag and length(deri(a)) = 1 and
                                 length(deri(b)) = 0) or
                                (not flipflag and length(deri(a)) = 0 and
                                 length(deri(b)) = 1)
                            ) and (sign:-sign) # 0
                        )
                    ) then

                block(
                    [tmp:deri(a)],
                    args[i]:funmake(op(a),
                                    append([covi(a),conti(a)],deri(b))),
                    args[j]:funmake(op(b),append([covi(b),conti(b)],tmp)),
                    if stop#[] then i:j:length(args)
                )
            )
        ),
        sign*funmake("*",args)
    )
    else exp
);


/* Always true symmetries */
decsym(ichr1,3,0,[sym(1,2)],[]);
decsym(ichr2,2,1,[sym(all)],[]);
decsym(icurvature,3,1,[anti(2,3)],[]);
/* decsym(ifb,3,0,[anti(2,3)],[]);	<-- not valid with torsion
 * decsym(icc1,3,0,[sym(1,2)],[]);
 * decsym(icc2,2,1,[sym(all)],[]);
 * decsym(ifc1,3,0,[sym(1,2)],[]);
 * decsym(ifc2,2,1,[sym(all)],[]);
 * decsym(ikt1,3,0,[sym(1,2)],[]);
 * decsym(ikt2,2,1,[sym(all)],[]);*/