/************************************************************************ July 29, 2019 dgmatrix3.mac is a package of Maxima functions which contains code for explicit Dirac spinor and matrix calculations and is part of the Dirac package. Maxima by Example Ch. 12, Dirac Algebra and Quantum Electrodynamics version = dirac3 Copyright (C) 2010, 2011,2017, 2018,2019 Edwin L. Woollett http://www.csulb.edu/~woollett This program is free software: you can redistribute it and/or modify it under the terms of the GNU GENERAL PUBLIC LICENSE, Version 2, June 1991, as published by the Free Software Foundation. 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. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.fsf.org/licensing/. ************************************************************************/ /* file dgmatrix3.mac for use with symbolic package dirac3.mac use comp_def to define aray components of four vectors (comp_def is defined in dgeval.mac) dgmatrix3 functions in file order: ----------------------------------- eps3 completely antisymmetric tensor in 3 dimensions Cross(j,La,Lb) j'th component of pair of vectors represented as lists. vprod (La, Lb) list of 3 components of the vector product of lists La and Lb cdouble, mdouble Sig[j],Sigb[j],I2,Z2,Cz2,sL2,sLb2 I4,Z4,CZ4,RZ4 ipm Avsq hc comm acomm Gam[mu] sL(p) P(sv,Sp) or P(sv) m_tr sbar mbar Chi UU VV mnctimesp, mnctimesp(a . b) -> true mncexptp, mncexptp (a^^2 ) --> true mdivp mdivp (a/x) --> true acon mcon mp_split rootcrunch trigcrunch */ /* to remove all global arrays currently existing, could use: if length (arrays) > 0 then apply ('remarray, arrays)$ */ /********* block matrix tools cdouble and mdouble *************/ /**** cdouble produces a 2-element column vector from two numbers, or produces a 4-element column vector from two 2-element column vectors *****************/ cdouble (ma,mb) := mat_unblocker ( matrix ([ma],[mb]) )$ /******** mdouble constructs 4 by 4 matrix from four 2 by 2 matrices ****/ mdouble(Ma,Mb,Mc,Md) := mat_unblocker (matrix ([Ma,Mb], [Mc,Md]) )$ /* pauli 2 x 2 matrices Sig[j] and Sigb[j] for j = 0,1,2,3 ref: Maggiore (M) 2.47, 2.67 */ /**************************************/ ( I2 : ident(2), Z2 : zeromatrix (2,2), CZ2 : matrix ( [0],[0] ), Sig[0] : ident (2), Sig[1] : matrix ([0,1],[1,0]), Sig[2] : matrix ([0,-%i],[%i,0]), Sig[3] : matrix ([1,0],[0,-1]), Sigb[0] : ident (2), for j thru 3 do Sigb[j] : - Sig[j], sL2(_q%) := sum ( _q%[jj]*Sig[jj],jj,1,3 ), sL2b(_q%) := sum ( _q%[jj]*Sigb[jj],jj,1,3 ) )$ /* Four component tools */ /* 4 x 4 unit matrix */ I4 : ident (4)$ /* 4 x 4 zero matrix */ Z4 : zeromatrix(4,4)$ /* 4 element column vector of zeros */ CZ4 : cdouble ( CZ2, CZ2)$ /* 4 element row vector of zeros */ RZ4 : matrix ( [0,0,0,0])$ /* print(151)$ */ /************ ipm(M1,M2) *******************************/ /* Lorentz invariant inner product of two matrices which each carry a single Lorentz index --> four vector dot product of two four vectors, each represented by a list of elements. */ ipm (_mp%, _mq%) := _mp%[0] . _mq%[0] - _mp%[1] . _mq%[1] - _mp%[2] . _mq%[2] - _mp%[3] . _mq%[3]$ /* print(164)$ */ /******************** misc tools for matrices *********************************/ /********** absolute value squared ********************/ Avsq(_e%) := expand (_e% * conj(_e%))$ /******** hermitian conjugate of a matrix ***********/ /* hc (_MM%) := transpose ( conj (_MM%))$ */ hc (_MM%) := (conj (_MM%), transpose (%%))$ /*********** commutator of two matrices *************/ comm (_M1%,_M2%) := _M1% . _M2% - _M2% . _M1% $ /* print(185)$ */ /* 3-vector utilities */ /* 3d completely antisymmetric unit tensor */ declare ( eps3, antisymmetric )$ tellsimpafter(eps3[1,2,3],1)$ /* (%i10) eps3[1,2,3]; (%o10) 1 (%i11) eps3[2,1,3]; (%o11) -1 (%i12) eps3[1,1,3]; (%o12) 0 */ /* 3 vector cross product of two vectors represented by lists */ Cross(ii,A,B) := sum (sum (eps3[ii,jj,kk]*A[jj]*B[kk],kk,1,3),jj,1,3)$ /* (%i14) p : [2,-2,1]; (%o14) [2,-2,1] (%i15) q : [-3,4,0]; (%o15) [-3,4,0] (%i16) p . p; (%o16) 9 (%i17) q . q; (%o17) 25 (%i18) q . p; (%o18) -14 (%i19) for i thru 3 do print(i," ",Cross(i,p,q))$ 1 -4 2 -3 3 2 */ /* 3-vector cross product of pair of vectors defined as lists; returns a list */ vprod (La, Lb) := [Cross(1,La,Lb), Cross(2,La,Lb), Cross(3,La,Lb)]$ /* MCross(i,Amatrix,Bmatrix) works with the 4 x 4 matrices Gam[j], j = 1,2 ,3 defined above */ MCross(ii,Am,Bm) := sum (sum (eps3[ii,jj,kk]*Am[jj] . Bm[kk],kk,1,3),jj,1,3)$ /* (%i24) %i*MCross(1,Gam,Gam)/2; (%o24) matrix([0,1,0,0],[1,0,0,0],[0,0,0,1],[0,0,1,0]) (%i25) %i*MCross(2,Gam,Gam)/2; (%o25) matrix([0,-%i,0,0],[%i,0,0,0],[0,0,0,-%i],[0,0,%i,0]) (%i26) %i*MCross(3,Gam,Gam)/2; (%o26) matrix([1,0,0,0],[0,-1,0,0],[0,0,1,0],[0,0,0,-1]) */ /* function which writes the expansion (a sum over sig[n3]) equivalent to the commutator comm ( Sig[n1], Sig[n2] ) making use of eps3 Sig_comm (n1,n2) := if eps3[n1,n2,nn] = 0 then Z2 else ( 2*%i* sum ( eps3 [n1,n2,n3]*Sig[n3],n3,1,3), ev(%%) )$ */ /*********** anticommutator of two matrices ************/ acomm (_M1%,_M2%) := _M1% . _M2% + _M2% . _M1% $ /* 4 x 4 matrices constructed from 2 x 2 Pauli matrices and Z2 First, definition of SIG[j] for j = 1,2,3, equivalent to the definition SIG-vec = (%i/2) * Gam-vec x Gam-vec (3 vec cross product), The helicity operator (in matrix form) for a Dirac particle is (1/2) times the three vector dot product of SIG-vec and p-hat-vec, where the latter is p-vec/magnitude(p-vec). */ for j thru 3 do SIG[j] : mdouble (Sig[j],Z2,Z2,Sig[j])$ /* print(285)$ */ /*********** Gam[mu] ***********************************************/ /* construct 4 x 4 gamma matrices in terms of the Pauli matrices Sig[j] and Z2 and I2 (P/S convention chiral rep) */ /*****************************************************/ ( /* Chiral def of dirac gamma matrices Gam[j], Peskin-Schroeder convention, also Majjiore */ Gam[0] : mdouble ( Z2,I2,I2,Z2 ), for j thru 3 do Gam[j] : mdouble (Z2,Sig[j],-Sig[j],Z2 ), Gam[5] : %i*Gam[0] . Gam[1] . Gam[2] . Gam[3] )$ /*********** four vector and spinor utilities *************/ /* to use Feynman slash sL and/or helicity projection matrices P, need to define four vector components as arrays via the funtion comp_def defined in dgeval3.mac */ /*********** Feynman slash matrix sL(p) and sL1(p) where p is an array or sum of arrays in sL and a single array in sL1; can't use declare(sL, linear) when args are arrays. *****************/ sL1 (_r%) := ipv (_r%,Gam)$ /* new sL function Dec. 12, 2018 */ sL (_e%) := block ([Mret : Z4, _ep%, badarg : false], if atom(_e%) then return (sL1(_e%) ) else if op (_e%) = "-" then return ( - sL1(- _e%)) else for j thru length (_e%) do ( _ep% : part(_e%, j), if atom (_ep%) then Mret : Mret + sL1(_ep%) else if op (_ep%) = "-" then Mret : Mret - sL1(- _ep%) else ( print ("cannot have constant multiplying array in sL"), badarg : true, return())), if badarg then Z4 else Mret)$ /***** explicit helicity projection matrix P(sv,Sp) or P(sv) with sv = -1 or 1, Sp1 is the right-handed spin 4-vector generated by a four vector p1 with components (E,px,py,pz) = (E,p3_vec) = (E,p * p_hat) in which case Sp1 has components (p/m, (E/m) * p_hat ) P(sv) --> I4 + sv*Gam[5]/2 , m = 0 case P(sv, Sp) --> I4 + sv*Gam[5] . sL (Sp)/2 , m > 0 case */ P([v]) := block ([lv], lv : length(v), if lv = 1 then (I4 + v[1]*Gam[5])/2 else if lv = 2 then (I4 + v[1]*Gam[5].sL(v[2]))/2 else ( disp ("P(sv) or P(sv,Sp); sv = 1 or -1"), Z4 ) )$ /* non commuting times a . b (%i12) ?mnctimesp(a . b); (%o12) true */ mnctimesp(e) := ?mnctimesp(e)$ /* non commuting power a^^2 = a . a (%i13) ?mncexptp(a^^2); (%o13) true */ mncexptp(e) := ?mncexptp(e)$ /* note also ?mtimesp behavior (%i4) ?mtimesp(a/b); (%o4) true (%i5) ?mtimesp(a*b); (%o5) true */ mdivp(e) := ?mtimesp(e) and op(e) = "/"$ /* (%i6) op(a/b); (%o6) "/" (%i8) mdivp(a*b); (%o8) false (%i9) mdivp(a/b); (%o9) true */ /* mp_split(asum) 6-10-10 returns the list [ pptL, mptL] , where pptL is a list of the part numbers of the multiple momenta, and mptL is a list with the single part number of the mass term, called by m_tr */ mp_split(asum) := block ([jj,pptL:[],mptL:[],pp ], for jj thru length (asum) do ( pp : part (asum,jj), if (not atom(pp) and op(pp) = "-") then pp : -pp, if (not atom(pp) and op(pp) = "*") then pptL : cons (jj,pptL), if (atom(pp) and massp(pp)) then mptL : cons (jj,mptL), if (atom(pp) and not massp(pp)) then pptL : cons (jj,pptL)), [pptL,mptL] )$ /* examples: (%i38) sm : c1*p1 + p2 +m$ (%i39) mp_split(sm); (%o39) [[2,1],[3]] (%i40) part(sm,2); (%o40) c1*p1 (%i41) part(sm,1); (%o41) p2 (%i42) part(sm,3); (%o42) m (%i43) sm : c1*p1 + p2 - m$ (%i44) mp_split(sm); (%o44) [[2,1],[3]] (%i45) part(sm,3); (%o45) -m */ /********** m_tr 2-28-11, 6-12-10, 4-21-10 , sept 2018 ***************************/ /* ability to handle multiple momenta args 6-12-10 (%i12) a . Gam[5] . 1; (%o12) a . Gam[5] */ /* syntax m_tr(p,p+m,q-m,S(1),mu,S(-1,Sp),nu,p1+p2+M,G5) where mu and nu (if meant as Lorentz indices) are either integers in the range 0 thru 3 or symbols on the list indexL. G5 represents Gam[5]. */ m_tr ([_v%]) := block ([_vv%,adummy,xx,_PP%,fac1,_PPA%,_pp1%,_pp2%,badarg,_mprod%, psL,mt,jj,kk], /* disp ("m_tr 2-28-11 "), */ /* display (_v%), */ if length (_v%) = 0 then return(1), if length (_v%) = 1 then if _v%[1] = 1 then return(4), /* we want to use the jj do loop number in returned error diagnostics */ /* Aug 22, 2018, reverse arg list */ _vv% : cons (adummy, reverse(_v%) ), /* build up a return product value */ _mprod% : 1, for jj thru length (_vv%) do ( xx : part (_vv%, jj), /* display (jj, xx), */ if xx # adummy then ( if atom (xx) then ( /* disp (" case atom(xx) "), */ if xx = G5 then _mprod% : Gam[5] . _mprod% else if (integerp(xx) or indexp(xx)) then ( /* disp ("case indexp = true "), */ _mprod% : Gam[xx] . _mprod% /*, display (xx,_mprod%) */ ) else ( /* disp ("case not indexp "), */ _mprod% : sL(xx) . _mprod% /* ,display (_mprod%) */ ) ) /* end case atom */ else if op(xx) = S then ( /* disp (" case op = S "), */ _PP% : apply ('P,args(xx)), _mprod% : _PP% . _mprod% /*,display (_mprod%) */ ) /* end case S */ /* case op(xx) = "+" */ else if mplusp (xx) then ( /* disp ("case op(xx) is plus "), display (xx), */ /* cases p + m, p - m ,p1+p2 +m,p1-p2-m, etc */ [psL,mt] : mp_split (xx), psum : 0, badarg : false, for kk thru length (psL) do ( _PPA% : part (xx,psL[kk]), [fac1,_PP%] : NDfac (_PPA%), if atom (_PP%) then psum : psum + fac1*sL (_PP%) /* case mtimesp */ else if mtimesp (_PP%) then ( _pp1% : part (_PP%,1), _pp2% : part (_PP%,2), if scalarp (_pp1%) then psum : psum + fac1* _pp1%*sL(_pp2%) else if scalarp (_pp2%) then psum : psum + fac1* _pp2%*sL(_pp1%) else ( disp ("m_tr: bad arg "), badarg : true, return ())) else ( disp ("m_tr: bad arg "), badarg : true, return ())), /* end kk do loop */ if badarg then return (), if length (mt) = 0 then _mprod% : (psum) . _mprod% else ( psum : psum + part (xx,mt[1])*I4, _mprod% : (psum) . _mprod%)) /* end mplusp(xx) case */ else ( print (" m_tr: property not found for jj = ",jj), badarg : true, return ()))), /* end construct jj contrib. to _mprod% */ if badarg then return (apply ('M_TR,_v%)), /* display (_mprod%), */ expand (mat_trace (_mprod%)) )$ /********** end m_tr 6-12-10, 4-21-10, sept 2018 ***************************/ /**** sbar: make barred dirac spinor ***************/ /* sbar (_uu%) := hc (_uu%) . Gam[0]$ */ /* sbar (_uu%) := block ([temp], temp : conj (_uu%), temp : transpose (_uu%), temp : temp . Gam[0], temp)$ */ sbar (_uu%) := ( conj (_uu%), transpose (%%), %% . Gam[0] )$ /******* barred dirac matrix mbar ***************/ mbar(_M%) := Gam[0] . hc (_M%) . Gam[0] $ /****** dirac four component spinors ********/ /** cdouble, used below, and defined above, produces a 2-element column vector from two numbers and produces a 4-element column vector from two 2-element column vectors. ****/ /** first define two component helicity eigenspinors Chi **/ Chi (_th%,_phi%,_sv%) := ( if _sv% = 1 then cdouble (cos (_th%/2),sin (_th%/2)*exp(%i*_phi%) ) else if _sv% = -1 then cdouble( -sin (_th%/2)*exp(-%i*_phi%), cos (_th%/2) ) else ( disp ("Chi: helicity quantum number sv can only be -1 or 1 "), CZ2))$ /* print(595)$ */ /* Four component particle and antiparticle dirac spinors with helicity sv = 1 or -1, energy E, three momentum magnitude p and spherical-polar angles th and phi. The range of the angle th (the angle between the 3-momentum vector and the positive z axis) is 0 to %pi radians, and the range of the azimuthal angle phi is 0 to 2*%pi radians */ /* old UU UU( _EE%,_pp%,_tth%,_pphi%,_ssv%) := (cdouble ( sqrt (_EE% - _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%), sqrt (_EE% + _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%) ), expand( matrixmap (rootscontract,%%) ) )$ */ /* new UU July 11, 2018, has radexpand locally set to false UU( _EE%,_pp%,_tth%,_pphi%,_ssv%) := block([radexpand:false], cdouble ( sqrt (_EE% - _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%), sqrt (_EE% + _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%) ))$ */ /* print(615)$ */ /* UU Sept. 21, 2018 UU( _EE%,_pp%,_tth%,_pphi%,_ssv%) := block([radexpand:false, halfangles:false], cdouble ( sqrt (_EE% - _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%), sqrt (_EE% + _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%) ), matrixmap (lambda ([x],ratsubst (cos (_tth%/2), sin ( (%pi - _tth%)/2 ), x)), %%), matrixmap (lambda ([x],ratsubst (sin (_tth%/2), cos ( (%pi - _tth%)/2 ), x)), %%), matrixmap (lambda ([x],ratsubst (sin (_tth%/2), cos ( (_tth% - %pi )/2 ), x)), %%), matrixmap (lambda ([x],ratsubst (- cos (_tth%/2), sin ( (_tth% - %pi )/2 ), x)), %%), rootscontract (%% ) )$ */ sub_all (acol, angle) := (list_matrix_entries (acol), map (lambda ([x], ratsubst (sin(angle/2), cos((angle - %pi)/2), x)), %%), map (lambda ([x], ratsubst (- cos(angle/2), sin((angle - %pi)/2), x)), %%), map (lambda ([x], ratsubst (cos(angle/2), sin((%pi - angle)/2), x)), %%), map (lambda ([x], ratsubst (sin(angle/2), cos((%pi - angle)/2), x)), %%), makelist ( [%%[i]], i, 1, length(acol)), apply ('matrix, %%))$ uu1( energy,mommag,pangle,azangle,helicity) := (cdouble ( sqrt (energy - helicity * mommag)*Chi (pangle,azangle,helicity), sqrt (energy + helicity * mommag)*Chi (pangle,azangle,helicity) ), rootscontract (%%))$ /* UU sept 22, 2018 UU( _EE%,_pp%,_tth%,_pphi%,_ssv%) := (uu1( _EE%,_pp%,_tth%,_pphi%,_ssv%), sub_all (%%, _tth%))$ */ /* print(662)$ */ /* UU sept 23, 2018 */ UU( _EE%,_pp%,_tth%,_pphi%,_ssv%) := (uu1( _EE%,_pp%,_tth%,_pphi%,_ssv%), subst ([cos((%pi - th)/2) = sin (th/2), sin((%pi - th)/2) = cos (th/2), cos((th - %pi)/2) = sin(th/2), sin((th - %pi)/2) = - cos(th/2)],%%))$ /* (cdouble ( sqrt (_EE% - _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%), sqrt (_EE% + _ssv% * _pp%)*Chi (_tth%,_pphi%,_ssv%) ), */ /* print(659)$ */ VV(_E%,_p%,_th%,_phi%,_sv%) := -Gam[5] . UU (_E%,_p%,_th%,_phi%, - _sv% )$ /* print(665)$ */ /* vv1(_E%,_p%,_th%,_phi%,_sv%) := -Gam[5] . uu1 (_E%,_p%,_th%,_phi%, - _sv% )$ */ /* VV(_E%,_p%,_th%,_phi%,_sv%) := block([radexpand:false, halfangles:false, l0, l1,l2,l3,l4,l5], l0 : UU (_E%,_p%,_th%,_phi%, - _sv% ), print ( " l0 = ",l0), l1 : -Gam[5] . l0, print ( " l1 = ",l1), l2 : matrixmap (lambda([x], ratsubst (cos (_th%/2), sin ( (%pi - _th%)/2 ), x)),l1), print ( " l2 = ", l2), l3 : matrixmap (lambda([x], ratsubst (sin (_th%/2), sin ( (%pi - _th%)/2 ), x)), l2), print ( " l3 = ", l3), l4 : matrixmap (lambda([x], ratsubst (sin (_th%/2), cos ( (_th% - %pi )/2 ), x)),l3), print ( " l4 = ", l4 ), l5 : matrixmap (lambda([x], ratsubst (- cos (_th%/2), sin ( (_th% - %pi )/2 ), x)),l4), print ( " l5 = ", l5), l5)$ */ /* print(688)$ */ /* code for contraction of a product of two antisymmetric tensors transferred from coneps2.mac */ /* kronecker delta: use core function: kron_delta(a,b) (%i31) kron_delta(0,0); (%o31) 1 (%i32) kron_delta(0,1); (%o32) 0 core properties of kron_delta: (%i1) constantp(kron_delta(a,b)); (%o1) false (%i2) constantp(kron_delta(a,a)); (%o2) true (%i3) scalarp(kron_delta(a,b)); (%o3) true (%i4) scalarp(kron_delta(a,a)); (%o4) true (%i5) kron_delta(n,n); (%o5) 1 (%i6) kron_delta(n1,n1); (%o6) 1 (%i7) kron_delta(a,a); (%o7) 1 */ /* following assumes the product eps4[...]*eps4L[...] is in canonical form: repeated indices are in first slots in both factors and in the same slots in both factors: if mu is a repeated index, it is in the same slot for both factors. */ /* all four indices summed over [0,3] */ eps4_prod4 : -24$ /* first three indices summed over [0,3] */ eps4_prod3(%n4,%m4) := -6*kron_delta(%n4,%m4)$ /* first two indices summed over */ eps4_prod2(%n3,%n4,%m3,%m4 ) := -2*kron_delta(%n3,%m3)*kron_delta(%n4,%m4) + 2*kron_delta(%n3,%m4)*kron_delta(%n4,%m3)$ /* first index summed over */ eps4_prod1 (%n2,%n3,%n4,%m2,%m3,%m4 ) := -kron_delta(%n2,%m2)*kron_delta(%n3,%m3)*kron_delta(%n4,%m4) + kron_delta(%n2,%m3)*kron_delta(%n3,%m2)*kron_delta(%n4,%m4) + kron_delta(%n2,%m2)*kron_delta(%n3,%m4)*kron_delta(%n4,%m3) - kron_delta(%n2,%m3)*kron_delta(%n3,%m4)*kron_delta(%n4,%m2) - kron_delta(%n2,%m4)*kron_delta(%n3,%m2)*kron_delta(%n4,%m3) + kron_delta(%n2,%m4)*kron_delta(%n3,%m3)*kron_delta(%n4,%m2)$ /* examples of syntax: one index contraction case: coneps2 (eps4[n1,n2,n3,n4]*eps4L[n1,m2,m3,m4], n1) two indices contraction case: coneps2 (eps4[n1,n2,n3,n4]*eps4L[n1,n2,m3,m4], n1,n2) syntax: coneps2( an_eps_product,n1,n2,....) contraction on repeated indices supplied. only allowed eps product is eps4*eps4L or eps4L*eps4 */ coneps2 (%epsprod,[%vL]) := block ([pp, %fac:1, Nr : 0, Nrepeat,aL1, aL2,vv1,rL:[], op1,op2, pp1,pp2,ff1,ff2,LL1,LL2,LL12 ], pp : %epsprod, Nrepeat : length(%vL), if Nrepeat = 0 then (print (" coneps2: no contraction indices provided "), return (done)), /* print (" pp = ",pp," %vL = ", %vL), */ /* print (" %fac = ",%fac," Nrepeat = ",Nrepeat), */ if op (pp) = "-" then (%fac : -1, pp : first (pp)), /* print (" %fac = ",%fac," pp = ",pp ), */ if length (pp) # 2 then (print(" coneps2: not just a product of two eps factors"), return (done)), op1 : op(first(pp)), op2 : op(second(pp)), /* print(" op1 = ",op1," op2 = ",op2), */ if op1 = op2 then (print("coneps2: syntax: coneps(eps4*eps4L,n1,..) or coneps(eps4L*eps4,n1,..)"), return (done)), if not inlist ([eps4, eps4L], op1) then (print("coneps2: syntax: coneps(eps4*eps4L,n1,..) or coneps(eps4L*eps4,n1,..)"), return (done)), if not inlist ([eps4, eps4L], op2) then (print("coneps2: syntax: coneps(eps4*eps4L,n1,..) or coneps(eps4L*eps4,n1,..)"), return (done)), aL1 : args(first(pp)), aL2 : args(second(pp)), /* print(" aL1 = ",aL1," aL2 = ",aL2), */ /* case: four index summation, in canonical order */ if Nrepeat=4 and aL1 = aL2 then return (%fac*eps4_prod4 ), /* Nr is the number of repeated indices found. rL is a list of repeated indices found in the order of appearance in the list aL1 */ for vv1 in aL1 do if not numberp(vv1) and inlist(aL2,vv1) then (Nr : Nr+1, rL : cons (vv1, rL)), rL : reverse(rL), /* print ( " Nr = ",Nr," rL = ",rL), */ if Nr # Nrepeat then (print( " coneps2: Nr # Nrepeat "), return(done)), for j thru Nrepeat do ( vv1 : rL[j], pp1 : pos (aL1,vv1), pp2 : pos (aL2, vv1), [ff1,aL1] : move (aL1,pp1,j), [ff2,aL2] : move (aL2,pp2,j), %fac : %fac*ff1*ff2), LL1 : rest(aL1,Nrepeat), LL2 : rest(aL2,Nrepeat), LL12 : flatten(cons (LL1,LL2)), /* print(" %fac = ",%fac," LL12 = ",LL12), */ if Nrepeat = 1 then %fac*apply('eps4_prod1,LL12) else if Nrepeat = 2 then %fac*apply('eps4_prod2,LL12) else %fac*apply('eps4_prod3,LL12))$ /* Nrepeat = 4 cases (%i2) coneps2 (eps4[n1,n2,n3,n4]*eps4L[n1,n2,n3,n4],n1,n2,n3,n4); (%o2) -24 (%i3) coneps2 (eps4[n2,n1,n3,n4]*eps4L[n1,n2,n3,n4],n1,n2,n3,n4); (%o3) 24 (%i4) coneps2 (eps4[n2,n3,n1,n4]*eps4L[n1,n2,n3,n4],n1,n2,n3,n4); (%o4) -24 (%i5) coneps2 (eps4[n2,n3,n4,n1]*eps4L[n1,n2,n3,n4],n1,n2,n3,n4); (%o5) 24 */ /* Nrepeat = 3 cases (%i6) coneps2 (eps4[n1,n2,n3,n4]*eps4L[n1,n2,n3,m4],n1,n2,n3); (%o6) -6*kron_delta(m4,n4) (%i7) coneps2 (eps4[n2,n1,n3,n4]*eps4L[n1,n2,n3,m4],n1,n2,n3); (%o7) 6*kron_delta(m4,n4) (%i8) coneps2 (eps4[n2,n3,n1,n4]*eps4L[n1,n2,n3,m4],n1,n2,n3); (%o8) -6*kron_delta(m4,n4) (%i13) coneps2 (eps4[n2,n3,n1,n4]*eps4[n1,n2,n3,m4],n1,n2,n3); coneps2: syntax: coneps(eps4*eps4L,n1,..) or coneps(eps4L*eps4,n1,..) (%o13) done */ /* Nrepeat = 2 cases (%i9) coneps2 (eps4[n1,n2,n3,n4]*eps4L[n1,n2,m3,m4],n1,n2); (%o9) 2*kron_delta(m3,n4)*kron_delta(m4,n3)-2*kron_delta(m3,n3) *kron_delta(m4,n4) (%i10) coneps2 (eps4[n2,n1,n3,n4]*eps4L[n1,n2,m3,m4],n1,n2); (%o10) 2*kron_delta(m3,n3)*kron_delta(m4,n4)-2*kron_delta(m3,n4) *kron_delta(m4,n3) */ /* following two functions assume signs and/or numerical factors have been removed from precursor of Aprod using NDfac */ strip_atoms(Aprod) := listToProd (prodToList (Aprod))$ strip_non_atoms(Aprod) := (strip_ops2(Aprod, eps4)[2], strip_ops2 (%%,eps4L)[2])$ /* (%i4) aprod : f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]; (%o4) f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i8) strip_non_atoms(aprod); (%o8) f1*f2 (%i10) strip_atoms(aprod); (%o10) eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i15) [num_fac, aprod] : NDfac(-2*f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]/3); (%o15) [-2/3,f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]] (%i16) num_fac; (%o16) -2/3 (%i17) aprod; (%o17) f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] */ /* new acon, acon1 and convertKD 6-13-2017 */ /* convertKD replaces KD with core symbol kron_delta */ convertKD(expr) := subst(kron_delta,KD,expr)$ acon1(myexpr,myindex) := (subst (0,myindex,myexpr) - subst (1,myindex,myexpr) - subst (2,myindex,myexpr) - subst (3,myindex,myexpr), ev (%%),expand(%%))$ /* acon 6-13-2017 allows dealing with possible symbolic kronecker delta symbol KD(n1,n2) in contractions of one term with respect to one index. */ acon(_ee%,_mu%) := block ([%fac,%rterm,%ee,KL, KargsL,Kargs,kargsL1,kde,sarg,Kreturn ], /* %fac is numerical factor which can include a minus sign */ /* %rterm is the rest of the expression */ if _ee% = 0 then return(0), %ee : _ee%, /* print (" %ee = ",%ee), */ [%fac,%rterm] : NDfac (%ee), /* display (%fac,%rterm), */ if atom (%rterm) then (print(" acon: %rterm is an atom"), print(" anerror "), return(false)), [ KL,%rterm] : strip_ops2 (%rterm,KD), /* print (" KL = ",KL," %rterm = ",%rterm), */ if length (KL) > 0 then (KargsL : map ('args,KL), /* print (" KargsL = ",KargsL), */ if length(KL) = 1 then ( /* print (" case term contains one factor of KD"), */ Kargs : KargsL[1], /* print (" Kargs = ",Kargs), */ if inlist (Kargs, _mu%) then ( /* print (" case contraction index contained in Kargs"), */ if length (delete (_mu%, Kargs)) = 0 then return (4*%fac*%rterm), sarg : part (delete (_mu%, Kargs),1), /* print (" sarg = ",sarg), */ return (%fac*subst (sarg,_mu%,%rterm)))) else ( /* print (" case expression contains two or more KD's "), */ for %k thru length(KL) do (Kargs : KargsL[%k], /* print (" %k = ",%k," Kargs = ",Kargs), */ if inlist (Kargs,_mu%) then ( kargsL1 : [], for j thru length(KargsL) do if KargsL[j] # Kargs then kargsL1 : cons(KargsL[j],kargsL1), kargsL1 : reverse(kargsL1), /* print (" kargsL1 = ",kargsL1), */ kde : listToProd( map ('lambda([x], apply('KD, x)),kargsL1)), /* print (" kde = ",kde), */ if length (delete (_mu%, Kargs)) = 0 then (Kreturn : 4*%fac*kde*%rterm, return()) else (sarg : part (delete (_mu%, Kargs),1), /* print (" sarg = ",sarg), */ Kreturn : %fac*subst (sarg,_mu%,kde*%rterm), return()))), return(Kreturn))), acon1(%ee, _mu%))$ /************************************ end of acon 6-13-2017 **********************************/ /* mcon passes one term _b% at a time to mcon1 */ /* old buggy version of mcon1 mcon1(_b%,[xL]) := block ([_qsum%,jj ], _qsum% : acon (_b%,xL[1]), for jj:2 thru length (xL) do _qsum% : acon (_qsum%,xL[jj]), expand (_qsum%))$ */ /* new mcon1 6/6/2018 example of call: mcon1 (expr, mu,nu) called by new mcon */ mcon1(_b%,[xL]) := block ([_qsum%,jj,zeroans:false ], if _b% = 0 then return(0), _qsum% : acon (_b%,xL[1]), if _qsum% = 0 then return(0), for jj:2 thru length (xL) do (_qsum% : acon (_qsum%,xL[jj]), if _qsum% = 0 then (zeroans:true, return())), if zeroans then 0 else expand (_qsum%))$ /* old version of mcon accepts eps4*eps4 product and doesn't use ratexpand see far below for new version of mcon mcon (_e%,[vL]) := block ([vp,_ec%,mysumL,myrsum,_a% ], /* disp (" mcon, 4-9-11 "), */ if _e% = 0 then return(0), _ec% : expand (_e%), mysumL : sumToList (_ec%), myrsum : 0, for _a% in mysumL do myrsum : myrsum + apply ('mcon1,cons(_a%,vL)), myrsum )$ */ /* find_eps4_prod(expr,n1,n2,...) returns [mult_fac,eps_prod,found_indicesL] if an expression eps4*eps4L is found with both factors containing at least one and the same of the contraction indices n1,n2,... , with mult_fac being the rest of the expression, and found_indicesL being a sublist of the contraction indices provided which are contained in each factor of the eps4*eps4L product. otherwise, returns [expr,false,[]], or, if fatal error found, returns [false,false,[] ]. wL = list of one or more symbolic contraction indices */ find_eps4_prod (anexpr,[wL]) := block ([nfac,rfac,opsL,eps4_list, the_rest,eps4L_list, eps4_arg_list, eps4L_arg_list, eps4_found:false, eps4_index:[],eps4_pos,eargL, eps4L_found:false, eps4L_index:[],eps4L_pos,common_index:[], eps4_product], if debug then print ("find_eps4_prod"), if debug then print ("anexpr = ",anexpr), if debug then print ("wL = ",wL), if length(wL) = 0 then (print("find_eps4_prod: no symbolic contraction indices found "), return ([false,false,[]])), /* fatal error */ /* nfac is numerical factor which may include a sign */ [nfac,rfac] : NDfac(anexpr), if debug then print ("nfac = ",nfac), if debug then print ("rfac = ",rfac), /* does rfac contain at least one eps4 and at least one eps4L? */ opsL : opList (prodToList (rfac)), if debug then print ("opsL = ", opsL), if not inlist (opsL,eps4) or not inlist (opsL,eps4L) then return ([anexpr,false,[]]), /* case: we may find suitable product */ if debug then print ("look at details..."), /* eps4_list contains all the eps4's in the product */ [eps4_list, the_rest] : strip_ops2(rfac,eps4), if debug then print (" eps4_list = ",eps4_list), if debug then print (" the_rest = ", the_rest), /* eps4L_list contains all the eps4L's in the product */ [eps4L_list, the_rest] : strip_ops2(the_rest,eps4L), if debug then print (" eps4L_list = ",eps4L_list), if debug then print (" the_rest = ", the_rest), eps4_arg_list : map ('args,eps4_list), if debug then print (" eps4_arg_list = ",eps4_arg_list), eps4_arg_list : map ('lambda ([xL],sublist (xL,symbolp)), eps4_arg_list), if debug then print ("eps4_arg_list = ",eps4_arg_list), eps4L_arg_list : map ('args,eps4L_list), if debug then print (" eps4L_arg_list = ",eps4L_arg_list), eps4L_arg_list : map ('lambda ([xL],sublist (xL,symbolp)), eps4L_arg_list), if debug then print ("eps4L_arg_list = ",eps4L_arg_list), /* find position in eps4_list of the eps4 which contains any of the given contraction indices */ for j thru length (eps4_arg_list) do (eargL : eps4_arg_list[j], for k in wL do if inlist (eargL,k) then (eps4_found : true, eps4_pos : j, eps4_index : cons (k, eps4_index)), if eps4_found then (eps4_index : reverse (eps4_index), return () )), if not eps4_found then return ([anexpr,false,[]]), /* case eps4_found = true */ if debug then print ( "eps4_index = ",eps4_index), if debug then print ( "eps4_pos = ", eps4_pos), /* find position in eps4L_list of the eps4L which contains any of the given contraction indices */ for j thru length (eps4L_arg_list) do (eargL : eps4L_arg_list[j], for k in wL do if inlist (eargL,k) then (eps4L_found : true, eps4L_pos : j, eps4L_index : cons (k, eps4L_index)), if eps4L_found then (eps4L_index : reverse (eps4L_index), return () )), if not eps4L_found then return ([anexpr,false,[]]), /* case eps4_found = true and eps4L_found = true */ /* largest number of contraction indices which appear in each list? */ if debug then print ( "eps4L_index = ",eps4L_index), if debug then print ( "eps4L_pos = ", eps4L_pos), for k in eps4_index do if inlist (eps4L_index,k) then common_index : cons (k, common_index), if length (common_index) = 0 then return ([anexpr,false,[]]), common_index : reverse (common_index), eps4_product : eps4_list[eps4_pos] * eps4L_list[eps4L_pos], if debug then print (" eps4_product = ",eps4_product), for j thru length (eps4_list) do if j # eps4_pos then the_rest : the_rest*eps4_list[j], for j thru length (eps4L_list) do if j # eps4L_pos then the_rest : the_rest*eps4L_list[j], if debug then print (" the_rest = ", the_rest), the_rest : the_rest*nfac, if debug then print (" the_rest = ", the_rest), [the_rest, eps4_product, common_index] )$ /* mcon (expr,n1,n2,...) new mcon 5/30/2017 calls find_eps4_prod and coneps2 to contract over a product of two eps factors (which should be of the form eps4*eps4L ) as part of the expression passed to mcon. (%i32) prodToList(-f1*f2*eps4[n1,n2,n3,n4]*eps4L[n1,n2,n3,m4]); (%o32) [eps4L[m4,n1,n2,n3],eps4[n1,n2,n3,n4]] */ mcon (_e%,[vL]) := block ([ _ec%,mysumL,myrsum,_a%, pfac,eprod,ep_indexL, anerr:false, ep_con ], if debug then print (" mcon, 5-30-2017 "), if _e% = 0 then return(0), _ec% : expand (_e%), mysumL : sumToList (_ec%), if debug then print (" mysumL = ",mysumL), myrsum : 0, if debug then print(" vL = ",vL," myrsum = ",myrsum), for _a% in mysumL do ( if debug then print (" _a% = ",_a%), [ pfac,eprod,ep_indexL ] : apply ('find_eps4_prod, cons (_a%,vL)), if debug then print(" pfac = ",pfac), if debug then print (" eprod = ",eprod), if debug then print (" ep_indexL = ",ep_indexL), if pfac = false then ( print (" mcon: fatal error"), anerr : true, return ()), /* get out of do loop */ if eprod = false then (if debug then print (" case eprod = false, just use mcon1 "), myrsum : myrsum + apply ('mcon1,cons(_a%,vL))) else if length(ep_indexL) = length(vL) then ( if debug then print (" case length(ep_indexL) = length(vL) "), ep_con : apply ('coneps2, cons(eprod,ep_indexL)), if debug then print (" ep_con = ",ep_con), myrsum : myrsum + pfac*ep_con) else (if debug then print (" case length(ep_indexL) < length(vL) "), for j in ep_indexL do vL : delete (j,vL), if debug then print (" new vL = ",vL), myrsum : myrsum + apply ('mcon1,cons(pfac*ep_con,vL)))), if debug then print (" out of do loop "), if debug then print (" myrsum = ",myrsum), if debug then print (" anerr = ",anerr), if anerr then done else expand (myrsum) )$ /* (%i2) mcon (KD(n1,n2)*eps4[n1,n3,n4,n5],n1); (%o2) eps4[n2,n3,n4,n5] (%i4) mcon (-4*KD(n1,n2)*eps4[n1,n3,n4,n5]/3,n1); (%o4) -(4*eps4[n2,n3,n4,n5])/3 (%i3) mcon (KD(n1,n2)*KD(n6,n7)*eps4[n1,n3,n4,n5],n1); (%o3) eps4[n2,n3,n4,n5]*KD(n6,n7) (%i5) mcon(f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4],n1,n2); (%o5) 2*f1*f2*kron_delta(m3,n4)*kron_delta(m4,n3) -2*f1*f2*kron_delta(m3,n3)*kron_delta(m4,n4) (%i6) mcon(f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]*KD(n5,n6),n1,n2); (%o6) 2*f1*f2*kron_delta(m3,n4)*kron_delta(m4,n3)*KD(n5,n6) -2*f1*f2*kron_delta(m3,n3)*kron_delta(m4,n4)*KD(n5,n6) tools used and practice: strip_atoms(Aprod) := listToProd (prodToList (Aprod))$ (%i2) strip_atoms(Aprod) := listToProd (prodToList (Aprod))$ (%i3) strip_atoms (f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]); (%o3) eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i4) aprod : f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]; (%i7) strip_non_atoms(Aprod) := (strip_ops2(Aprod, eps4)[2], strip_ops2 (%%,eps4L)[2])$ (%i8) strip_non_atoms(aprod); (%o8) f1*f2 (%i9) aprod; (%o9) f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i10) strip_atoms(aprod); (%o10) eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i2) myterm : -f1*f2*eps4[n1,n2,n3,n4]*eps4L[n1,n2,m3,m4]; (%o2) -f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i3) [fac,en] : NDfac(myterm); (%o3) [-1,f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]] (%i6) NDfac(2*myterm/3); (%o6) [-2/3,f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4]] /* prodToList ignores atoms , prodToList2 includes atoms */ (%i5) en1 : prodToList(en); (%o5) [eps4L[m3,m4,n1,n2],eps4[n1,n2,n3,n4]] (%i8) oopL : opList(en1); (%o8) [eps4L,eps4] (%i11) ne1 : pos (oopL,eps4L); (%o11) 1 (%i12) ne2 : pos(oopL,eps4); (%o12) 2 (%i13) pos(oopL,cos); (%o13) 0 (%i16) remL1(oopL,1); (%o16) [eps4] (%i17) remL1(oopL,2); (%o17) [eps4L] (%i18) remL2(oopL,1,2); (%o18) [] using strip_ops2 (%i36) en; (%o36) f1*f2*eps4L[m3,m4,n1,n2]*eps4[n1,n2,n3,n4] (%i37) strip_ops2(en,eps4); (%o37) [[eps4[n1,n2,n3,n4]],f1*f2*eps4L[m3,m4,n1,n2]] (%i38) strip_ops2(f1*f2*eps4L[m3,m4,n1,n2],eps4L); (%o38) [[eps4L[m3,m4,n1,n2]],f1*f2] */ /* useful for getting zero from roots of roots warning: rootcrunch squares the argument!! so if arg is not acually zero, the answer returned will be the square, thus possibly losing signs. */ rootcrunch(_x%):= (expand (_x%), expand (%%*conj(%%)),rootscontract (%%), factor (%%))$ trigcrunch1(expr,x) := (trigsimp (expand (demoivre (expr))),to_half_angle(%%,x), rootcrunch (%%) )$ trigcrunch(expr,x) := (trigsimp (expand (demoivre (expr))),to_ao2(%%,x), rootcrunch (%%) )$