Applications, Examples and Libraries

Share your work here

Playing mini-golf recently, I realized that my protractor can only help me so far since it can't calculate the speed of the swing needed.  I decided a more sophisticated tool was needed and modeled a trick-shot in MapleSim.

To start, I laid out the obstacles, the ball and club, the ground, and some additional visualizations in the MapleSim environment.

 

When running the simulation, my first result wasn't even close to the hole (similar to when I play in real life!).

 

The model clearly needed to be optimized. I went to the Optimization app in MapleSim (this can be found under Add Apps or Templates  on the left hand side).

 

Inside the app I clicked "Load System" then selected the parameters I wanted to optimize.

 

For this case, I'm optimizing 's' (the speed of the club) and 'theta' (the angle of the club). For the Objective Function I added a Relative Translation Sensor to the model and attached a probe to the Vector Norm of the output.

 

Inside the app, I switched to the Objective Function section.  Selecting Probes, I added the new probe as the Objective Function by giving it a weight of 1.

 

 

Scrolling down to "Execute Parameter Optimization", I checked the "Use Global Optimization Toolbox" checkbox, and clicked Run Parameter Optimization.

 

Following a run time of 120 seconds, the app returns the graph of the objective function. 

 

Below the plot, optimal values for the parameters are given. Plugging these back into the parameter block for the simulation we see that the ball does in fact go into the hole. Success!

 

 

Mini_golf_Global_Optimization.msim

Can you guess what P() produces, without executing it?

P:=proc(N:=infinity) local q,r,t,k,n,l,h, f;
q,r,t,k,n,l,h := 1,0,1,1,3,3,0:
while h<N do 
   if 4*q+r-t < n*t
   then f:=`if`(++h mod 50=0,"\n",`if`(h mod 10=0," ","")); printf("%d"||f,n);   
        q,r,t,k,n,l := 10*q,10*(r-n*t),t,k,iquo(10*(3*q+r),t)-10*n,l
   else q,r,t,k,n,l := q*k,(2*q+r)*l,t*l,k+1,iquo(q*(7*k+2)+r*l,t*l),l+2
   fi
od: NULL
end:

I hope you will like it (maybe after execution).

 

Feynman Diagrams
The scattering matrix in coordinates and momentum representation

  

Mathematical methods for particle physics was one of the weak spots in the Physics package. There existed a FeynmanDiagrams command, but its capabilities were too minimal. People working in the area asked for more functionality. These diagrams are the cornerstone of calculations in particle physics (collisions involving from the electron to the Higgs boson), for example at the CERN. As an introduction for people curious, not working in the area, see "Why Feynman Diagrams are so important".

  

This post is thus about a new development in Physics: a full rewriting of the FeynmanDiagrams command, now including a myriad of new capabilities (mainly a. b. and c. in the Introduction), reversing the previous status of things entirely. This is work in collaboration with Davide Polvara from Durham University, Centre for Particle Theory.

  


The complexity of this material is high, so the introduction to the presentation below is as brief as it can get, emphasizing the examples instead. This material is reproducible in Maple 2019.2 after installing the Physics Updates, v.598 or higher.

  

 

  

At the end they are attached the worksheet corresponding to this presentation and a PDF version of it, as well as the new FeynmanDiagrams help page with all the explanatory details.

Introduction

  

A scattering matrix S relates the initial and final states, `#mfenced(mrow(mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")` and `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")`, of an interacting system. In an 4-dimensional spacetime with coordinates X, S can be written as:

S = T(exp(i*`#mrow(mo("&int;"),mi("L"),mo("&ApplyFunction;"),mfenced(mi("X")),mo("&DifferentialD;"),msup(mi("X"),mn("4")))`))

  

where i is the imaginary unit  and L is the interaction Lagrangian, written in terms of quantum fields  depending on the spacetime coordinates  X. The T symbol means time-ordered. For the terminology used in this page, see for instance chapter IV, "The Scattering Matrix", of ref.[1] Bogoliubov, N.N., and Shirkov, D.V. Quantum Fields.

  

This exponential can be expanded as

S = 1+S[1]+S[2]+S[3]+`...`

  

where

S[n] = `#mrow(mo("&ApplyFunction;"),mfrac(msup(mi("i"),mi("n")),mrow(mi("n"),mo("&excl;")),linethickness = "1"),mo("&InvisibleTimes;"),mo("&int;"),mi("&hellip;"),mo("&InvisibleTimes;"),mo("&int;"),mi("T"),mo("&ApplyFunction;"),mfenced(mrow(mi("L"),mo("&ApplyFunction;"),mfenced(mi("\`X__1\`")),mo("&comma;"),mi("&hellip;"),mo("&comma;"),mi("L"),mo("&ApplyFunction;"),mfenced(mi("\`X__n\`")))),mo("&InvisibleTimes;"),mo("&DifferentialD;"),msup(mi("\`X__1\`"),mn("4")),mo("&InvisibleTimes;"),mi("&hellip;"),mo("&InvisibleTimes;"),mo("&DifferentialD;"),msup(mi("\`X__n\`"),mn("4")))`

  

and T(L(X[1]), `...`, L(X[n])) is the time-ordered product of n interaction Lagrangians evaluated at different points. The S matrix formulation is at the core of perturbative approaches in relativistic Quantum Field Theory.

  

In connection, the FeynmanDiagrams  command has been rewritten entirely for Maple 2020. In brief, the new functionality includes computing:

a. 

The expansion S = 1+S[1]+S[2]+S[3]+`...` in coordinates representation up to arbitrary order (the limitation is now only your hardware)

b. 

The S-matrix element `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` in momentum representation up to arbitrary order for given number of loops and initial and final particles (the contents of the `#mfenced(mrow(mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")` and `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")` states); optionally, also the transition probability density, constructed using the square of the scattering matrix element abs(`#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")`)^2, as shown in formula (13) of sec. 21.1 of ref.[1].

c. 

The Feynman diagrams (drawings) related to the different terms of the expansion of S or of its matrix elements `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")`.

  

Interaction Lagrangians involving derivatives of fields, typically appearing in non-Abelian gauge theories, are also handled, and several options are provided enabling restricting the outcome in different ways, regarding the incoming and outgoing particles, the number of loops, vertices or external legs, the propagators and normal products, or whether to compute tadpoles and 1-particle reducible terms.

 

Examples

 

For illustration purposes set three coordinate systems , and set phi to represent a quantum operator

with(Physics)

Setup(mathematicalnotation = true, coordinates = [X, Y, Z], quantumoperators = phi)

`Systems of spacetime coordinates are:`*{X = (x1, x2, x3, x4), Y = (y1, y2, y3, y4), Z = (z1, z2, z3, z4)}

 

_______________________________________________________

 

[coordinatesystems = {X, Y, Z}, mathematicalnotation = true, quantumoperators = {phi}]

(1.1)

Let L be the interaction Lagrangian

L := lambda*phi(X)^4

lambda*Physics:-`^`(phi(X), 4)

(1.2)

The expansion of S in coordinates representation, computed by default up to order = 3 (you can change that using the option order = n), by definition containing all possible configurations of external legs, displaying the related Feynman Diagrams, is given by

%eval(S, `=`(order, 3)) = FeynmanDiagrams(L, diagrams)

 

 

 

%eval(S, order = 3) = 1+%FeynmanIntegral(lambda*_GF(_NP(phi(X), phi(X), phi(X), phi(X))), [[X]])+%FeynmanIntegral(16*lambda^2*_GF(_NP(phi(X), phi(X), phi(X), phi(Y), phi(Y), phi(Y)), [[phi(X), phi(Y)]])+96*lambda^2*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)], [phi(X), phi(Y)]])+72*lambda^2*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)]]), [[X], [Y]])+%FeynmanIntegral(1728*lambda^3*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y), phi(Z), phi(Z)), [[phi(X), phi(Z)], [phi(X), phi(Y)], [phi(Z), phi(Y)]])+2592*lambda^3*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y)), [[phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Z), phi(Y)], [phi(Z), phi(Y)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y), phi(Z), phi(Z)), [[phi(X), phi(Y)], [phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+3456*lambda^3*_GF(_NP(phi(X), phi(X)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+576*lambda^3*_GF(_NP(phi(X), phi(X), phi(X), phi(Y), phi(Y), phi(Z), phi(Z), phi(Z)), [[phi(X), phi(Y)], [phi(Y), phi(Z)]]), [[X], [Y], [Z]])

(1.3)


The expansion of S  in coordinates representation to a specific order shows in a compact way the topology of the underlying Feynman diagrams. Each integral is represented with a new command, FeynmanIntegral , that works both in coordinates and momentum representation. To each term of the integrands corresponds a diagram, and the correspondence is always clear from the symmetry factors.

In a typical situation, one wants to compute a specific term, or scattering process, instead of the S matrix up to some order with all possible configurations of external legs. For example, to compute only the terms of this result that correspond to diagrams with 1 loop use numberofloops = 1 (for tree-level, use numerofloops = 0)

%eval(S, [`=`(order, 3), `=`(loops, 1)]) = FeynmanDiagrams(L, numberofloops = 1, diagrams)

%eval(S, [order = 3, loops = 1]) = %FeynmanIntegral(72*lambda^2*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)]]), [[X], [Y]])+%FeynmanIntegral(1728*lambda^3*_GF(_NP(phi(X), phi(X), phi(Y), phi(Y), phi(Z), phi(Z)), [[phi(X), phi(Z)], [phi(X), phi(Y)], [phi(Z), phi(Y)]]), [[X], [Y], [Z]])

(1.4)


In the result above there are two terms, with 4 and 6 external legs respectively.

A scattering process with matrix element `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` in momentum representation, corresponding to the term with 4 external legs (symmetry factor = 72), could be any process where the total number of incoming + outgoing parties is equal to 4. For example, one with 2 incoming and 2 outgoing particles. The transition probability for that process is given by

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;",mathcolor = "olive")),open = "&lang;",close = "&rang;")` = FeynmanDiagrams(L, incomingparticles = [phi, phi], outgoingparticles = [phi, phi], numberofloops = 1, diagrams)

 

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&comma;",mathcolor = "olive"),mi("&phi;",fontstyle = "normal",mathcolor = "olive"),mo("&InvisibleTimes;",mathcolor = "olive")),open = "&lang;",close = "&rang;")` = %FeynmanIntegral((9/8)*lambda^2*Dirac(-P__3-P__4+P__1+P__2)/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1-P__2-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2*Dirac(-P__3-P__4+P__1+P__2)/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__3-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2*Dirac(-P__3-P__4+P__1+P__2)/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__4-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.5)

When computing in momentum representation, only the topology of the corresponding Feynman diagrams is shown (i.e. the diagrams associated to the corresponding Feynman integral in coordinates representation).

The transition matrix element `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` is related to the transition probability density dw (formula (13) of sec. 21.1 of ref.[1]) by

dw = (2*Pi)^(3*s-4)*n__1*`...`*n__s*abs(F(p[i], p[f]))^2*delta(sum(p[i], i = 1 .. s)-(sum(p[f], f = 1 .. r)))*` d `^3*p[1]*` ...`*`d `^3*p[r]

where n__1*`...`*n__s represent the particle densities of each of the s particles in the initial state `#mfenced(mrow(mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&verbar;",close = "&rang;")`, the delta (Dirac) is the expected singular factor due to the conservation of the energy-momentum and the amplitude F(p[i], p[f])is related to `#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` via

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = F(p[i], p[f])*delta(sum(p[i], i = 1 .. s)-(sum(p[f], f = 1 .. r)))

To directly get the probability density dw instead of`#mfenced(mrow(mo("&InvisibleTimes;"),mi("f"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("i"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")`use the option output = probabilitydensity

FeynmanDiagrams(L, incomingparticles = [phi, phi], outgoingparticles = [phi, phi], numberofloops = 1, output = probabilitydensity)

