/************************************************************************ dirac3.mac is a package of Maxima functions which contains code for high energy physics calculations and is part of the Dirac package. Maxima by Example Ch 12, Dirac Algebra and Quantum Electrodynamics Copyright (C) 2010, 2011,2017,2018 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/. ************************************************************************/ /* dirac3.mac 5/5/2017 */ /* Nlast initialization. declare D and Gm to be symmetric. use Mindex to get a default set of index symbols. definitions for: conj divout div_tb dummyp eps4 eps4L fr_ao2 gmet pullfac pfactor str_to_int take_parts to_half_angle to_ao2 ts */ print("dirac3.mac"); kill(all); reset(); /* (%i2) apply('f,[a,b]); (%o2) f(a,b) (%i3) a: %pi; (%o3) %pi (%i4) b : %pi/2; (%o4) %pi/2 (%i5) [a,b]; (%o5) [%pi,%pi/2] (%i6) kill(a,b); (%o6) done (%i7) [a,b]; (%o7) [a,b] */ if not lfreeof (rules, eps4rule1) then apply ('kill, [eps4,eps4L] )$ /* remove all global arrays currently existing */ if length (arrays) > 0 then apply ('remarray, arrays)$ /* load("new-timedate.lisp")$ mydate : apply ('sconcat, rest (charlist (timedate(+7)),-15))$ _binfo% : "Maxima 5.24.0"$ */ load ("simplifying-new.lisp")$ print ("simplifying-new.lisp")$ load(dgtrace3)$ print ("dgtrace3.mac")$ load (dgcon3)$ print ("dgcon3.mac")$ load(dgeval3)$ print ("dgeval3.mac")$ load (dgmatrix3)$ print ( "dgmatrix3.mac ")$ /* load(mandelstam)$ print("mandelstam.mac")$ */ /* load(qcd1); disp ("qcd1.mac"); */ /* test_expr : [a,b,c]$ display ( test_expr)$ */ /******************************************/ /* str_to_int (string) 6-9-10, see examples below a reminder of some string manipulations and called by dummyp */ str_to_int(sval) := block ([cs,rsum:0,pp,jj ], cs : reverse (charlist (sval)), for jj thru length (cs) do ( pp : cint (cs[jj]) - 48, rsum : rsum + pp*10^(jj -1)), rsum )$ /* (%i1) display2d:false$ (%i2) sn : string(Q572); (%o2) "Q572" (%i3) cs : charlist(sn); (%o3) ["Q","5","7","2"] (%i4) sn1 : sremove ("Q",sn); (%o4) "572" (%i14) str_to_int ("572"); (%o14) 572 (%i15) str_to_int ("52"); (%o15) 52 (%i2) display2d:false$ (%i3) sremove("N","N20"); (%o3) "20" (%i4) sremove("N","N20N"); (%o4) "20" (%i5) sremovefirst ("N","N20N"); (%o5) "20N" */ /*** dummyp 3-25-11 *******************/ /* dummy summation variables generated by traces of gamma5 times other dirac matrices are automatically of the form [N1,N2,...] depending on how many are needed. We want dummyp (N2345) --> true for example. dummyp is used by sum_eps1 in dgeval2.mac. */ dummyp(expr) := block ([se,dval:false], se : string (expr), if first (charlist (se)) = "N" then if integerp (str_to_int (sremovefirst ("N",se))) then dval : true, dval)$ /* (%i25) dummyp (N1); (%o25) true (%i26) dummyp (N21); (%o26) true (%i27) dummyp (n1); (%o27) false (%i28) dummyp (n20); (%o28) false */ /* complex conjugation for expressions containg %i's and real numbers */ conj(_e%) := subst (-%i,%i,_e%)$ /* replace square of momentum symbol by symbols for energy and mass 5-5-2017 */ psq_to_Esq (expr,%p,%E,%m) := expand (ratsubst (%E^2 - %m^2, %p^2, expr))$ /* p_Em is old name for psq_to_Esq, is already used in some of our old batch files for a more specific replacement. p_Em (expr,%p,%E,%m) := expand (ratsubst (%E^2 - %m^2, %p^2, expr))$ */ Esq_to_psq (expr,%p,%E,%m) := expand (ratsubst (%m^2 + %p^2 , %E^2, expr))$ /* Ep_to_m (expr,p,E,m) replaces sqrt(E - p) * sqrt(E + p) by m */ Ep_to_m(eexpr,%p,%E,%m) := expand( ratsubst (%m, sqrt(%E - %p)*sqrt(%E + %p), eexpr) )$ /* Ep_to_mM(e,p,E,m,M) applies to m1=m2=m3=m4=m 2 to 2 reaction in center of momentum frame, effective system mass M = 2*E, replaces E by M/2, p^2 by M^2/4 - m^2 */ Ep_to_mM(eexpr,%p,%E,%m,%M) := ( expand (ratsubst (%M^2/4 - %m^2, %p^2,eexpr)), expand ( ratsubst (%M/2, %E, %%)))$ /***** take_parts take_parts(expr,n1,n2) returns the sum of the parts of expr from part(expr,n1) thru part(expr,n2). *********/ take_parts (expr,nni,nnf) := block ([rsum:0,jj], for jj:nni thru nnf do rsum : rsum + part (expr,jj), rsum )$ /***************/ /* div_tb = divide numerator (top) and denominator (bottom) by same symbol */ /* transferred from tempalg78.mac 8-6-10 */ div_tb(expr,_x%) := block ([ top,bottom], top : num(expr), bottom : denom(expr), top : expand (top/_x%), /* display (top), */ bottom : expand (bottom/_x%), /* display (bottom), */ top/bottom)$ /* if trigsimp is not working, then use */ ts(e,v) := ratsubst(1,cos(v)^2 + sin(v)^2, e)$ /* convert expression containing an angle a into an equivalent expression containing instead the half angle a/2 */ to_half_angle(expr,a) := (ratsubst (2*sin(a/2)*cos(a/2),sin(a),expr), ratsubst (cos(a/2)^2 - sin(a/2)^2,cos(a),%%), trigsimp (ratsimp (expand (%%))))$ /* to_ao2(expr, angle) convert an expression which is the sum of expressions to half angle form. "to a over 2" calls to_half_angle */ to_ao2(expr,aa) := block ([pp,jj,rsum], if atom(expr) then return(expr), if op(expr) = "+" then ( rsum : 0, for jj thru length (expr) do ( pp : part (expr,jj), rsum : rsum + expand (to_half_angle (pp,aa)) ) ) else rsum : expand ( to_half_angle(expr,aa)), rsum )$ /* convert from a/2 trig args to a. almost the same as ev (expr, trigexpand, halfangles) */ fr_ao2(expr,a) := ( ratsubst (sin(a)/2,cos(a/2)*sin(a/2),expr), ratsubst (sqrt ( (1+cos(a))/2 ),cos (a/2),%% ), ratsubst (sqrt ( (1-cos(a))/2 ),sin (a/2),%% ), rootscontract (%%), trigsimp(%%), expand (%%) )$ divout(eexpr,_aa) := map (lambda([_x], ratsimp(_x/_aa)),eexpr)$ pullfac1(eexpr,_n%) := if op(eexpr) = "+" then _n%*divout (eexpr,_n%) else eexpr$ /* suggested version by Stavros Macrakis */ pullfac(eexpr,_n%) := if op(eexpr) = "+" then _n%*map(lambda([zz],zz/_n%), eexpr) else eexpr$ /* partial factorization pulls out factor var^pp */ pfactor(expr,var,pp) := block ( [exprc,cexpr,subexpr,restexpr ], exprc : copy(expr), cexpr : coeff(exprc,var,pp), /* display (cexpr), */ subexpr : var^pp*cexpr, /* display (subexpr), */ restexpr : exprc - expand(subexpr), subexpr + restexpr )$ /* (%i95) t2a; (%o95) cos(th)*E^2+E^2+m^2*cos(th)-m^2 (%i109) pfactor (t2a,E,2); (%o109) (cos(th)+1)*E^2+m^2*cos(th)-m^2 (%i110) pfactor (%,m,2); (%o110) (cos(th)+1)*E^2+m^2*(cos(th)-1) (%i111) (pfactor (t2a,E,2), pfactor (%%,m,2)); (%o111) (cos(th)+1)*E^2+m^2*(cos(th)-1) */ /***************************************/ uv_simp1(%expr,%p,%E,%m) := (trigsimp (%expr),rootscontract(%%), psq_to_Esq (%%, %p,%E,%m))$ uv_simp (%expr,%p,%E,%m) := trigsimp (uv_simp1(%expr,%p,%E,%m))$ /* uv_simp (%expr,%p,%E,%m) := (trigsimp (%expr),rootscontract(%%), psq_to_Esq (%%, %p,%E,%m), trigsimp(%%))$ */ Nlast : 0$ declare ( [D,Gm], symmetric )$ /* old method, before definition of Mindex: declare ([n1,n2,n3,n4,n5,n6,n7,n8,n9,n10], index); */ /* Mindex(..) is defined in dgtrace2.mac and used here */ /* greek symbols, for example, Mindex (la,mu,nu,rh,si,ta,al,be,ga,de,ep)$ */ /* or numbered indices Mindex (n1,n2,n3,n4,n5,n6,n7,n8,n9,n10)$ */ /* or both */ Mindex (n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,la,mu,nu,rh,si,ta,al,be,ga,de,ep)$ Mmass (m,M)$ array (gmet, fixnum, 3,3)$ fillarray (gmet, [1,0,0,0, 0,-1,0,0, 0,0,-1,0, 0,0,0,-1])$ declare (eps4,antisymmetric)$ tellsimpafter ( eps4[0,1,2,3], 1 )$ declare (eps4L,antisymmetric)$ tellsimpafter ( eps4L[0,1,2,3], -1 )$ /* declare (epsU,antisymmetric)$ tellsimpafter ( epsU[0,1,2,3], 1 )$ declare (epsL,antisymmetric)$ tellsimpafter ( epsL[0,1,2,3], -1 )$ */ /* interactive check on correct contraction values: (%i13) declare (epsU,antisymmetric)$ (%i14) tellsimpafter ( epsU[0,1,2,3], 1 )$ (%i15) declare (epsL,antisymmetric)$ (%i16) tellsimpafter ( epsL[0,1,2,3], -1 )$ (%i17) epsU[0,1,2,3]; (%o17) 1 (%i18) epsL[0,1,2,3]; (%o18) -1 (%i19) sum( sum( sum(epsU[n1,n2,n3,0]*epsL[n1,n2,n3,0],n1,0,3),n2,0,3),n3,0,3); (%o19) -6 (%i20) sum( sum( sum( sum(epsU[n1,n2,n3,n4]*epsL[n1,n2,n3,n4],n1,0,3),n2,0,3),n3,0,3),n4,0,3); (%o20) -24 */ /**********************************************/ /* note: the function Mscalar (in dgtrace3.mac) displays scalarL */ /* display2d:false$ */ print ("current symbol assignments")$ /*printf (true, " ~& ~a ", " current symbol assignments...")$ */ Mscalar (c1,c2,c3,c4,c5,c6,c7,c8,c9,c10)$ /* display (indexL,massL,Nlast )$ */ printf (true, " ~& ~a ~a "," indexL = ",indexL)$ printf (true, " ~& ~a ~a ","massL = ",massL)$ printf (true, " ~& ~a ~a "," Nlast = ",Nlast)$ print (" reserved program capital letter name use: ")$ print (" Chi, Con, D, Eps, EpsL, G, G(1), G5, G5p ")$ print (" Gam, Gm, Gtr, KD, LI, Nlast, UI, P, S ")$ print (" Sig, Sigb, SIG, UU, VP, VV ")$ print (" I2, Z2, CZ2, I4, Z4, CZ4, RZ4, N1,N2,... ")$ print (" reserved array names: gmet, eps4, eps4L ")$ printf (true, " ~& ~a ~a "," invar_flag = ",invar_flag)$ printf (true, " ~& ~a ~a "," stu_flag = ",stu_flag)$ /* collect_terms (expr, fac) makes use of Maxima's collectterms */ collect_terms (_e%, _f%) := block ( [_a%], ratsubst (_a%, _f%, _e%), collectterms (%%, _a%), ev (%%, _a% = _f%))$ /* (%i33) expr : 1/(-(k*cos(th))/m+k/m+1); (expr) 1/(-(k*cos(th))/m+k/m+1) (%i34) expr1 : denom (expr); (expr1) -(k*cos(th))/m+k/m+1 (%i35) collect_terms (expr1, k/m); (%o35) (k*(1-cos(th)))/m+1 */ display2d: true$ trigexpand : true$ halfangles : false$