MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • Maple 2020 offers many improvements motivated and driven by our users.

    Every single update in a new release has a story behind it. It might be a new function that a customer wants, a response to some feedback about usability, or an itch that a developer needs to scratch.

    I’ll end this post with a story about acoustic guitars and how they drove improvements in signal and audio processing. But first, here are some of my personal favorites from Maple 2020.

    Graph theory is a big focus of Maple 2020. The new features include more control over visualization, additional special graphs, new analysis functions, and even an interactive layout tool.

    I’m particularly enamoured by these:

    • We’ve introduced new centrality measures - these help you determine the most influential vertices, based on their connections to other vertices
    • You now have more control over the styling of graphs – for example, you can vary the size or color of a nodebased on its centrality

    I’ve used these two new features to identify the most influential MaplePrimes users. Get the worksheet here.

    @Carl Love – looks like you’re the biggest mover and shaker on MaplePrimes (well, according to the eigenvector centrality of the MaplePrimes interaction graph).

    We’ve also started using graph theory elsewhere in Maple. For example, you can generate static call graph to visualize dependencies between procedures calls in a procedure

    You now get smoother edges for 3d surfaces with non-numeric values. Just look at the difference between Maple 2019 and 2020 for this plot.

    Printing and PDF export has gotten a whole lot better.  We’ve put a lot of work into the proper handling of plots, tables, and interactive components, so the results look better than before.

    For example, plots now maintain their aspect ratio when printed. So your carefully constructed psychrometric chart will not be squashed and stretched when exported to a PDF.

    We’ve overhauled the start page to give it a cleaner, less cluttered look – this is much more digestible for new users (experienced users might find the new look attractive as well!). There’s a link to the Maple Portal, and an updated Maple Fundamentals guide that helps new users learn the product.

    We’ve also linked to a guide that helps you choose between Document and Worksheet, and a link to a new movie.

    New messages also guide new users away from some very common mistakes. For example, students often type “e” when referring to the exponential constant – a warning now appears if that is detected

    We’re always tweaking existing functions to make them faster. For example, you can now compute the natural logarithm of large integers much more quickly and with less memory.

    This calculation is about 50 times faster in Maple 2020 than in prior versions:

    Many of our educators have asked for this – the linear algebra tutorials now return step by step solutions to the main document, so you have a record of what you did after the tutor is closed.

    Continuing with this theme, the Student:-LinearAlgebra context menu features several new linear algebra visualizations to the Student:-LinearAlgebra Context Menu. This, for example, is an eigenvector plot.

    Maple can now numerically evaluate various integral transforms.

    The numerical inversion of integral transforms has application in many branches of science and engineering.

    Maple is the world’s best tool for the symbolic solution of ODEs and PDEs, and in each release we push the boundary back further.

    For example, Maple 2020 has improved tools for find hypergeometric solutions for linear PDEs.

    This might seem like a minor improvement that’s barely worth mentions, but it’s one I now use all the time! You can now reorder worksheet tabs just by clicking and dragging.

    The Hough transform lets you detect straight lines and line segments in images.

    Hough transforms are widely used in automatic lane detection systems for autonomous driving. You can even detect the straight lines on a Sudoku grid!

    The Physics package is always a pleasure to write about because it's something we do far better than the competition.

    The new explore option in TensorArray combines two themes in Maple - Physics and interactive components. It's an intuitive solution to the real problem of viewing the contents of higher dimensional tensorial expressions.

    There are many more updates to Physics in Maple 2020, including a completely rewritten FeynmanDiagrams command.

    The Quantum Chemistry Toolbox has been updated with more analysis tools and curriculum material.

    There’s more teaching content for general chemistry.

    Among the many new analysis functions, you can now visualize transition orbitals.

    I promised you a story about acoustic guitars and Maple 2020, didn’t I?

    I often start a perfectly innocuous conversation about Maple that descends into several weeks of intense, feverish work.

    The work is partly for me, but mostly for my colleagues. They don’t like me for that.

    That conversation usually happens on a Friday afternoon, when we’re least prepared for it. On the plus side, this often means a user has planted a germ of an idea for a new feature or improvement, and we just have to will it into existence.

    One Friday afternoon last year, I was speaking to a user about acoustic guitars. He wanted to synthetically generate guitar chords with reverb, and export the sound to a 32-bit Wave file. All of this, in Maple.

    This started a chain of events that that involved least-square filters, frequency response curves, convolution, Karplus-Strong string synthesis and more. We’ll package up the results of this work, and hand it over to you – our users – over the next one or two releases.

    Let me tell you what made it into Maple 2020.

    Start by listening to this:

    It’s a guitar chord played twice, the second time with reverb, both generated with Maple.

    The reverb was simulated with convolving the artificially generated guitar chord with an impulse response. I had a choice of convolution functions in the SignalProcessing and AudioTools packages.

    Both gave the same results, but we found that SignalProcessing:-Convolution was much faster than its AudioTools counterpart.

    There’s no reason for the speed difference, so R&D modified AudioTools:-Convolution to leverage SignalProcessing:-Convolution for the instances for which their options are compatible. In this application, AudioTools:-Convolution is 25 times faster in Maple 2020 than Maple 2019!

    We also discovered that the underlying library we use for the SignalProcessing package (the Intel IPP) gives two options for convolution that we were previously not using; a method which use an explicit formula and a “fast” method that uses FFTs. We modified SignalProcessing:-Convolution to accept both options (previously, we used just one of the methods),

    That’s the story behind two new features in Maple 2020. Look at the entirety of what’s new in this release – there’s a tale for each new feature. I’d love to tell you more, but I’d run out of ink before I finish.

    To read about everything that’s new in Maple 2020, go to the new features page.

    Today we celebrated International Women's Day at Maplesoft. As part of our celebration, we had a panel of 5 successful women from within the community share their experiences and insights with us. 

    Hearing these women speak has given me the courage to share my personal experience and advice to women in technology. If what I write here helps even one woman, then I will have accomplished something great today. 


    What do you do at Maplesoft?

    My name is Karishma. I'm the Director, Product Management - Academic. 


    Where did you grow up and where did you go to school (Diploma/degree)?

    I was born and raised in Montreal to parents of Indian descent. Like most Indian parents, they “encouraged" me to pursue a career in either Law, Medicine, or Engineering, despite my true calling to pursue a career in theatre (at least that's what I believed it to be at the time)

    Given that I had no siblings to break the ice, and that rebelling wasn't my Modus Operandi (that came much later), I did what any obedient teenager would do: I pursued a career in Electrical Engineering at McGill University. In my mind, this was the fastest way to landing a job and fleeing the proverbial nest. 

    Electrical Engineering was far from glamorous, and after two years, I was ready to switch. It was due to the sheer insistence of my mother that I completed the degree. 

    So how did I end up pursuing a graduate degree in Biomedical Engineering at McGill University? It wasn't the future I envisioned, but the economic downturn in 2001-2002 saw a massive decrease in hiring, and the job that I had held-out patiently for during those four years became a far-off dream. So I did the thing I never imagined I would: I accepted the offer to pursue a Master's and the very generous stipend that came with it. In case you are wondering, I only applied because my father nagged me into submission. (Insistence and nagging are two innate traits of Indian parents)

    Contrary to what I expected, I loved my Master's degree! It gave me the freedom to immerse myself wholly in a topic I found exciting and allowed me to call the shots on my schedule, which led to my involvement in student government as VP Internal. But apart from the research and the independence, pursuing a master's degree opened doors to opportunities that I couldn't have imagined, such as an internship with the International Organization for Migration in Kenya, a job offer in Europe, and the chance to work at Maplesoft. (I guess my parents did know what was best for me.) 


    What is the best part of your job?

    It's figuring out how to solve problems our users have as well as the ones they might not realize they have. 

    At Maplesoft, I work with some most brilliant minds I've ever encountered to build a product that makes math more accessible to our users, whether they be a student, researcher, scientist, or engineer. 

    Some of the aspects of my role that I love the most include: 

    • speaking to and learning from our customers, 
    • interpreting the meaning behind their words, facial expressions, vocal intonations, and body language, and
    • collaborating with the sales, marketing, and development teams to turn what was 'said' into tangible actions that will enhance the product and user experience. 

    Most nights, when I leave work, I do so with a sense of excitement because I know my actions and the actions of those I work with will help our users achieve their goals and ambitions. There's no better high. 


    What advice do you have for young women interested in a career in your field? 
    Throughout my career, I've had the privilege to work with some amazing women and men who've given me advice that I wish I had known when I was an undergraduate student. If you are a woman pursuing a STEM degree or starting your first job in a tech firm, here are three tips that may help you: 

    1.   Don't be scared of the 'N' word. 
    Don't be scared of NETWORKING. I know it can be intimidating, but it truly is the best way to land a job, advance your career, or meet the person you admire most. Remember that networking can take place anywhere - it's not exclusive to networking events. Some advice that I received that helped me overcome my fear of networking: 

    • Smile - Before you approach a person or enter a networking session, force yourself to smile. It will help you diffuse any tension you are holding and will make you appear more approachable. 
    • Research - Take the time to research the person(s) you would like to meet. Find out as much as you can about them and their company. Prepare some icebreaker questions and other questions to help carry the conversation forward ahead of time. Remember that people like to talk about themselves and their experiences. 
    • Don't take it personally - The person you approach may find networking equally tricky. So if they seem disinterested or aloof, don't take it personally. 
    • Just do it - Networking gets easier with practice. Don't let a failed attempt set you back. The worse thing that will happen is that you don't make a connection. 


     2.   It's ok to ask for help.
    If you are a woman in an environment that is dominated by men, you might hesitate to ask for help. DON'T! There's nothing wrong with asking for help. That said, many women ask for help in a way that undermines their confidence and thus erodes others’ perception of them. Next time you need help, have a question or require clarification, take a moment to phrase your request, so you don't inadvertently put yourself down. 


    3.   Play to your strengths
    Don't think you need to know everything. Nobody expects it. If you landed a new job or co-op placement, and you are finding yourself doing things you've never done and don't come naturally to you yet, don't let your brain convince you that you don't deserve it. Remember that you earned it because of your qualities and strengths. 

    When discussing Maple programming, we often refer to for-loops, while-loops, until-loops, and do-loops (the latter being an infinite loop). But under the hood, Maple has only two kinds of loop, albeit very flexible and powerful ones that can combine the capabilities of any or all of the above, making it possible to write very concise code in a natural way.

    Before looking at some actual examples, here is the formal definition of the loops' syntax, expressed in Wirth Syntax Notation, where "|" denotes alternatives, "[...]" denotes an optional part, "(...)" denotes grouping, and Maple keywords are in boldface:

    [ for  ] [ from  ] [ by  ] [ to  ]
        [ while  ]
    ( end do | until  )
    [ for  [ , variable ] ] in 
        [ while  ]
    ( end do | until  )

    In the first form, every part of the loop syntax is optional, except the do keyword before the body of the loop, and either end do or an until clause after the body. (For those who prefer it, end do can also be written as od.) In the second form, only the in clause is required.

    The simplest loop is just:

    end do

    This will repeat the forever, unless a break or return statement is executed, or an error occurs.

    One or two loop termination conditions can be added:

    • A while clause can be written before the do, specifying a condition that is tested before each iteration begins. If the condition evaluates to false, the loop ends.
    • An until clause can be written instead of the end do, specifying a condition that is tested after each iteration finishes. If the condition evaluates to true, the loop ends.

    A so-called for-loop is just a loop to which iteration clauses have been added. These can take one of two forms:

    • Any combination of for (with a single variable), from, by, and to clauses. The last three can appear in any order.
    • A for clause with one or two variables, followed by an in clause.

    The following for-loop executes 10 times:

    for  from 1 to 10 do
    end do

    However, if the doesn't depend on the value of , both the for and from clauses can be omitted:

    to 10 do
    end do

    In this case, Maple supplies an implicit for clause (with an inaccessible internal variable), as well as an implicit "from 1" clause. In fact, all of the clauses are optional, and the infinite loop shown earlier is understood by Maple in exactly the same way as:

    for  from 1 by 1 to infinity while true do
    until false

    When looping over the contents of a container, such as a one-dimensional array A, there are several possible approaches. The one closest to how it would be done in most other programming languages is (this example and those that follow can be copied and pasted into a Maple session):

     := Array([,"foo",42]);
    for  from lowerbound() to upperbound() do
    end do;

    If only the entries in the container are of interest, it is not necessary to loop over the indices. Instead, one can write:

     := Array([,"foo",42]);
    for  in  do
    end do;

    If both the indices and values are needed, one can write:

     := Array([,"foo",42]);
    for ,  in  do
    end do;

    For a numerically indexed container such as an Array, this is equivalent to the for-from-to example. However, this method also works with arbitrarily indexed containers such as a Matrix or table:

     := LinearAlgebra:-RandomMatrix(2,3);
    for ,  in  do
    end do;
     := table({1="one","hello"="world",=42});
    for ,  in eval() do
    end do;

    (The second example requires the call to eval due to last-name evaluation of tables in Maple, a topic for another post.)

    As with a simple do-loop, a while and/or until clause can be added. For example, the following finds the first negative entry, if any, in a Matrix (traversing the Matrix in storage order):

     := LinearAlgebra:-RandomMatrix(2,3);
    for ,  in  do
        # nothing to do here
    until  < 0;
    if  < 0 then
    end if;

    Notice that the test, < 0, is written twice, since it is possible that the Matrix has no negative entry. Another way to write the same loop but only perform the test once is as follows:

     := LinearAlgebra:-RandomMatrix(2,3);
    for ,  in  do
        if  < 0 then
        end if;
    end do;

    Here, we perform the test within the loop, perform the desired processing on the found value (just printing in this case), and use a break statement to terminate the loop.

    Sometimes, it is useful to abort the current iteration of the loop and move on to the next one. The next statement does exactly that. The following loop prints all the indices but only the positive values in a Matrix:

     := LinearAlgebra:-RandomMatrix(2,3);
    for ,  in  do
        if  < 0 then
        end if;
    end do;

    (Note that a simple example like this would be better written by enclosing the printing of the value in an if-statement instead of using next. The latter is generally only used if the former is not possible.)

    Maple's loop statements are very flexible and powerful, making it possible to write loops with complex combinations of termination conditions in a concise yet readable way. The ability to use while and/or until in conjunction with for means that break statements are often unnecessary, further improving clarity.

    The binary search algorithm is used to obtain the index of a given number by dividing the search bound in half over iteration. If the value entered in the array a message pop up telling that ''value is not present in the array". Please see the code. 

    restart; with(ArrayTools); AA := Array(1 .. 10, [20, 2, 30, 4, 50, 7, 60, 8, 90, 100]); AA := sort(AA); KEYVALUE := 200; DUP_KEYVALUE := infinity; low := 1; high := NumElems(AA); while DUP_KEYVALUE <> KEYVALUE do mid := floor((low+high)*(1/2)); if AA[mid] = KEYVALUE then DUP_KEYVALUE := KEYVALUE; printf("%s\n %a\n", "the index is ", mid) elif AA[mid] < KEYVALUE then low := mid+1 elif AA[mid] > KEYVALUE then high := mid-1 end if; if low > high then printf("%s\n", "value not present in the array"); break end if end do




    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!




    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
    od: NULL

    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.



    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]+`...`



    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:


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


    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].


    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.




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


    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}]


    Let L be the interaction Lagrangian

    L := lambda*phi(X)^4

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


    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]])


    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]])


    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]])


    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]]))


    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]])


    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]]))


    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]])


    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]])


    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]])


    Consider the interaction Lagrangian of Quantum Electrodynamics (QED). To formulate this problem on the worksheet, start defining the vector field 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)}


    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]


    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]


    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]])


    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]])


    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]])


    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]])


    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}]


    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]])


    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]])


    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]


    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])


    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))


    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*`' `




    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)))


    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}], {}


    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`


    Setup(quantumoperators = {W, Z})

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


    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)}


    CompactDisplay((W, Z)(X))

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


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


    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])


    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])


    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])))


    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 .




    [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.



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

    The ideas here are to allow 3D plotting commands such as plot3d to handle a `size` option similarly to how 2D plotting commands do so, and for the plots:-display command to also handle it for 3D plots.

    The size denotes the dimensions of the inlined plotting window, and not the relative lengths of the three axes.

    I'd be interested in any new problems introduced with this, eg. export, etc.


    # Using ToInert/FromInert
    # This might go in an initialzation file.
      if __ver>=18.0 and __ver<=2019.2 then
        if :-has(:-op([5,2,2,2,1],__KK),:-_Inert_PARAM(__NN)) then
          :-print("3D size patch done");
          :-print("3D size patch not appropriate; possibly already done");
        end if;
        :-print(sprintf("3D size patch not appropriate for version %a"),__ver);
      end if;
      :-print("3D size patch failed");
    end try:

    "3D size patch done"


    P := plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, size=[150,150],
                font=[Times,5], labels=["","",""]):

    plots:-display(P, size=[300,300], font=[Times,10]);

    # inherited from the contourplot3d (the plot3d is unset).
      plots:-contourplot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1,
                           thickness=3, contours=20, size=[800,800]),
      plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, color="Gray",
             transparency=0.1, style=surface)

    # Some options should still act as 2D-plot-specific.
    try plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, legend="Q");
        print("Not OK");
    if StringTools:-FormatMessage(lastexception[2..-1])
       ="the legend option is not available for 3-D plots"
    then print("OK"); else print("Not OK"); error; end if; end try;




    If this works fine then it might be a candidate for inclusion in an initialization file, so that it's
    automatically available.

    I was trying to display a Physics[Vectors] vector name in a 3dplot with an up arrow
    on it. I found that this old 2008 trick still works in MAPLE 2018.





    # Using MAPLE 2018.2


    t:= textplot3d([-1.1,1.1,1,v_]):




    # I found this on an old 2008 post
    t:= textplot3d([-1.1,1.1,1,typeset(`#mover(mi(` || v ||  `),mo("→"))`)]):




    Application of MapleSim in Science and Engineering: a simulationbased approach

    In this research work I show the methods of embedded components together with modeling and simulation carried out with Maple and MapleSim for the main areas of science and engineering (mathematics, physics, civil, mechanical etc); These two latest scientific softwares belonging to the company Maplesoft. Designed to be generated and used by teachers of education, as well as by university teachers and engineers; the results are highly optimal since the times saved in calculations are invested in analyzes and interpretations; among other benefits; in this way we can use our applications in the cloud since web technology supports Maple code with procedural and component syntax.


    Lenin AC

    Ambassador of Maple