@mmcdara I see, using union was just an experiment.

It's very handy that the HTML entities begin with ampersand (&), so that they can be used in this way. I showed you `≅`, but as dharr mentioned you can used others, many available from the Palettes. Your choice.

A few more notes, just in case. I'm still not quite sure whether you might have wanted to use T, literally. If so, please see the attached.

note. The infix notation is not a consequence of being a binary operator, and that's true mathematically as well as in Maple. Being a binary operator brings extra conditions to a binary function, eg. same domain for each argument, and (depending on who you read) same codomain. But infix notation is not a requirement. (eg. exponentiation). It just happens that some common algebraic operations have infix notation for operators.

Download infix_note.mw

note: I don't see & neutral operators in use so much any more. Physics has its sophisticated own stuff. Once upon a time, evalm used &* commonly for denoting noncommutative matrix multiplication.

@mmcdara Here are a few notes (my interpretations, right or wrong) on trying to force the harder order of the integrations, with nested 1D methods.

Using Maple 2024.1.

Download integration_issue_ac3.mw

Changing the order of the bounds (ie. which is the "inner" integral) makes a huge difference here.  The other order took so long I interrupted it.

Download Integration_issue_ac2.mw

@Ditteon That first link is just to a youtube page. You don't need any Maplesoft login to watch it.

My Chrome browser lets me watch it without signing in to Google.

Is there a special reason for not using a bare piecewise?

Download charfcn_plot_pw.mw

So that we can compare and double check suggestions, what version of Maple are you using?

Also, how many correct decimal digits do you need in the answer?

@segfault You should provide the actual example.

You might look at some of the approaches in these responses to these old postings: 1, 2, 3

note: Considerable portions of those approaches can now be done more straightforwardly, ie. in modern Maple. For example, using the background option of the plot command, using convert(plot(...),Image) rather than exporting a plot to an external temp file, etc. When I first used an image here, in the year 2012, even the ImageTools:-Embed command did not yet exist, so I had to first insert a Label Component manually.

How are you currently exporting your plots, for use in LaTeX?

What specific steps (mouse+GUI, or programmatic commands) are you using? What target file format? What resolution?

@paulmcquad For your examples, you can also use mmcdara's straightforward Answer.

@Kitonum I gave it a vote, and marked it as a favorite.

@erik10 Yes, the integrand has a singularity where z^2/(p^2*v__lin(z)^2) - 1, and is complex-valued for z making that expression negative.

For your positive data, at z=p*v__lin(z) .

Seismology_acc_ev_ac_ev_ac.mw

@erik10 So what is happening in your new sheet is that you are passing x(r) and y(r) to the plot command, and those evaluate to expressions involving theta(r).

But your procedure assigned to theta is not prepared to be called with symbolic r (ie. not an actual number). If you were to issue just theta(r) as a separate statement then you'd see the same error. So, those function calls passed to plot are evaluating prematurely.

One way around that is to use the so-called operator-form calling sequnce of the plot command, eg,
    plot([x, y, r__0 .. R], ...options)

Another way is to instead construct the theta procedure so that it returns unevaluated when called with a symbolic argument such as just a name. You could do that using MakeProcedure (unapply) with its numeric option, as I showed before for your Z__spline, etc. Or you could write it all out, which I do in the attachment.

I also add the numeric integration option epsilon=1e-5 here, which relaxes the accuracy tolerance for the theta computation (the outer integration). That makes your nested numeric integration faster, as I mentioned earlier that I suspected. It's now fast enough that you can get smooth adaptive plotting without too long a wait.

Seismology_acc_ev_ac.mw

@erik10 The adaptive=false option instructs the plot command to not do adaptive sampling of the domain. It makes it use an evenly spaced number of values in the domain.

The numpoints=N option instructs the plot command to use N values in the domain, initially. If adaptive=false then that means N evenly spread values, and otherwise it means N values initially in whichever adaptive sampling scheme it's using.

In other words, if one wants only N data points, evenly spaced, then one can supply,
   adaptive=false, numpoints=N
This is a way to know exactly how many evaluations of the function will be done, so is handy when it might be time consuming. When going this route I experiment to find the smallest number of values which make the curve appear smooth "enough". (Sometimes I start with N=3, as sanity-check experiment that things are working at all...)

Of course, an alternative is to use seq or Vector to get a specific number of evenly spaced data points. But I prefer to use the above options since if everything goes well or if I speed it up then I can just remove those options.

@erik10 You could use the MakeFunction (unapply) approach that I'd shown above:

Seismology_ac.mw

Or you could use an operator-form for your new integrand, as variantion on what I'd shown before for that kind of idea. (This variant now gains a slight inefficiency of a separate procedure construction for each call to theta, though in practice that can be relatively reasonable if the numeric quadrature is itself expensive enough.)

Seismology_acc.mw

[edit] I suspect that the final plot (of theta) might be computed faster if the epsilon tolerance were loosened on one of the (nested) numeric integrations. But for that you'd need to give up the pretty 2D integral look. I left it with non-adaptive plotting and a prettier integral, in these two worksheet attachments, and not too much of a wait.

ps. Often it's more useful to supply the background example up front.

pps. It's a good question, though, and has come up before. The "evalf-Int" routines fail to deal with the constructed interpolating function (an object) in as versatile a way as they do regular appliable procedures. I should make a(nother) bug report on it.