Physics:-FeynmanDiagrams:-ProbabilityDensity(4*Pi^2*%mul(n[i], i = 1 .. 2)*abs(F)^2*Dirac(-P__3-P__4+P__1+P__2)*%mul(dP_[f]^3, f = 1 .. 2), F = %FeynmanIntegral((9/8)*lambda^2/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1-P__2-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__3-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])+%FeynmanIntegral((9/8)*lambda^2/(Pi^6*(E__1*E__2*E__3*E__4)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__4-p__2)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]]))

(1.6)

In practice, the most common computations involve processes with 2 or 4 external legs. To restrict the expansion of the scattering matrix in coordinates representation to that kind of processes use the numberofexternallegs option. For example, the following computes the expansion of S up to order = 3, restricting the outcome to the terms corresponding to diagrams with only 2 external legs

%eval(S, [`=`(order, 3), `=`(legs, 2)]) = FeynmanDiagrams(L, numberofexternallegs = 2, diagrams)

%eval(S, [order = 3, legs = 2]) = %FeynmanIntegral(96*lambda^2*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Y)], [phi(X), phi(Y)]]), [[X], [Y]])+%FeynmanIntegral(3456*lambda^3*_GF(_NP(phi(X), phi(X)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]]), [[X], [Y], [Z]])

(1.7)


This result shows two Feynman integrals, with respectively 2 and 3 loops, the second integral with two terms. The transition probability density in momentum representation for a process related to the first integral (1 term with symmetry factor = 96) is then

FeynmanDiagrams(L, incomingparticles = [phi], outgoingparticles = [phi], numberofloops = 2, diagrams, output = probabilitydensity)

Physics:-FeynmanDiagrams:-ProbabilityDensity((1/2)*%mul(n[i], i = 1 .. 1)*abs(F)^2*Dirac(-P__2+P__1)*%mul(dP_[f]^3, f = 1 .. 1)/Pi, F = %FeynmanIntegral(%FeynmanIntegral(((3/8)*I)*lambda^2/(Pi^7*(E__1*E__2)^(1/2)*(p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__3^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1-p__2-p__3)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]]), [[p__3]]))

(1.8)

In the above, for readability, the contracted spacetime indices in the square of momenta entering the amplitude F (as denominators of propagators) are implicit. To make those indices explicit, use the option putindicesinsquareofmomentum

F = FeynmanDiagrams(L, incoming = [phi], outgoing = [phi], numberofloops = 2, indices)

`* Partial match of  '`*indices*`' against keyword '`*putindicesinsquareofmomentum*`' `

 

F = %FeynmanIntegral(%FeynmanIntegral(((3/8)*I)*lambda^2*Dirac(-P__2[`~kappa`]+P__1[`~kappa`])/(Pi^7*(E__1*E__2)^(1/2)*(p__2[mu]*p__2[`~mu`]-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__3[nu]*p__3[`~nu`]-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1[beta]-p__2[beta]-p__3[beta])*(-P__1[`~beta`]-p__2[`~beta`]-p__3[`~beta`])-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]]), [[p__3]])

(1.9)

This computation can also be performed to higher orders. For example, with 3 loops, in coordinates and momentum representations, corresponding to the other two terms and diagrams in (1.7)

%eval(S[3], [`=`(legs, 2), `=`(loops, 3)]) = FeynmanDiagrams(L, legs = 2, loops = 3)

`* Partial match of  '`*legs*`' against keyword '`*numberoflegs*`' `

 

`* Partial match of  '`*loops*`' against keyword '`*numberofloops*`' `

 

%eval(S[3], [legs = 2, loops = 3]) = %FeynmanIntegral(3456*lambda^3*_GF(_NP(phi(X), phi(X)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]])+10368*lambda^3*_GF(_NP(phi(X), phi(Y)), [[phi(X), phi(Y)], [phi(X), phi(Z)], [phi(X), phi(Z)], [phi(Y), phi(Z)], [phi(Y), phi(Z)]]), [[X], [Y], [Z]])

(1.10)

A corresponding S-matrix element in momentum representation:

%eval(%Bracket(phi, S[3], phi), `=`(loops, 3)) = FeynmanDiagrams(L, incomingparticles = [phi], outgoingparticles = [phi], numberofloops = 3)

%eval(%Bracket(phi, S[3], phi), loops = 3) = %FeynmanIntegral(%FeynmanIntegral(%FeynmanIntegral((9/32)*lambda^3*Dirac(-P__2+P__1)/(Pi^11*(E__1*E__2)^(1/2)*(p__3^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__4^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__5^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-p__3-p__4-p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+P__2+p__3+p__4+p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__3]]), [[p__4]]), [[p__5]])+2*%FeynmanIntegral(%FeynmanIntegral(%FeynmanIntegral((9/32)*lambda^3*Dirac(-P__2+P__1)/(Pi^11*(E__1*E__2)^(1/2)*(p__3^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__4^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__5^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-p__3-p__4-p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+p__4+p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__3]]), [[p__4]]), [[p__5]])+%FeynmanIntegral(%FeynmanIntegral((1/2048)*lambda*Dirac(-P__2+P__1)*%FeynmanIntegral(576*lambda^2/((p__2^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-p__2-p__4-p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])/(Pi^11*(E__1*E__2)^(1/2)*(p__4^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*(p__5^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((-P__1+p__4+p__5)^2-m__phi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__4]]), [[p__5]])

(1.11)

Consider the interaction Lagrangian of Quantum Electrodynamics (QED). To formulate this problem on the worksheet, start defining the vector field A[mu].

Define(A[mu])

`Defined objects with tensor properties`

 

{A[mu], Physics:-Dgamma[mu], P__1[mu], P__2[mu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], p__1[mu], p__2[mu], p__3[mu], p__4[mu], p__5[mu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X), Physics:-SpaceTimeVector[mu](Y), Physics:-SpaceTimeVector[mu](Z)}

(1.12)

Set lowercase Latin letters from i to s to represent spinor indices (you can change this setting according to your preference, see Setup ), also the (anticommutative) spinor field will be represented by psi, so set psi as an anticommutativeprefix, and set A and psi as quantum operators

Setup(spinorindices = lowercaselatin_is, anticommutativeprefix = psi, op = {A, psi})

`* Partial match of  '`*op*`' against keyword '`*quantumoperators*`' `

 

_______________________________________________________

 

[anticommutativeprefix = {psi}, quantumoperators = {A, phi, psi}, spinorindices = lowercaselatin_is]

(1.13)

The matrix indices of the Dirac matrices  are written explicitly and use conjugate  to represent the Dirac conjugate conjugate(psi[j])

L__QED := alpha*conjugate(psi[j](X))*Dgamma[mu][j, k]*psi[k](X)*A[mu](X)

alpha*Physics:-`*`(conjugate(psi[j](X)), psi[k](X), A[mu](X))*Physics:-Dgamma[`~mu`][j, k]

(1.14)

Compute S[2], only the terms with 4 external legs, and display the diagrams: all the corresponding graphs have no loops

%eval(S[2], `=`(legs, 4)) = FeynmanDiagrams(L__QED, numberofvertices = 2, numberoflegs = 4, diagrams)

%eval(S[2], legs = 4) = %FeynmanIntegral(-2*alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(psi[k](X), A[mu](X), conjugate(psi[i](Y)), A[alpha](Y)), [[psi[l](Y), conjugate(psi[j](X))]])+alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(conjugate(psi[j](X)), psi[k](X), conjugate(psi[i](Y)), psi[l](Y)), [[A[mu](X), A[alpha](Y)]]), [[X], [Y]])

(1.15)

The same computation but with only 2 external legs results in the diagrams with 1 loop that correspond to the self-energy of the electron and the photon (page 218 of ref.[1])

%eval(S[2], `=`(legs, 2)) = FeynmanDiagrams(L__QED, numberofvertices = 2, numberoflegs = 2, diagrams)

 

 

%eval(S[2], legs = 2) = %FeynmanIntegral(-2*alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(psi[k](X), conjugate(psi[i](Y))), [[A[mu](X), A[alpha](Y)], [psi[l](Y), conjugate(psi[j](X))]])-alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(A[mu](X), A[alpha](Y)), [[psi[l](Y), conjugate(psi[j](X))], [psi[k](X), conjugate(psi[i](Y))]]), [[X], [Y]])

(1.16)

where the diagram with two spinor legs is the electron self-energy. To restrict the output furthermore, for example getting only the self-energy of the photon, you can specify the normal products you want:

%eval(S[2], [`=`(legs, 2), `=`(products, _NP(A, A))]) = FeynmanDiagrams(L__QED, numberofvertices = 2, numberoflegs = 2, normalproduct = _NP(A, A))

`* Partial match of  '`*normalproduct*`' against keyword '`*normalproducts*`' `

 

%eval(S[2], [legs = 2, products = _NP(A, A)]) = %FeynmanIntegral(alpha^2*Physics:-Dgamma[`~mu`][j, k]*Physics:-Dgamma[`~alpha`][i, l]*_GF(_NP(A[mu](X), A[alpha](Y)), [[conjugate(psi[j](X)), psi[l](Y)], [psi[k](X), conjugate(psi[i](Y))]]), [[X], [Y]])

(1.17)

The corresponding S-matrix elements in momentum representation

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = FeynmanDiagrams(L__QED, incomingparticles = [psi], outgoing = [psi], numberofloops = 1, diagrams)

 

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("&psi;",fontstyle = "normal"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = -%FeynmanIntegral((1/8)*Physics:-FeynmanDiagrams:-Uspinor[psi][i](P__1_)*conjugate(Physics:-FeynmanDiagrams:-Uspinor[psi][l](P__2_))*(-Physics:-g_[alpha, nu]+p__2[nu]*p__2[alpha]/m__A^2)*alpha^2*Physics:-Dgamma[`~alpha`][l, m]*Physics:-Dgamma[`~nu`][n, i]*((P__1[beta]+p__2[beta])*Physics:-Dgamma[`~beta`][m, n]+m__psi*Physics:-KroneckerDelta[m, n])*Dirac(-P__2+P__1)/(Pi^3*(p__2^2-m__A^2+I*Physics:-FeynmanDiagrams:-epsilon)*((P__1+p__2)^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.18)


In this result we see u[psi] spinor (see ref.[2]), and the propagator of the field A[mu] with a mass m[A]. To indicate that this field is massless use

Setup(massless = A)

`* Partial match of  '`*massless*`' against keyword '`*masslessfields*`' `

 

_______________________________________________________

 

[masslessfields = {A}]

(1.19)

Now the propagator for A[mu] is the one of a massless vector field:

FeynmanDiagrams(L__QED, incoming = [psi], outgoing = [psi], numberofloops = 1)

-%FeynmanIntegral(-(1/8)*Physics:-FeynmanDiagrams:-Uspinor[psi][i](P__1_)*conjugate(Physics:-FeynmanDiagrams:-Uspinor[psi][l](P__2_))*Physics:-g_[alpha, nu]*alpha^2*Physics:-Dgamma[`~alpha`][l, m]*Physics:-Dgamma[`~nu`][n, i]*((P__1[beta]+p__2[beta])*Physics:-Dgamma[`~beta`][m, n]+m__psi*Physics:-KroneckerDelta[m, n])*Dirac(-P__2+P__1)/(Pi^3*(p__2^2+I*Physics:-FeynmanDiagrams:-epsilon)*((P__1+p__2)^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.20)

The self-energy of the photon:

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = FeynmanDiagrams(L__QED, incomingparticles = [A], outgoing = [A], numberofloops = 1)

`#mfenced(mrow(mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("S"),mo("&InvisibleTimes;"),mo("&verbar;"),mo("&InvisibleTimes;"),mi("A"),mo("&InvisibleTimes;")),open = "&lang;",close = "&rang;")` = -%FeynmanIntegral((1/16)*Physics:-FeynmanDiagrams:-PolarizationVector[A][nu](P__1_)*conjugate(Physics:-FeynmanDiagrams:-PolarizationVector[A][alpha](P__2_))*(m__psi*Physics:-KroneckerDelta[l, n]+p__2[beta]*Physics:-Dgamma[`~beta`][l, n])*alpha^2*Physics:-Dgamma[`~alpha`][n, i]*Physics:-Dgamma[`~nu`][m, l]*((P__1[tau]+p__2[tau])*Physics:-Dgamma[`~tau`][i, m]+m__psi*Physics:-KroneckerDelta[i, m])*Dirac(-P__2+P__1)/(Pi^3*(E__1*E__2)^(1/2)*(p__2^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)*((P__1+p__2)^2-m__psi^2+I*Physics:-FeynmanDiagrams:-epsilon)), [[p__2]])

(1.21)

where epsilon[A] is the corresponding polarization vector.

When working with non-Abelian gauge fields, the interaction Lagrangian involves derivatives. FeynmanDiagrams  can handle that kind of interaction in momentum representation. Consider for instance a Yang-Mills theory with a massless field B[mu, a] where a is a SU2 index (see eq.(12) of sec. 19.4 of ref.[1]). The interaction Lagrangian can be entered as follows

Setup(su2indices = lowercaselatin_ah, massless = B, op = B)

`* Partial match of  '`*massless*`' against keyword '`*masslessfields*`' `

 

`* Partial match of  '`*op*`' against keyword '`*quantumoperators*`' `

 

_______________________________________________________

 

[masslessfields = {A, B}, quantumoperators = {A, B, phi, psi, psi1}, su2indices = lowercaselatin_ah]

(1.22)

Define(B[mu, a], quiet)

F__B[mu, nu, a] := d_[mu](B[nu, a](X))-d_[nu](B[mu, a](X))

Physics:-d_[mu](B[nu, a](X), [X])-Physics:-d_[nu](B[mu, a](X), [X])

(1.23)

L := (1/2)*g*LeviCivita[a, b, c]*F__B[mu, nu, a]*B[mu, b](X)*B[nu, c](X)+(1/4)*g^2*LeviCivita[a, b, c]*LeviCivita[a, e, f]*B[mu, b](X)*B[nu, c](X)*B[mu, e](X)*B[nu, f](X)

(1/2)*g*Physics:-LeviCivita[a, b, c]*Physics:-`*`(Physics:-d_[mu](B[nu, a](X), [X])-Physics:-d_[nu](B[mu, a](X), [X]), B[`~mu`, b](X), B[`~nu`, c](X))+(1/4)*g^2*Physics:-LeviCivita[a, b, c]*Physics:-LeviCivita[a, e, f]*Physics:-`*`(B[mu, b](X), B[nu, c](X), B[`~mu`, e](X), B[`~nu`, f](X))

(1.24)

The transition probability density at tree-level for a process with two incoming and two outgoing B particles is given by

FeynmanDiagrams(L, incomingparticles = [B, B], outgoingparticles = [B, B], numberofloops = 0, output = probabilitydensity, factor, diagrams)

`* Partial match of  '`*factor*`' against keyword '`*factortreelevel*`' `

(1.25)

 

 

Physics:-FeynmanDiagrams:-ProbabilityDensity(4*Pi^2*%mul(n[i], i = 1 .. 2)*abs(F)^2*Dirac(-P__3[`~sigma`]-P__4[`~sigma`]+P__1[`~sigma`]+P__2[`~sigma`])*%mul(dP_[f]^3, f = 1 .. 2), F = (((1/8)*I)*Physics:-LeviCivita[a1, a3, h]*((-P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics:-g_[`~lambda`, `~tau`]+(P__1[`~lambda`]+P__2[`~lambda`]+P__3[`~lambda`])*Physics:-g_[`~kappa`, `~tau`]-Physics:-g_[`~kappa`, `~lambda`]*(P__3[`~tau`]-P__4[`~tau`]))*Physics:-LeviCivita[a2, d, g]*((P__1[`~beta`]+(1/2)*P__2[`~beta`])*Physics:-g_[`~alpha`, `~sigma`]+(-(1/2)*P__1[`~sigma`]+(1/2)*P__2[`~sigma`])*Physics:-g_[`~alpha`, `~beta`]-(1/2)*Physics:-g_[`~beta`, `~sigma`]*(P__1[`~alpha`]+2*P__2[`~alpha`]))*Physics:-g_[sigma, tau]*Physics:-KroneckerDelta[a2, a3]/((-P__1[chi]-P__2[chi])*(-P__1[`~chi`]-P__2[`~chi`])+I*Physics:-FeynmanDiagrams:-epsilon)-((1/16)*I)*((-P__1[`~beta`]+P__3[`~beta`]-P__4[`~beta`])*Physics:-g_[`~lambda`, `~tau`]+(P__1[`~lambda`]-P__2[`~lambda`]-P__3[`~lambda`])*Physics:-g_[`~beta`, `~tau`]+Physics:-g_[`~beta`, `~lambda`]*(P__2[`~tau`]+P__4[`~tau`]))*Physics:-LeviCivita[a1, a3, g]*((P__1[`~sigma`]+P__3[`~sigma`])*Physics:-g_[`~alpha`, `~kappa`]+(-2*P__1[`~kappa`]+P__3[`~kappa`])*Physics:-g_[`~alpha`, `~sigma`]+Physics:-g_[`~kappa`, `~sigma`]*(P__1[`~alpha`]-2*P__3[`~alpha`]))*Physics:-LeviCivita[a2, d, h]*Physics:-g_[sigma, tau]*Physics:-KroneckerDelta[a2, a3]/((-P__1[chi]+P__3[chi])*(-P__1[`~chi`]+P__3[`~chi`])+I*Physics:-FeynmanDiagrams:-epsilon)-((1/16)*I)*((-P__1[`~beta`]-P__3[`~beta`]+P__4[`~beta`])*Physics:-g_[`~kappa`, `~tau`]+(P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics:-g_[`~beta`, `~tau`]+Physics:-g_[`~beta`, `~kappa`]*(P__2[`~tau`]+P__3[`~tau`]))*Physics:-LeviCivita[a3, g, h]*((P__1[`~sigma`]+P__4[`~sigma`])*Physics:-g_[`~alpha`, `~lambda`]+(P__1[`~alpha`]-2*P__4[`~alpha`])*Physics:-g_[`~lambda`, `~sigma`]-2*Physics:-g_[`~alpha`, `~sigma`]*(P__1[`~lambda`]-(1/2)*P__4[`~lambda`]))*Physics:-LeviCivita[a1, a2, d]*Physics:-g_[sigma, tau]*Physics:-KroneckerDelta[a2, a3]/((-P__1[chi]+P__4[chi])*(-P__1[`~chi`]+P__4[`~chi`])+I*Physics:-FeynmanDiagrams:-epsilon)-((1/16)*I)*(Physics:-KroneckerDelta[g, h]*Physics:-KroneckerDelta[a1, d]*(Physics:-g_[`~alpha`, `~beta`]*Physics:-g_[`~kappa`, `~lambda`]+Physics:-g_[`~alpha`, `~kappa`]*Physics:-g_[`~beta`, `~lambda`]-2*Physics:-g_[`~alpha`, `~lambda`]*Physics:-g_[`~beta`, `~kappa`])+Physics:-KroneckerDelta[d, h]*(Physics:-g_[`~alpha`, `~beta`]*Physics:-g_[`~kappa`, `~lambda`]-2*Physics:-g_[`~alpha`, `~kappa`]*Physics:-g_[`~beta`, `~lambda`]+Physics:-g_[`~alpha`, `~lambda`]*Physics:-g_[`~beta`, `~kappa`])*Physics:-KroneckerDelta[a1, g]-2*(Physics:-g_[`~alpha`, `~beta`]*Physics:-g_[`~kappa`, `~lambda`]-(1/2)*Physics:-g_[`~beta`, `~kappa`]*Physics:-g_[`~alpha`, `~lambda`]-(1/2)*Physics:-g_[`~alpha`, `~kappa`]*Physics:-g_[`~beta`, `~lambda`])*Physics:-KroneckerDelta[d, g]*Physics:-KroneckerDelta[a1, h]))*g^2*conjugate(Physics:-FeynmanDiagrams:-PolarizationVector[B][kappa, h](P__3_))*conjugate(Physics:-FeynmanDiagrams:-PolarizationVector[B][lambda, a1](P__4_))*Physics:-FeynmanDiagrams:-PolarizationVector[B][alpha, d](P__1_)*Physics:-FeynmanDiagrams:-PolarizationVector[B][beta, g](P__2_)/(Pi^2*(E__1*E__2*E__3*E__4)^(1/2)))

(1.26)

To simplify the repeated indices, us the option simplifytensorindices. To check the indices entering a result like this one use Check ; there are no free indices, and regarding the repeated indices:

Check(Physics[FeynmanDiagrams]:-ProbabilityDensity(4*Pi^2*%mul(n[i], i = 1 .. 2)*abs(F)^2*Dirac(-P__3[`~sigma`]-P__4[`~sigma`]+P__1[`~sigma`]+P__2[`~sigma`])*%mul(dP_[f]^3, f = 1 .. 2), F = (((1/8)*I)*Physics[LeviCivita][a1, a3, h]*((-P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics[g_][`~lambda`, `~tau`]+(P__1[`~lambda`]+P__2[`~lambda`]+P__3[`~lambda`])*Physics[g_][`~kappa`, `~tau`]-Physics[g_][`~kappa`, `~lambda`]*(P__3[`~tau`]-P__4[`~tau`]))*Physics[LeviCivita][a2, d, g]*((P__1[`~beta`]+(1/2)*P__2[`~beta`])*Physics[g_][`~alpha`, `~sigma`]+(-(1/2)*P__1[`~sigma`]+(1/2)*P__2[`~sigma`])*Physics[g_][`~alpha`, `~beta`]-(1/2)*Physics[g_][`~beta`, `~sigma`]*(P__1[`~alpha`]+2*P__2[`~alpha`]))*Physics[g_][sigma, tau]*Physics[KroneckerDelta][a2, a3]/((-P__1[chi]-P__2[chi])*(-P__1[`~chi`]-P__2[`~chi`])+I*Physics[FeynmanDiagrams]:-epsilon)-((1/16)*I)*((-P__1[`~beta`]+P__3[`~beta`]-P__4[`~beta`])*Physics[g_][`~lambda`, `~tau`]+(P__1[`~lambda`]-P__2[`~lambda`]-P__3[`~lambda`])*Physics[g_][`~beta`, `~tau`]+Physics[g_][`~beta`, `~lambda`]*(P__2[`~tau`]+P__4[`~tau`]))*Physics[LeviCivita][a1, a3, g]*((P__1[`~sigma`]+P__3[`~sigma`])*Physics[g_][`~alpha`, `~kappa`]+(-2*P__1[`~kappa`]+P__3[`~kappa`])*Physics[g_][`~alpha`, `~sigma`]+Physics[g_][`~kappa`, `~sigma`]*(P__1[`~alpha`]-2*P__3[`~alpha`]))*Physics[LeviCivita][a2, d, h]*Physics[g_][sigma, tau]*Physics[KroneckerDelta][a2, a3]/((-P__1[chi]+P__3[chi])*(-P__1[`~chi`]+P__3[`~chi`])+I*Physics[FeynmanDiagrams]:-epsilon)-((1/16)*I)*((-P__1[`~beta`]-P__3[`~beta`]+P__4[`~beta`])*Physics[g_][`~kappa`, `~tau`]+(P__1[`~kappa`]-P__2[`~kappa`]-P__4[`~kappa`])*Physics[g_][`~beta`, `~tau`]+Physics[g_][`~beta`, `~kappa`]*(P__2[`~tau`]+P__3[`~tau`]))*Physics[LeviCivita][a3, g, h]*((P__1[`~sigma`]+P__4[`~sigma`])*Physics[g_][`~alpha`, `~lambda`]+(P__1[`~alpha`]-2*P__4[`~alpha`])*Physics[g_][`~lambda`, `~sigma`]-2*Physics[g_][`~alpha`, `~sigma`]*(P__1[`~lambda`]-(1/2)*P__4[`~lambda`]))*Physics[LeviCivita][a1, a2, d]*Physics[g_][sigma, tau]*Physics[KroneckerDelta][a2, a3]/((-P__1[chi]+P__4[chi])*(-P__1[`~chi`]+P__4[`~chi`])+I*Physics[FeynmanDiagrams]:-epsilon)-((1/16)*I)*(Physics[KroneckerDelta][g, h]*Physics[KroneckerDelta][a1, d]*(Physics[g_][`~alpha`, `~beta`]*Physics[g_][`~kappa`, `~lambda`]+Physics[g_][`~alpha`, `~kappa`]*Physics[g_][`~beta`, `~lambda`]-2*Physics[g_][`~alpha`, `~lambda`]*Physics[g_][`~beta`, `~kappa`])+Physics[KroneckerDelta][d, h]*(Physics[g_][`~alpha`, `~beta`]*Physics[g_][`~kappa`, `~lambda`]-2*Physics[g_][`~alpha`, `~kappa`]*Physics[g_][`~beta`, `~lambda`]+Physics[g_][`~alpha`, `~lambda`]*Physics[g_][`~beta`, `~kappa`])*Physics[KroneckerDelta][a1, g]-2*(Physics[g_][`~alpha`, `~beta`]*Physics[g_][`~kappa`, `~lambda`]-(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[g_][`~beta`, `~kappa`]-(1/2)*Physics[g_][`~alpha`, `~kappa`]*Physics[g_][`~beta`, `~lambda`])*Physics[KroneckerDelta][d, g]*Physics[KroneckerDelta][a1, h]))*g^2*conjugate(Physics[FeynmanDiagrams]:-PolarizationVector[B][kappa, h](P__3_))*conjugate(Physics[FeynmanDiagrams]:-PolarizationVector[B][lambda, a1](P__4_))*Physics[FeynmanDiagrams]:-PolarizationVector[B][alpha, d](P__1_)*Physics[FeynmanDiagrams]:-PolarizationVector[B][beta, g](P__2_)/(Pi^2*(E__1*E__2*E__3*E__4)^(1/2))), all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`, the free indices are: `*{`...`}

 

[{a1, a2, a3, alpha, beta, chi, d, g, h, kappa, lambda, sigma, tau}], {}

(1.27)


This process can be computed with 1 or more loops, in which case the number of terms increases significantly. As another interesting non-Abelian model, consider the interaction Lagrangian of the electro-weak part of the Standard Model

Coordinates(clear, Z)

`Unaliasing `*{Z}*` previously defined as a system of spacetime coordinates`

(1.28)

Setup(quantumoperators = {W, Z})

[quantumoperators = {A, B, W, Z, phi, psi, psi1}]

(1.29)

Define(W[mu], Z[mu])

`Defined objects with tensor properties`

 

{A[mu], B[mu, a], Physics:-Dgamma[mu], P__1[mu], P__2[mu], P__3[alpha], P__4[alpha], Physics:-Psigma[mu], W[mu], Z[mu], Physics:-d_[mu], Physics:-g_[mu, nu], p__1[mu], p__2[mu], p__3[mu], p__4[mu], p__5[mu], psi[j], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X), Physics:-SpaceTimeVector[mu](Y)}

(1.30)

CompactDisplay((W, Z)(X))

` W`(X)*`will now be displayed as`*W

 

` Z`(X)*`will now be displayed as`*Z

(1.31)

F__W[mu, nu] := d_[mu](W[nu](X))-d_[nu](W[mu](X))

Physics:-d_[mu](W[nu](X), [X])-Physics:-d_[nu](W[mu](X), [X])

(1.32)

F__Z[mu, nu] := d_[mu](Z[nu](X))-d_[nu](Z[mu](X))

Physics:-d_[mu](Z[nu](X), [X])-Physics:-d_[nu](Z[mu](X), [X])

(1.33)

L__WZ := I*g*cos(`&theta;__w`)*((Dagger(F__W[mu, nu])*W[mu](X)-Dagger(W[mu](X))*F__W[mu, nu])*Z[nu](X)+W[nu](X)*Dagger(W[mu](X))*F__Z[mu, nu])

I*g*cos(theta__w)*(Physics:-`*`(Physics:-`*`(Physics:-d_[mu](Physics:-Dagger(W[nu](X)), [X])-Physics:-d_[nu](Physics:-Dagger(W[mu](X)), [X]), W[`~mu`](X))-Physics:-`*`(Physics:-Dagger(W[mu](X)), Physics:-d_[`~mu`](W[nu](X), [X])-Physics:-d_[nu](W[`~mu`](X), [X])), Z[`~nu`](X))+Physics:-`*`(W[nu](X), Physics:-Dagger(W[mu](X)), Physics:-d_[`~mu`](Z[`~nu`](X), [X])-Physics:-d_[`~nu`](Z[`~mu`](X), [X])))

(1.34)

This interaction Lagrangian contains six different terms. The S-matrix element for the tree-level process with two incoming and two outgoing W particles is shown in the help page for FeynmanDiagrams .

NULL

References

 

[1] Bogoliubov, N.N., and Shirkov, D.V. Quantum Fields. Benjamin Cummings, 1982.

[2] Weinberg, S., The Quantum Theory Of Fields. Cambridge University Press, 2005.

 

FeynmanDiagrams_and_the_Scattering_Matrix.PDF

FeynmanDiagrams_and_the_Scattering_Matrix.mw

FeynmanDiagrams_-_help_page.mw


Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

with(plots):
S:=cat("Happy New Year 2020!   "$3):
N:=length(S): a:=0.77*Pi: h:=2*Pi/N:
display(seq(textplot([cos(a-k*h), sin(a-k*h),S[k+1]], 
        rotation=-Pi/2+a-k*h, 'font'=["times","roman",24]), k=0..N-4), axes=none);

When plotted, these parametric equations say "happy new year" (and were constructed with this worksheet)

x := piecewise(t <= 58, -15.0*sin(1.43 + 0.650*t) - 14.8*sin(-1.64 + 0.703*t) - 1.28*sin(-2.97 + 1.25*t) - 11.9*sin(-1.58 + 0.540*t) - 1.07*sin(-1.60 + 1.35*t) - 3.85*sin(-2.09 + 1.41*t) - 7.13*sin(1.13 + 1.73*t) - 4.40*sin(1.32 + 1.30*t) - 26.3*sin(1.53 + 0.380*t) - 9.42*sin(-4.65 + 0.433*t) - 3.43*sin(1.42 + 2.06*t) - 7.57*sin(-1.77 + 2.11*t) - 2.65*sin(-4.34 + 0.323*t) - 1.95*sin(-4.57 + 2.54*t) - 5.39*sin(-1.38 + 2.60*t) - 49.2*sin(1.52 + 0.487*t) - 0.754*sin(-4.38 + 2.87*t) - 9.67*sin(-1.58 + 2.65*t) - 7.88*sin(-4.59 + 1.95*t) - 2.39*sin(-1.67 + 2.71*t) - 15.1*sin(1.53 + 0.108*t) - 39.0*sin(1.47 + 0.757*t) - 1.80*sin(1.37 + 2.22*t) - 4.22*sin(-1.95 + 0.973*t) - 7.72*sin(-1.44 + 2.17*t) - 8.80*sin(-1.66 + 0.813*t) - 3.59*sin(1.13 + 1.57*t) - 15.4*sin(-1.64 + 1.62*t) - 6.70*sin(1.36 + 1.19*t) - 791.*sin(-1.57 + 0.0540*t) - 2.55*sin(-1.55 + 1.89*t) - 6.92*sin(-1.87 + 1.68*t) - 3.95*sin(1.17 + 1.08*t) - 44.1*sin(-1.67 + 1.14*t) - 25.8*sin(1.51 + 0.597*t) - 31.4*sin(1.42 + 1.46*t) - 96.8*sin(-1.59 + 0.162*t) - 18.7*sin(1.53 + 0.217*t) - 7.87*sin(-4.66 + 2.98*t) - 4.99*sin(1.22 + 3.03*t) - 6.92*sin(1.43 + 2.44*t) - 48.3*sin(1.47 + 1.03*t) - 24.2*sin(1.48 + 1.52*t) - 9.58*sin(1.43 + 2.49*t) - 4.29*sin(1.33 + 2.27*t) - 6.34*sin(1.22 + 2.33*t) - 12.0*sin(1.45 + 2.00*t) - 0.388*sin(-1.25 + 2.92*t) - 2.74*sin(-1.43 + 1.79*t) - 6.71*sin(-1.66 + 1.84*t) - 0.713*sin(-3.63 + 2.38*t) - 43.1*sin(-1.59 + 0.271*t) - 2.51*sin(1.12 + 2.76*t) - 1.29*sin(-3.92 + 2.82*t) - 21.3*sin(-1.70 + 0.867*t) - 12.4*sin(1.50 + 0.920*t), 58 < t and t <= 84, -500 - 321.*sin(-8.608 + 0.121*t) - 7.18*sin(-12.408 + 0.241*t) - 57.1*sin(-22.608 + 0.361*t) - 21.9*sin(-26.682 + 0.484*t) - 21.3*sin(-33.474 + 0.603*t) - 50.2*sin(-43.800 + 0.725*t) - 20.6*sin(-50.760 + 0.845*t) - 41.5*sin(-54.756 + 0.967*t) - 9.74*sin(-61.89 + 1.09*t) - 41.1*sin(-72.03 + 1.21*t) - 2.49*sin(-78.88 + 1.33*t) - 3.30*sin(-83.227 + 1.45*t) - 6.73*sin(-89.99 + 1.57*t) - 5.88*sin(-96.59 + 1.69*t) - 16.4*sin(-106.99 + 1.81*t) - 1.61*sin(-111.8982 + 1.93*t) - 1.84*sin(-117.970 + 2.05*t) - 0.464*sin(-127.83 + 2.17*t) - 1.64*sin(-134.90 + 2.30*t) - 3.94*sin(-142.37 + 2.41*t) - 2.35*sin(-149.22 + 2.54*t) - 2.72*sin(-154.3362 + 2.66*t) - 8.41*sin(-160.453 + 2.78*t) - 4.39*sin(-171.17 + 2.90*t), 84 < t, -300 - 2.66*sin(-205.04 + 2.41*t) - 1.26*sin(-207.397 + 2.46*t) - 2.21*sin(-196.59 + 2.31*t) - 2.31*sin(-166.83 + 1.96*t) - 48.9*sin(-39.688 + 0.452*t) - 0.697*sin(-252.158 + 3.01*t) - 2.51*sin(-179.22 + 2.11*t) - 1.57*sin(-222.14 + 2.66*t) - 0.745*sin(-226.24 + 2.71*t) - 49.4*sin(-10.020 + 0.100*t) - 0.289*sin(-159.628 + 1.91*t) - 95.9*sin(-32.358 + 0.402*t) - 60.0*sin(-43.928 + 0.502*t) - 3.76*sin(-73.736 + 0.854*t) - 3.05*sin(-183.97 + 2.16*t) - 0.629*sin(-158.50 + 1.86*t) - 9.25*sin(-49.272 + 0.603*t) - 4.46*sin(-74.716 + 0.904*t) - 10.4*sin(-79.040 + 0.955*t) - 2.65*sin(-103.67 + 1.21*t) - 1.99*sin(-145.57 + 1.71*t) - 1.52*sin(-197.315 + 2.36*t) - 0.685*sin(-258.12 + 3.06*t) - 1.04*sin(-247.58 + 2.96*t) - 64.8*sin(-18.514 + 0.201*t) - 68.5*sin(-31.278 + 0.352*t) - 579.*sin(-5.8068 + 0.0502*t) - 6.52*sin(-95.20 + 1.11*t) - 5.03*sin(-96.28 + 1.16*t) - 0.396*sin(-211.620 + 2.51*t) - 7.28*sin(-150.00 + 1.76*t) - 2.42*sin(-153.92 + 1.81*t) - 10.4*sin(-112.11 + 1.31*t) - 24.8*sin(-85.95 + 1.00*t) - 3.91*sin(-124.83 + 1.46*t) - 1.69*sin(-185.369 + 2.21*t) - 1.18*sin(-189.238 + 2.26*t) - 16.6*sin(-56.662 + 0.653*t) - 1.33*sin(-222.31 + 2.61*t) - 0.593*sin(-238.70 + 2.81*t) - 1.88*sin(-233.58 + 2.76*t) - 3.91*sin(-133.01 + 1.56*t) - 4.94*sin(-134.16 + 1.61*t) - 9.59*sin(-128.89 + 1.51*t) - 1.02*sin(-240.2714 + 2.86*t) - 2.15*sin(-247.83 + 2.91*t) - 5.52*sin(-90.85 + 1.06*t) - 3.83*sin(-171.25 + 2.01*t) - 0.523*sin(-171.66 + 2.06*t) - 0.284*sin(-141.80 + 1.66*t) - 23.2*sin(-11.174 + 0.151*t) - 1.58*sin(-114.615 + 1.36*t) - 2.67*sin(-120.75 + 1.41*t) - 5.83*sin(-19.524 + 0.251*t) - 13.7*sin(-23.774 + 0.301*t) - 14.8*sin(-107.89 + 1.26*t) - 15.5*sin(-60.842 + 0.703*t) - 37.7*sin(-65.176 + 0.754*t) - 2.02*sin(-217.95 + 2.56*t) - 13.2*sin(-69.466 + 0.804*t) - 37.7*sin(-45.052 + 0.553*t)):

y := piecewise(t <= 58, -28.1*sin(1.45 + 1.62*t) - 2.23*sin(-2.39 + 1.89*t) - 17.8*sin(-1.51 + 1.19*t) - 4.85*sin(-1.61 + 2.38*t) - 2.52*sin(1.55 + 1.95*t) - 20.0*sin(1.55 + 2.11*t) - 24.8*sin(-1.62 + 2.00*t) - 19.9*sin(-1.81 + 2.06*t) - 4.22*sin(-0.422 + 2.60*t) - 6.94*sin(1.47 + 2.87*t) - 61.1*sin(1.49 + 0.323*t) - 13.9*sin(-4.68 + 0.540*t) - 3.97*sin(0.00256 + 2.33*t) - 69.8*sin(1.53 + 0.487*t) - 59.6*sin(1.50 + 0.813*t) - 132.*sin(-1.65 + 0.867*t) - 26.7*sin(-1.76 + 1.52*t) - 53.1*sin(1.40 + 1.57*t) - 139.*sin(1.57 + 0.0540*t) - 3.75*sin(-2.34 + 3.03*t) - 8.03*sin(1.24 + 1.73*t) - 22.9*sin(-4.61 + 0.217*t) - 16.7*sin(-1.67 + 0.703*t) - 23.3*sin(-1.82 + 1.68*t) - 78.9*sin(-4.70 + 0.271*t) - 2.72*sin(-2.38 + 2.49*t) - 3.45*sin(1.10 + 2.54*t) - 2.07*sin(-0.489 + 2.22*t) - 13.1*sin(-1.82 + 2.27*t) - 60.6*sin(-1.62 + 1.08*t) - 5.27*sin(1.55 + 2.44*t) - 4.17*sin(1.46 + 2.82*t) - 33.1*sin(-1.80 + 1.46*t) - 2.15*sin(-1.58 + 0.757*t) - 3.94*sin(-3.86 + 2.65*t) - 8.88*sin(1.51 + 1.79*t) - 9.97*sin(1.52 + 1.84*t) - 105.*sin(1.48 + 1.03*t) - 15.2*sin(-4.67 + 1.25*t) - 101.*sin(1.51 + 0.380*t) - 11.0*sin(-4.59 + 0.433*t) - 86.7*sin(1.50 + 0.973*t) - 170.*sin(1.53 + 0.597*t) - 41.2*sin(1.51 + 0.650*t) - 20.4*sin(-1.67 + 1.30*t) - 47.9*sin(-1.70 + 1.35*t) - 15.8*sin(-1.66 + 2.71*t) - 8.61*sin(-1.71 + 2.76*t) - 25.7*sin(-1.64 + 0.108*t) - 70.9*sin(1.55 + 0.162*t) - 0.668*sin(-2.42 + 2.92*t) - 4.78*sin(-4.60 + 2.98*t) - 106.*sin(1.49 + 0.920*t) - 17.6*sin(1.53 + 1.41*t) - 8.82*sin(1.05 + 2.17*t) - 113.*sin(-1.67 + 1.14*t), t <= 84, -800 - 7.30*sin(-171.17 + 2.90*t) - 3.28*sin(-6.550 + 0.121*t) - 1.46*sin(-17.878 + 0.241*t) - 20.4*sin(-22.438 + 0.361*t) - 28.9*sin(-29.862 + 0.484*t) - 9.13*sin(-36.364 + 0.603*t) - 45.3*sin(-40.650 + 0.725*t) - 97.4*sin(-50.770 + 0.845*t) - 13.1*sin(-54.916 + 0.967*t) - 80.8*sin(-61.97 + 1.09*t) - 39.1*sin(-71.92 + 1.21*t) - 42.8*sin(-78.87 + 1.33*t) - 108.*sin(-85.97 + 1.45*t) - 10.6*sin(-92.80 + 1.57*t) - 49.8*sin(-99.94 + 1.69*t) - 15.4*sin(-103.75 + 1.81*t) - 24.2*sin(-113.90 + 1.93*t) - 8.96*sin(-123.18 + 2.05*t) - 1.59*sin(-127.14 + 2.17*t) - 14.1*sin(-137.59 + 2.30*t) - 6.51*sin(-142.35 + 2.41*t) - 7.98*sin(-145.83 + 2.54*t) - 6.40*sin(-153.721 + 2.66*t) - 1.23*sin(-164.36 + 2.78*t), 84 < t, -1400 - 128.*sin(-32.358 + 0.402*t) - 68.5*sin(-43.928 + 0.502*t) - 2.55*sin(-242.18 + 2.86*t) - 6.86*sin(-219.136 + 2.61*t) - 5.76*sin(-222.904 + 2.66*t) - 2.39*sin(-226.835 + 2.71*t) - 101.*sin(-11.164 + 0.151*t) - 8.69*sin(-231.548 + 2.76*t) - 146.*sin(-31.268 + 0.352*t) - 8.30*sin(-179.37 + 2.11*t) - 2.68*sin(-261.69 + 3.06*t) - 10.4*sin(-162.98 + 1.91*t) - 30.1*sin(-73.606 + 0.854*t) - 24.1*sin(-77.946 + 0.904*t) - 10.0*sin(-146.01 + 1.71*t) - 72.5*sin(-69.416 + 0.804*t) - 8.91*sin(-85.97 + 1.00*t) - 8.58*sin(-175.51 + 2.06*t) - 27.4*sin(-109.01 + 1.31*t) - 16.8*sin(-113.17 + 1.36*t) - 162.*sin(-5.7968 + 0.0502*t) - 3.69*sin(-205.52 + 2.41*t) - 7.62*sin(-207.006 + 2.46*t) - 131.*sin(-53.522 + 0.653*t) - 95.3*sin(-60.882 + 0.703*t) - 8.53*sin(-197.627 + 2.36*t) - 1.74*sin(-247.32 + 2.91*t) - 27.2*sin(-121.51 + 1.46*t) - 51.7*sin(-49.332 + 0.603*t) - 8.81*sin(-104.925 + 1.26*t) - 10.2*sin(-100.703 + 1.21*t) - 9.35*sin(-183.90 + 2.16*t) - 7.82*sin(-188.20 + 2.21*t) - 42.8*sin(-26.964 + 0.301*t) - 16.8*sin(-48.312 + 0.553*t) - 15.2*sin(-9.980 + 0.100*t) - 213.*sin(-18.524 + 0.201*t) - 39.4*sin(-19.584 + 0.251*t) - 6.28*sin(-87.85 + 1.06*t) - 3.71*sin(-117.623 + 1.41*t) - 4.92*sin(-196.77 + 2.31*t) - 1.25*sin(-255.21 + 3.01*t) - 5.13*sin(-248.529 + 2.96*t) - 8.69*sin(-141.43 + 1.66*t) - 11.5*sin(-167.26 + 1.96*t) - 13.0*sin(-171.19 + 2.01*t) - 4.12*sin(-159.23 + 1.86*t) - 3.66*sin(-212.23 + 2.51*t) - 0.810*sin(-83.380 + 0.955*t) - 3.11*sin(-65.516 + 0.754*t) - 1.38*sin(-139.34 + 1.61*t) - 9.07*sin(-188.885 + 2.26*t) - 52.6*sin(-39.678 + 0.452*t) - 6.81*sin(-125.917 + 1.51*t) - 24.7*sin(-130.128 + 1.56*t) - 4.16*sin(-215.362 + 2.56*t) - 11.8*sin(-92.283 + 1.11*t) - 16.6*sin(-96.32 + 1.16*t) - 6.39*sin(-147.108 + 1.76*t) - 7.61*sin(-154.46 + 1.81*t) - 4.28*sin(-235.566 + 2.81*t)):

plot( [ x, y, t = 0 .. 146 ], scaling = constrained, discont = [ usefdiscont ], axes = boxed, thickness = 5, size = [600, 600]);

 

Here is a little animation to wish all of you a Merry Christmas

FireWorks.mw


Here we simulate the motion of a container with a flat bottom that can slide on a horizontal surface subject to dry friction (Coulomb friction).  Installed inside the container is an ordinary mass/spring/damper system where the mass slides horizontally.  We impart an initial velocity to the container.  That sets the mass into motion which then affects the container's motion.  Under certain conditions the container will undergo a stick-slip motion which is evident in the simulation.

This simulation very roughly approximates the motion of a partially filled bucket of water that slides on the floor when kicked.  The idea arose in a discussoin with Carl Love and mmcdara:
https://www.mapleprimes.com/posts/211677-Mass-Spring-Conveyor-Belt-And-Dry-Friction

In the animation below, the container is shown in dark color when it slides against the floor, and light color when it sticks.

Worksheet: slosh.mw

 

Here is an animation of a mass-spring system where the mass slides horizontally on a steadily moving conveyor belt.

The contact between the block of mass and the belt is of the dry friction kind (Coulomb's friction). Consequently the block periodically sticks to the belt and moves forward with it until the force of the stretching spring overcomes the force of friction and yanks it back, making it to slip against the belt. In the animation the block is shown in a dark color while slipping, and a light color while sticking.

The fully executed Maple worksheet can be slow to load and requires a good deal of memory. Therefore I have attached two versions which are identical in all respects except that in one of them I have removed the Maple output to make is easy to load if your computer has limitations.

Download worksheet (no Maple output) block-sliding-on-conveyor-belt-stripped.mw

Doiwnload worksheet (with Maple output) (sorry, exceeds MaplePrime's size limit)

Here is two solutions with Maple of the problem A2 of  Putnam Mathematical Competition 2019 . The first solution is entirely based on the use of the  geometry  package; the second solution does not use this package. Since the triangle is defined up to similarity, without loss of generality, we can set its vertices  A(0,0) , B(1,0) , C(x0,y0)  and then calculate the parameters  x0, y0  using the conditions of the problem. 


 

The problem

A2: In the triangle ∆ABC, let G be the centroid, and let I be the center of the
inscribed circle. Let α and β be the angles at the vertices A and B, respectively.
Suppose that the segment IG is parallel to AB and that  β = 2*arctan(1/3).  Find α.
 

# Solution 1 with the geometry package
restart;

# Calculation

with(geometry):
local I:
point(A,0,0): point(B,1,0): point(C,x0,y0):
assume(y0>0,-y0*(-1+x0-((1-x0)^2+y0^2)^(1/2))+y0*((x0^2+y0^2)^(1/2)+x0) <> 0):
triangle(t,[A,B,C]):
incircle(ic,t, 'centername'=I):
Cn:=coordinates(I):
centroid(G,t):
CG:=coordinates(G):
a:=-expand(tan(2*arctan(1/3))):
solve({Cn[2]=CG[2],y0/(x0-1)=a}, explicit);
point(C,eval([x0,y0],%)):
answer=FindAngle(line(AB,[A,B]),line(AC,[A,C]));

# Visualization (G is the point of medians intersection)

triangle(t,[A,B,C]):
incircle(ic,t, 'centername'=I):
centroid(G,t):
segment(s,[I,G]):
median(mB,B,t): median(mC,C,t):
draw([A(symbol=solidcircle,color=black),B(symbol=solidcircle,color=black),C(symbol=solidcircle,color=black),I(symbol=solidcircle,color=green),G(symbol=solidcircle,color=blue),t(color=black),s(color=red,thickness=2),ic(color=green),mB(color=blue,thickness=0),mC(color=blue,thickness=0)], axes=none, size=[800,500], printtext=true,font=[times,20]);

I

 

Warning, The imaginary unit, I, has been renamed _I

 

Warning, solve may be ignoring assumptions on the input variables.

 

{x0 = 0, y0 = 3/4}

 

answer = (1/2)*Pi

 

 


# Solution 2 by a direct calculation

# Calculation

restart;
local I;
sinB:=y0/sqrt(x0^2+y0^2):
cosB:=x0/sqrt(x0^2+y0^2):
Sol1:=eval([x,y],solve({y=-(x-1)/3,y=(sinB/(1+cosB))*x}, {x,y})):
tanB:=expand(tan(2*arctan(1/3))):
Sol2:=solve({y0/3=Sol1[2],y0=-tanB*(x0-1)},explicit);
A:=[0,0]: B:=[1,0]: C:=eval([x0,y0],Sol2[2]):
AB:=<(B-A)[]>: AC:=<(C-A)[]>:
answer=arccos(AB.AC/sqrt(AB.AB)/sqrt(AC.AC));

# Visualization

with(plottools): with(plots):
ABC:=curve([A,B,C,A]):
I:=simplify(eval(Sol1,Sol2[2]));
c:=circle(I,eval(Sol1[2],Sol2[2]),color=green):
G:=(A+B+C)/~3;
IG:=line(I,G,color=red,thickness=2):
P:=pointplot([A,B,C,I,G], color=[black$3,green,blue], symbol=solidcircle):
T:=textplot([[A[],"A"],[B[],"B"],[C[],"C"],[I[],"I"],[G[],"G"]], font=[times,20], align=[left,below]):
M:=plot([[(C+t*~((A+B)/2-C))[],t=0..1],[(B+t*~((A+C)/2-B))[],t=0..1]], color=blue, thickness=0):
display(ABC,c,IG,P,T,M, scaling=constrained, axes=none,size=[800,500]);

I

 

Warning, The imaginary unit, I, has been renamed _I

 

{x0 = 1, y0 = 0}, {x0 = 0, y0 = 3/4}

 

answer = (1/2)*Pi

 

[1/4, 1/4]

 

[1/3, 1/4]

 

 

 


 

Download Putnam.mw

 

Maple can easily solve the B4 problem of the Putnam Mathematical Competition 2019  link

 

B4.  Let F be the set of functions f(x,y) that are twice continuously differentiable for x≥1, y≥1 and that satisfy the following two equations:
    x*(diff(f(x, y), x))+y*(diff(f(x, y), y)) = x*y*ln(x*y)

x^2*(diff(f(x, y), x, x))+y^2*(diff(f(x, y), y, y)) = x*y

 

For each f2F, let

 

"m(f) = min[s>=1]  (f(s+1,s+1)-f(s+1,s)-f(s,s+1)+f(s,s))"

 

Determine m(f), and show that it is independent of the choice of f.


 

# Solution

pdsolve({
x*diff(f(x,y),x)+y*diff(f(x,y),y) = x*y*ln(x*y),
x^2*diff(f(x,y),x,x)+y^2*diff(f(x,y),y,y) = x*y
});

{f(x, y) = (1/2)*(x*y+2*_C1)*ln(x*y)-(1/2)*x*y-2*_C1*ln(x)+_C2}

(1)

f:=unapply(rhs(%[]), x,y);

proc (x, y) options operator, arrow; (1/2)*(y*x+2*_C1)*ln(y*x)-(1/2)*y*x-2*_C1*ln(x)+_C2 end proc

(2)

h := f(s+1, s+1) - f(s+1, s) - f(s, s+1) + f(s, s);

(1/2)*((s+1)^2+2*_C1)*ln((s+1)^2)-(1/2)*(s+1)^2-(s*(s+1)+2*_C1)*ln(s*(s+1))+s*(s+1)+(1/2)*(s^2+2*_C1)*ln(s^2)-(1/2)*s^2

(3)

minimize(h, s=1..infinity);

(4+2*_C1)*ln(2)-1/2-(2+2*_C1)*ln(2)

(4)

answer = simplify(%);

answer = 2*ln(2)-1/2

(5)

 


Download putnam2019-b4.mw

Hi, 

As an amusement,  I decided several months ago to develop some procedures to fill a simple polygon* by hatches or simple textures.

* A simple polygon is a polygon  whose sides either do not intersect or either have a common vertex.

This work started with the very simple observation that if I was capable to hatch or texture a triangle, I will be too to hatch or texture any simple polygon once triangulated.
I also did some work to extend this work to non-simple polygons but there remains some issues to fix (which explains while it is not deliverd here).

The main ideat I used for hatching and texturing is based upon the description of each triangles by a set of 3 inequalities that any interior point must verify.
A hatch of this triangle is thius a segment whose each point is interior.
The closely related idea is used for texturing. Given a simple polygon, periodically replicated to form the texture, the set of points of each replicate that are interior to a given triangle must verify a set of inequalities (the 3 that describe the triangle, plus N if the pattern of the texture is a simple polygon with N sides).

Unfortunately I never finalise this work.
Recently @Christian Wolinski asked a question about texturing that reminded me this "ancient" work of mine.
So I decided to post it as it is, programatically imperfect, lengthy to run, and most of all french-written for a large part.
I guess it is a quite unorthodox way to proceed but some here could be interested by this work to take it back and improve it.

The module named "trianguler" (= triangulate) is a translation into Maple of Frederic Legrand's Python code (full reference given in the worksheet).
I added my own procedure "hachurer" (= hatching) to this module.
The texturing part is not included in this module for it was still in development.

A lot of improvements can be done that I could list, but I prefer not to be too intrusive in your evaluation of this work. So make your own idea about it and do not hesitate to ask me any informations you might need (in particular translation questions).


PS: this work has been done with Maple 2015.2
 

restart:


Reference: http://www.f-legrand.fr/scidoc/docmml/graphie/geometrie/polygone/polygone.html
                    (in french)
                    reference herein : M. de Berg, O. Cheong, M. van Kreveld, M. Overmars,  
                                                 Computational geometry,  (Springer, 2010)

Direct translation of the Frederic Legrand's Python code


Meaning of the different french terms

voisin_sommet  (n, i, di)
        let L the list [1, ..., n] where n is the number of vertices
        This procedure returns the index of the neighbour (voisin) of the vertex (sommet) i when L is rotated by di

equation_droite  (P0, P1, M)
        Let P0 and P1 two vertices and M an arbitrary point.
        Let (P0, P1) the vector originated at P0 and ending at P1 (idem for (P0, M)) and u__Z the unitary vector in the Z direction.
        This procedure returns (P0, P1) o (P0, M) . u__Z

point_dans_triangle  (triangle, M) P1, P2]
        This procedure returns "true" if point M is within (strictly) the  triangle "triangle" and "false" if not.

sommet_distance_maximale  (polygone, P0, P1, P2, indices)    
        Given a polygon (polygone) threes vertices P0, P1 and P2 and a list of indices , this procedure returns
        the vertex of the polygon "polygone" which verifies: 1/ this vertex is strictly within
        the triangle [P0, P1, P2] and 2/ it is the farthest from side [P1, P2] amid all the vertices that verifies point 1/.
        If there is no such vertex the procedure returns NULL.

sommet_gauche (polygone)
        This procedure returns the index of the leftmost ("gauche" means "left") vertex in polygon "polygone".
        If more than one vertices have the same minimum abscissa value then only the first one is returned.

nouveau_polygone(polygone,i_debut,i_fin)
        This procedure creates a new polygon from index i_debut (could be translated by i_first) to i_end (i_last)

trianguler_polygone_recursif(polygone)
        This procedure recursively divides a polygon in two parts A and B from its leftmost vertex.
         If A (B) is a triangle the list "liste_triangles" (mening "list of triangles") is augmented by A (B);
         if not the procedure executes recursively on A and B.

trianguler_polygone(polygone)
         This procedure triangulates the polygon "polygon"

hachurer(shapes, hatch_angle, hatch_number, hatch_color)
         This procedure generates stes of hatches of different angles, colors and numbers


Limitations:
   1/ "polygone" is a simply connected polygon
   2/  two different sides S and S', either do not intersect or either have a common vertex

trianguler := module()
export voisin_sommet, equation_droite, interieur_forme, point_dans_triangle, sommet_distance_maximale,
       sommet_gauche, nouveau_polygone, trianguler_polygone_recursif, trianguler_polygone, hachurer:

#-------------------------------------------------------------------
voisin_sommet := (n, i, di) -> ListTools:-Rotate([$1..n], di)[i]:



#-------------------------------------------------------------------
equation_droite := proc(P0, P1, M) (P1[1]-P0[1])*(M[2]-P0[2]) - (P1[2]-P0[2])*(M[1]-P0[1]) end proc:


#-------------------------------------------------------------------
interieur_forme := proc(forme, M)
  local N:
  N := numelems(forme);
  { seq( equation_droite(forme[n], forme[piecewise(n=N, 1, n+1)], M) >= 0, n=1..N) }
end proc:


#-------------------------------------------------------------------
point_dans_triangle := proc(triangle, M)
  `and`(
          is( equation_droite(triangle[1], triangle[2], M) > 0 ),
          is( equation_droite(triangle[2], triangle[3], M) > 0 ),
          is( equation_droite(triangle[3], triangle[1], M) > 0 )
       )
end proc:



#-------------------------------------------------------------------
sommet_distance_maximale := proc(polygone, P0, P1, P2, indices)
  local n, distance, j, i, M, d;

  n        := numelems(polygone):
  distance := 0:
  j        := NULL:

  for i from 1 to n do
    if `not`(member(i, indices)) then
      M := polygone[i];
      if point_dans_triangle([P0, P1, P2], M) then
        d := abs(equation_droite(P1, P2, M)):
        if d > distance then
          distance := d:
          j := i
        end if:
      end if:
    end if:
  end do:
  return j:
end proc:


#-------------------------------------------------------------------
sommet_gauche := polygone -> sort(polygone, key=(x->x[1]), output=permutation)[1]:



#-------------------------------------------------------------------
nouveau_polygone := proc(polygone, i_debut, i_fin)
  local n, p, i:

  n := numelems(polygone):
  p := NULL:
  i := i_debut:

  while i <> i_fin do
    p := p, polygone[i]:
    i := voisin_sommet(n, i, 1)
  end do:
  p := [p, polygone[i_fin]]
end proc:



#-------------------------------------------------------------------
trianguler_polygone_recursif := proc(polygone)
  local n, j0, j1, j2, P0, P1, P2, j, polygone_1, polygone_2:
  global liste_triangles:
  n  := numelems(polygone):
  j0 := sommet_gauche(polygone):
  j1 := voisin_sommet(n, j0, +1):
  j2 := voisin_sommet(n, j0, -1):
  P0 := polygone[j0]:
  P1 := polygone[j1]:
  P2 := polygone[j2]:
  j  := sommet_distance_maximale(polygone, P0, P1, P2, [j0, j1, j2]):

  if `not`(j::posint) then
    liste_triangles := liste_triangles, [P0, P1, P2]:
    polygone_1      := nouveau_polygone(polygone,j1,j2):
    if numelems(polygone_1) = 3 then
      liste_triangles := liste_triangles, polygone_1:
    else
      thisproc(polygone_1)
    end if:

  else
    polygone_1 := nouveau_polygone(polygone, j0, j ):
    polygone_2 := nouveau_polygone(polygone, j , j0):
    if numelems(polygone_1) = 3 then
      liste_triangles := liste_triangles, polygone_1:
    else
      thisproc(polygone_1)
    end if:
    if numelems(polygone_2) = 3 then
      liste_triangles := liste_triangles, polygone_2:
    else
      thisproc(polygone_2)
    end if:
  end if:

  return [liste_triangles]:
end proc:


#-------------------------------------------------------------------
trianguler_polygone := proc(polygone)
  trianguler_polygone_recursif(polygone):
  return liste_triangles:
end proc:


#-------------------------------------------------------------------
hachurer := proc(shapes, hatch_angle::list, hatch_number::list, hatch_color::list)

local A, La, Lp;
local N, P, _sides, L_sides, Xshape, ch, rel, p_rel, n, sol, p_range:
local AllHatches, window, p, _segment:
local NT, ka, N_Hatches, p_range_t, nt, shape, p_hatches, WhatHatches:

#-----------------------------------------------------------------
# internal functions:
#
# La(x, y, alpha, p) is the implicit equation of a straight line of angle alpha relatively
#                    to the horizontal axis and intercept p
#
# Lp(x, y, P) is the implicit equation of a straight line passing through points P[1] and P[2]
#
# interieur_triangle(triangle, M)

La := (x, y, alpha, p) -> cos(alpha)*x - sin(alpha)*y + p;
Lp := proc(x, y, P::list) (x-P[1][1])*(P[2][2]-P[1][2]) - (y-P[1][2])*(P[2][1]-P[1][1] = 0) end proc;


p_range    := [+infinity, -infinity]:
NT         := numelems(shapes):

AllHatches := NULL:

for ka from 1 to numelems(hatch_angle) do
  A         := hatch_angle[ka]:
  N_Hatches := hatch_number[ka]:
  p_range_t := NULL:
  _sides    := []:
  L_sides   := []:
  rel       := []:
  for nt from 1 to NT do

    shape := shapes[nt]:
    # _sides  : two points description of the sides of the shape
    # L_sides : implicit equations of the straight lines that support the sides

    N        := [1, 2, 3];
    P        := [2, 3, 1];
    _sides   := [ _sides[] , [ seq([shape[n], shape[P[n]]], n in N) ] ];
    L_sides  := [ L_sides[], Lp~(x, y, _sides[-1]) ];

    # Inequalities that define the interior of the shape

    rel := [ rel[], trianguler:-interieur_forme(shape, [x, y]) ];
  
    # Given the orientation of the hatches we search here the extreme values of
    # the intercept p for wich a straight line of equation La(x, y, alpha, p)
    # cuts the shape.
    
    p_rel := NULL:
    
    for n from 1 to numelems(L_sides[-1]) do
      sol   := solve({La(x, y, A, q), lhs(L_sides[-1][n])} union rel[-1], [x, y]);
      p_rel := p_rel, `if`(sol <> [], [rhs(op(1, %)), rhs(op(3, %))], [+infinity, -infinity]);
    end do:
    p_range_t := p_range_t, evalf(min(op~(1, [p_rel]))..max(op~(2, [p_rel])));
    p_range   := evalf(min(op~(1, [p_rel]), op(1, p_range))..max(op~(2, [p_rel]), op(2, p_range)));

  end do: # end of the loop over triangles

  p_range_t := [p_range_t]:
  p_hatches := [seq(p_range, (op(2, p_range)-op(1, p_range))/N_Hatches)]:
  # Building of the hatches
  #
  # This construction is far from being optimal.
  # Here again the main goal was to obtain the hatches with a minimum effort
  # if algorithmic development.

  window      := min(op~(1..shape))..max(op~(1..shape)):
  WhatHatches := map(v -> map(u -> if verify(u, v, 'interval'('closed') ) then u end if, p_hatches), p_range_t):

  for nt from 1 to NT do
    for p in WhatHatches[nt] do
      _segment := []:
      for n from 1 to numelems(L_sides[nt]) do
         _segment := _segment, evalf( solve({La(x, y, A, p), lhs(L_sides[nt][n])} union rel[nt], [x, y]) );
      end do;
      map(u -> u[], [_segment]);
      AllHatches := AllHatches, plot(map(u -> rhs~(u), %), color=hatch_color[ka]):
    end do:
  end do;

end do: # end of loop over hatch angles

plots:-display(
  PLOT(POLYGONS(polygone, COLOR(RGB, 1$3))),
  AllHatches,
  scaling=constrained
)

end proc:

end module:

 

Legrand's example (see reference above)

 

 

global liste_triangles:
liste_triangles := NULL:

polygone := [[0,0],[0.5,-1],[1.5,-0.2],[2,-0.5],[2,0],[1.5,1],[0.3,0],[0.5,1]]:

trianguler:-trianguler_polygone(polygone):

PLOT(seq(POLYGONS(u, COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12)), u in liste_triangles), VIEW(0..2, -2..2))

 

trianguler:-hachurer([liste_triangles], [-Pi/4, Pi/4], [40, 40], [red, blue])
 

 

F := (P, a, b) -> map(p -> [p[1]+a, p[2]+b], P):

MOTIF  := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]];
motifs := [ seq(seq(F(MOTIF, 0+i*0.075, 0+j*0.075), i=0..26), j=-14..13) ]:

plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):
  
for ref in liste_triangles do
  ref;
  #
  # the three lines below are used to define REF counter clockwise
  #
  g           := trianguler:-sommet_gauche(ref):
  bas         := sort(op~(2, ref), output=permutation);
  REF         := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(u -> plots:-inequal(rel_ref union u, x=0..2, y=-1..1, color=blue, 'nolines'), rel_motifs):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

[[0, 0], [0.5e-1, 0], [0.5e-1, 0.5e-1], [0, 0.5e-1]]

 

 

MOTIF  := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]];
motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.05, 0.1)+i*0.1, 0+j*0.05), i=-0.2..20), j=-20..20) ]:
plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):
  
for ref in liste_triangles do
  ref;
  g := trianguler:-sommet_gauche(ref):
  bas := sort(op~(2, ref), output=permutation);
  REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(u -> plots:-inequal(rel_ref union u, x=0..2, y=-1..1, color=blue, 'nolines'), rel_motifs):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

[[0, 0], [0.5e-1, 0], [0.5e-1, 0.5e-1], [0, 0.5e-1]]

 

 

MOTIF  := [[0, 0], [0.4, 0], [0.4, 0.14], [0, 0.14]]:
motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.4, 0.2)+i*0.4, 0+j*0.14), i=-1..4), j=-8..7) ]:


plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

palettes := ColorTools:-PaletteNames():
ColorTools:-GetPalette("HTML"):

couleurs := [SandyBrown, PeachPuff, Peru, Linen, Bisque, Burlywood, Tan, AntiqueWhite,      NavajoWhite, BlanchedAlmond, PapayaWhip, Moccasin, Wheat]:

nc   := numelems(couleurs):
roll := rand(1..nc):

motifs_nb      := numelems(motifs):
motifs_couleur := [ seq(cat("HTML ", couleurs[roll()]), n=1..motifs_nb) ]:

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):
  
for ref in liste_triangles do
  ref;
  g := trianguler:-sommet_gauche(ref):
  bas := sort(op~(2, ref), output=permutation);
  REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(n -> plots:-inequal(rel_ref union rel_motifs[n], x=0..2, y=-1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

 

MOTIF  := [ seq(0.1*~[cos(Pi/6+Pi/3*i), sin(Pi/6+Pi/3*i)], i=0..5) ]:
motifs := [ seq(seq(F(MOTIF, i*0.2*cos(Pi/6)+piecewise(j::odd, 0, 0.08), j*0.3*sin(Pi/6)), i=0..12), j=-6..6) ]:


plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):


motifs_nb      := numelems(motifs):
motifs_couleur := [ seq(`if`(n::even, yellow, brown) , n=1..motifs_nb) ]:

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):
  
for ref in liste_triangles do
  ref;
  g := trianguler:-sommet_gauche(ref):
  bas := sort(op~(2, ref), output=permutation);
  REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(n -> plots:-inequal(rel_ref union rel_motifs[n], x=0..2, y=-1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

 

 


 

Download Triangulation_Hatching_Texturing.mw

We recently had a question about using some of the plotting commands in Maple to draw things. We were feeling creative and thought why not take it a step further and draw something in 3D.

Using the geom3d, plottools, and plots packages we decided to make a gingerbread house.

To make the base of the house we decided to use 2 cubes, as these would give us additional lines and segments for the icing on the house.

point(p__1,[2,3,2]):
point(p__2,[3,3,2]):
cube(c1,p__1,2):
cube(c2,p__2,2):
base:=draw([c1,c2],color=tan);

Using the same cubes but changing the style to be wireframe and point we made some icing lines and decorations for the gingerbread house.

base_decor1:=draw([c1,c2],style=wireframe,thickness=3,color=red,transparency=0.2):
base_decor2:=draw([c1,c2],style=wireframe,thickness=10,color=green,linestyle=dot):
base_decor3:=draw([c1,c2],style=point,thickness=2,color="Silver",symbol=sphere):
base_decor:=display(base_decor1,base_decor2,base_decor3);

To create the roof we found the vertices of the cubes and used those to find the top corners of the base.

v1:=vertices(c1):
v2:=vertices(c2):
pc1:=seq(point(pc1||i,v1[i]),i=1..nops(v1)):
pc2:=seq(point(pc2||i,v2[i]),i=1..nops(v2)):
topCorners:=[pc1[5],pc1[6],pc2[1],pc2[2]]:
d1:=draw(topCorners):

Using these top corners we found the midpoints (where the peak of the roof would be) and added the roof height to the coordinates.

midpoint(lc1,topCorners[1],topCorners[2]):
detail(lc1);

point(cc1,[-(2*sqrt(3))/3 + 2, (2*sqrt(3))/3 + 3+1, 2]):
d3:=draw(cc1):

midpoint(lc2,topCorners[3],topCorners[4]):
detail(lc2);

point(cc2,[(2*sqrt(3))/3 + 3, (2*sqrt(3))/3 + 3+1, 2]):
d4:=draw(cc2):

With the midpoints and vertices at the front and rear of the house we made two triangles for the attic of the gingerbread house.

triangle(tf,[topCorners[1],topCorners[2],cc1]):
front:=draw(tf,color=brown):

triangle(tb,[topCorners[3],topCorners[4],cc2]):
back:=draw(tb,color=tan):

Using these same points again we made more triangles to be the roof.

triangle(trl,[cc1,cc2,pc1[5]]):
triangle(trh,[pc2[2],pc1[6],cc1]):
triangle(tll,[cc1,cc2,pc2[2]]):
triangle(tlh,[pc2[1],pc1[5],cc2]):
roof:=draw([trl,trh,tll,tlh],color="Chocolate");

Our gingerbread house now had four walls, a roof, and icing, but no door. Creating the door was as easy as making a parallelepiped, but what is a door without more icing?

door:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="DarkRed"):
door_decor1:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="Gold",style=line):
door_decor2:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="Silver", style=line,linestyle=dot,thickness=5):
door_decor:=display(door_decor1,door_decor2):

Now having a door we could have left it like this, but what better way to decorate a gingerbread house than with candy canes? Naturally, if we were going to have one candy cane we were going to have lots of candy canes. To facilitate this we made a candy cane procedure.

candy_pole:=proc(c:=[0,0,0], {segR:=0.1}, {segH:=0.1}, {segn:=7}, {tilt_theta:=0}, {theta:=0}, {arch:=true}, {flip:=false})
local cane1,cane2,cane_s,cane_c,cane0,cane,i,cl,cd,ch, cane_a,tmp,cane_ac,cane_a1,cane00,cane01,cane02,cane_a1s,tmp2,cane_a2s:
uses plots,geom3d:
cl:=c[1]:
cd:=c[2]:
ch:=c[3]:
cane1:=plottools:-cylinder([cd, ch, cl], segR, segH,style=surface):

cane2:=display(plottools:-rotate(cane1,Pi/2,[[cd,ch,cl],[cd+1,ch,cl]])):
cane_s:=[cane2,seq(display(plottools:-translate(cane2,0,segH*i,0)),i=1..segn-1)]:
cane_c:=seq(ifelse(type(i,odd),red,white),i=1..segn):

cane0:=display(cane_s,color=[cane_c]):

if arch then

cane_a:=plottools:-translate(cane2,0,segH*segn-segH/2,0):
tmp:=i->plottools:-rotate(cane_a,i*Pi/24, [ [cd,ch+segH*segn-segH/2,cl+segR*2] , [cd+1,ch+segH*segn-segH/2,cl+segR*2] ] ):

cane_ac:=seq(ifelse(type(i,odd),red,white),i=1..24):

                cane_a1s:=seq(plottools:-translate(tmp(i),0,segH*i/12,segR*i/4),i=1..12):

tmp2:=i->plottools:-rotate(cane_a1s[12],i*Pi/24,[[cd,ch+segH*segn-0.05,cl+segR*2],[cd+1,ch+segH*segn-0.05,cl+segR*2]]):

cane_a2s:=seq(plottools:-translate(tmp2(i),0,-segH*i/500,segR*i/4),i=1..12):
cane_a1:=display(cane_a1s,cane_a2s,color=[cane_ac]):
cane00:=display(cane0,cane_a1);

                if flip then

cane01:=plottools:-rotate(cane00,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):
cane02:=plottools:-rotate(cane01,theta,[[cd,ch,cl],[cd,ch+1,cl]]):
cane:=plottools:-reflect(cane01,[[cd,ch,cl],[cd,ch+1,cl]]):

                else

cane01:=plottools:-rotate(cane00,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):
cane:=plottools:-rotate(cane01,theta,[[cd,ch,cl],[cd,ch+1,cl]]):

                end if:

                return cane:

else

                cane:=plottools:-rotate(cane0,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):

                return cane:

end if:

end proc:

With this procedure we decided to add candy canes to the front, back, and sides of the gingerbread house. In addition we added two candy poles.

Candy Canes in front of the house:

cane1:=candy_pole([1.2,0,2],segn=9,arch=false):
cane2:=candy_pole([2.8,0,2],segn=9,arch=false):
cane3:=candy_pole([2.7,0.8,3.3],segn=9,segR=0.05,tilt_theta=-Pi/7):
cane4:=candy_pole([1.3,0.8,3.3],segn=9,segR=0.05,tilt_theta=-Pi/7,flip=true):
front_canes:=display(cane1,cane2,cane3,cane4):

Candy Canes at the back of the house:

caneb3:=candy_pole([1.5,4.2,2.5],segn=15,segR=0.05,tilt_theta=-Pi/3,flip=true):
caneb4:=candy_pole([2.5,4.2,2.5],segn=15,segR=0.05,tilt_theta=-Pi/3):}
back_canes:=display(caneb3,caneb4):

Candy Canes at the left of the house:

canel1:=candy_pole([0.8,1.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
canel2:=candy_pole([0.8,2.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
canel3:=candy_pole([0.8,4,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
left_canes:=display(canel1,canel2,canel3):

Candy Canes at the right of the house:

caner1:=candy_pole([3.2,1.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
caner2:=candy_pole([3.2,2.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
caner3:=candy_pole([3.2,4,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
right_canes:=display(caner1,caner2,caner3):

canes:=display(front_canes,back_canes,right_canes,left_canes):

With these canes in place all that was left was to create the ground and display our Gingerbread House.

ground:=display(plottools:-parallelepiped([5,0,0],[0,0.5,0],[0,0,4],[0,1.35,0]),color="DarkGreen"):

display([door,door_decor,d1,base,base_decor,d3,d4,front,back,roof,ground,canes],orientation=[-100,0,95]);

You can download the full worksheet creating the gingerbread house below:

Geometry_Gingerbread.mw

There was a discussion regarding Maple and two light sources a couple of days ago (unfortunately the content of that conversation has been removed by the original poster so if any of you thought you were getting early signs of alzheimers, rest assured the questions/posts were indeed there and you can all put your memories at ease for a little).

I did indeed go back into Maple to see how far it was that we lost the two light source capability.  I can tell you that as far back as Maple 10.06 dual (or more) light sources was not working.  So I worked my way up - of course the OP of the deleted questions mentioned Maple Vr4 (I think), so I started there.  Maple 6 and Maple 7 is ok.  I was unable to check Maple 8 and 9 as my computers with that software are currently missing (probably buried in my shed somewhere - as long as the wife hasn't got rid of them.. anyways)

So I'll share with you what once was.  From Maple 7, and an animated gif (which hopefully works) to illustrate two light sources (one brown light, and one mauve or purple light) and what it looked like from Maple 7.

with(plots):
setoptions3d(style=patchnogrid,projection=.9):
a:=plot3d(-x^2-y^2-x*y,x=-1..1,y=-1..1,color=[.5,.9,.9],grid=[15,15]):
b:=n->display(a,orientation=[n,60]):
display([seq(b(10*n),n=0..35)],light=[90,-80,1,.2,.1],light=[90,40,1,.5,.8],insequence=true);

The animation is not working so the uploaded gif file is here for you to check out. .. Couldn't upload the gif file so I've zipped it - that worked

twolightsource.zip

Program of Transfer matrix method for solving the vibrational behavior of a cracked beam. For more details and comprehension check the paper entitled ''A transfer matrix method for free vibration analysis and crack identification of stepped beams with multiple edge cracks and different boundary conditions ''
 

restart; with(LinearAlgebra)NULLNULL

W11 := A[1]*cos(nu*x)+A[2]*sin(nu*x)+A[3]*cosh(nu*x)+A[4]*sinh(nu*x);

theta11 := -A[1]*nu*sin(nu*x)+A[2]*nu*cos(nu*x)+A[3]*nu*sinh(nu*x)+A[4]*nu*cosh(nu*x)  NULL

M11 := EI*(-A[1]*nu^2*cos(nu*x)-A[2]*nu^2*sin(nu*x)+A[3]*nu^2*cosh(nu*x)+A[4]*nu^2*sinh(nu*x))
S11 := EI*(A[1]*nu^3*sin(nu*x)-A[2]*nu^3*cos(nu*x)+A[3]*nu^3*sinh(nu*x)+A[4]*nu^3*cosh(nu*x))
 

MD11 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, W11); MD12 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, W11); MD13 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, W11); MD14 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, W11)
NULL

MD21 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, theta11); MD22 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, theta11); MD23 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, theta11); MD24 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, theta11)
 

MD31 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, M11); MD32 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, M11); MD33 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, M11); MD34 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, M11)
 

MD41 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, S11); MD42 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, S11); MD43 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, S11); MD44 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, S11)
 

 

TM := Matrix(4, 4, [[MD11, MD12, MD13, MD14], [MD21, MD22, MD23, MD24], [MD31, MD32, MD33, MD34], [MD41, MD42, MD43, MD44]])

C := Matrix(4, 4, [[1, 0, 0, 0], [0, 1, c44, 0], [0, 0, 1, 0], [0, 0, 0, 1]])NULL

NULL    TM2 := subs(x = 0, TM); TM3 := subs(x = L, TM)

with(MTM)

TM4 := inv(TM)NULLNULL 

TM5 := inv(TM2)NULL

Y11 := MatrixMatrixMultiply(TM3, TM4)

    Y22 := MatrixMatrixMultiply(C, TM)   

Y33 := MatrixMatrixMultiply(Y11, Y22)

Y44 := MatrixMatrixMultiply(Y33, TM5)NULL

BB11 := Y44[3, 3]

BB12 := Y44[3, 4]

BB21 := Y44[4, 3]

BB22 := Y44[4, 4] 

BB := Matrix(2, 2, [[BB11, BB12], [BB21, BB22]]) 

NULL

R11 := det(BB) 

NULL

L := .18; L1 := 1; h := 0.5e-2; b := 0.2e-1; rho := 957.5; area = b.h; m := rho*h*b; EI := 0.2682e10*b*h^3mu := ((m.(omega^2))*L^4/EI)^(1/4); x := .5; c44 := .1; c11 := 0NULLNULL

plot(R11, omega = 1 .. 100)

 

 

 

``

NULL

NULL


 

Download transfer.mw

In the applications I am working on, the information are often represented by hierarchical tables (that is tables where some entries can also be tables, and so on).
To help people to understand how this information is organized, I have thought to representent this hierarchical table as a tree graph.
Once this graph built, it becomes very simple to find where a "terminal leaf", that is en entry which is no longer a table, is located in the original table (by location I mean the sequence of indices for which the entry is this "terminal leaf".

The code provided here is pretension free and I do not doubt a single second  that people here will be able to improve it.
I published it for i thought other people could face the same kind of problems that I do.


 

restart

with(GraphTheory):
interface(version);

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

gh := proc(T)
  global s, counter, types:
  local  i:
  if type(T, table) then
    for i in [indices(T, nolist)] do
      if type(T[i], table) then
         s := s, op(map(u -> [i, u], [indices(T[i], nolist)] ));
      else
         counter := counter+1:
         types   := types, _Z_||counter = whattype(T[i]);
         s       := s, [i, _Z_||counter];
      end if:
      thisproc(T[i]):
    end do:
  else
    return s
  end if:
end proc:

t := table([a1=[alpha=1, beta=2], a2=table([a21=2, a22=table([a221=x, a222=table([a2221={1, 2, 3}, a2222=Matrix(2, 2), a2223=u3, a2224=u4])])]), a3=table([a31=u, a32=v])]);

global s, counter, types:
s       := NULL:
counter := 0:
types   := NULL:

ghres := gh(t):
types := [types]:

t := table([a1 = [alpha = 1, beta = 2], a3 = table([a32 = v, a31 = u]), a2 = table([a22 = table([a222 = table([a2222 = (Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0})), a2223 = u3, a2221 = {1, 2, 3}, a2224 = u4]), a221 = x]), a21 = 2])])

(2)


These 3 lines determine the set of edges of the form ['t', v], that are not been captured by procedure h.
They correspond to "first level" indices of table t (v in {a1, a2, a3} in the example above)

L := convert(op~(1, [ghres]), set):     
R := convert(op~(2, [ghres]), set):
FirstLevelEdges := map(u -> ['t', u], L union R minus R):


Complete the set of the edges, build the graph representation TG of table t and draw TG.

edges := convert~({ghres, FirstLevelEdges[]}, set):
TG := Graph(edges):

HighlightVertex(TG, Vertices(TG), white):
p := DrawGraph(TG, style=tree, root='t'):
 


The first line is used to change the the "terminal leaves" of names  _Z_n by their type.

eval(t);

p       := subs(types, p):
enlarge := plottools:-transform((x,y) -> [3*x, y]):

plots:-display(enlarge(p), size=[1000, 400]);

table([a1 = [alpha = 1, beta = 2], a3 = table([a32 = v, a31 = u]), a2 = table([a22 = table([a222 = table([a2222 = (Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0})), a2223 = u3, a2221 = {1, 2, 3}, a2224 = u4]), a221 = x]), a21 = 2])])

 

 


This procedure is used to find the "indices path" to a terminal leaf.
FindLeaf is then applied to all the terminal leaves.

FindLeaf := proc(TG, leaf)
   local here:
   here := GraphTheory:-ShortestPath(TG, 't', leaf)[1..-2]:
   here := cat(convert(here[1], string), convert(here[2..-1], string)):
   here := StringTools:-SubstituteAll(here, ",", "]["):
   here := parse(here);
end proc:

# where is a2221

printf("%a\n", FindLeaf(TG, a2221));

t[a2][a22][a222]

 

 


 

Download Table_Unfolding_2.mw

 

First 19 20 21 22 23 24 25 Last Page 21 of 76