@C_R The options that I used are documented on the pages ?plot,axis, ?plot,tickmarks, and ?plot,typesetting.

@Fereydoon_Shekofte I don't see the connection. The matrices in this Question don't have any arithmetic properties, their entries aren't necessarily numbers, they aren't necessarily square matrices, and only rowwise, not columnwise, properties are considered.

@planetmknzm 

Let's say that you've plotted 2 spacecurves (on the unit sphere) and extracted their data arrays as explained at the beginning of this thread. Let's say those arrays are A1 and A2. The following code will give you the 9 (arbitrary number; use however many you want) pairs of indices [i, j] such that the angle between A1[i] and A2[j] is as close to 30 degrees as can be achieved with the discrete (i.e., finite) number of points:

COS:= evalf(cos(30*(Pi/180))): #cos(30 degrees)
[lhs]~(sort([op(2, abs~(A1.A2^%T -~ COS))[]], key= rhs))[..9];
   [[178, 187], [14, 23], [62, 121], [80, 139], [65, 127],
    [74, 136], [15, 191], [10, 186], [29, 39]]

#Check angle of 1st pair:
arccos(A1[178].A2[187]^%T)*evalf(180/Pi);
                          30.00128868

The above computation takes 0.7 seconds on my computer.

I don't know if you've upgraded your Maple yet. Depending on your version, the key option to sort may not work. If that's the case, I can give you an easy workaround.

@planetmknzm Which are these scenarios is closest to your situation?

1. You have a finite list or set of points on the sphere, let's say the points generated for the 2 space curves. You want all pairs of these points that are close to 30 degrees apart.
2. You have a finite list or set of points on the sphere, the points generated for the 2 space curves. You want all pairs of these points that are close to 30 degrees apart containing 1 point from each curve.
3. You have two curves that can be parameterized by a real parameter, let's say C1(t) and C2(s). Given a particular value of t, you want to find a value or values of s such that C1(t) and C2(s) are exactly 30 degrees apart.

Options 1 and 2 are crude. I would only use them if I couldn't get the curves continuously parameterized.

@vv Thanks for the correction.

@planetmknzm Given two points on the unit sphere, the central angle (angle at the origin) between them is the arccos of their dot product. If the points are represented as lists P and Q, that's 

arccos(<P>^%T  . <Q>)

There is no picture attached to your Question.

If you have a Maple worksheet, it'd be more helpful if you uploaded that. Use the green up arrow on the toolbar.

@planetmknzm Thank you. I finished the explanation. Let me know if you have any questions.

@planetmknzm 

The geometric simplicity of the surfaces being plotted allows some of the commands that I used to become deceptively simple. The first such simplicity is:

• Given a point P on a sphere centered at the origin, if P is viewed as a vector v, then v is a normal vector to the sphere at P and hence also a normal vector to the plane tangent to the sphere at P.

I'll detail the other simplicity later.

Let's look at the details of an animate command:

plots:-animate(PP, [A(t)], t= a..b, frames= n)

• PP is a procedure that returns a plot.
• A(t) is the arguments to PP parameterized by t.
• a and b are real numbers, the range of values for t.
• frames is a keyword, i.e., it appears literally.
• n is a positive integer, the number of evenly spaced values of t that'll be used.

My animate command was 

plots:-animate(Plane, [C(t)], t= -Pi..Pi, frames= 100, background= Static)

C(t) is my command that generates points on the curve on the sphere. So, 100 such points are passed to Plane. Each point is a list of three real numbers.

Now let's look at Plane:

Plane:= proc(P)
local v:= <P>^%T, B:= LinearAlgebra:-NullSpace(v)^~%T;
    plots:-polygonplot3d(<v+B[1], v+B[2], v-B[1], v-B[2]>)
end proc

A list can be turned into a column vector by enclosing it in angle brackets: <P>. The transpose operator ^%T converts column vectors into row vectors and vice versa. So, v is a row vector, the normal vector of the plane. The distinction between lists, column vectors, and row vectors is of very little mathematical significance; these conversions are just for syntactic reasons.

When NullSpace is passed a row vector v, it returns a basis for the subspace of vectors that are orthogonal to v. In this case, that means the vectors in the tangent plane. The basis is returned as a set of column vectors. I want them to be row vectors for the next command. By adding the elementwise symbol ~ to the transpose operator ^%T (so that the operator is now ^~%T), the entire set is converted to row vectors. So B is now a set of 2 row vectors orthogonal to v.

Now here is the second hidden simplifying detail:

• The NullSpace command returns an orthonormal basis, that is, the basis vectors are orthogonal to each other and they have length 1.

So, adding and subtracting the basis vectors to the point of tangency gives the four corners of a square. That square has diagonal length 2, thus side length sqrt(2). This seemed to be a good size for the animated squares, but you could change it by multiplying the vectors in B by a scaling factor.

The command polygonplot3d accepts a three-column matrix each row of which represents a vertex of the polygon. (The rows are viewed as points, not vectors.) The angle bracket syntax <..., ..., ..., ...> stacks the row vectors into a column, which is the same thing as a matrix.

@Ronan In the following example, note how the invisible function application operator distributes over the comma sequencing operator between the f and g:

(f,g)(x,y);
                        f(x, y), g(x, y)

@ivar_kjelberg You wrote:

• I got the impression that "init" was the old way, and "ModuleLoad" the more recent implementation.

It is true that init is older than ModuleLoad, but, as I said, they are independent. They do different things, and ModuleLoad is not a replacement for init; nor has init been labeled as deprecated. They are different because ModuleLoad() is triggered by the module (which is not necessarily a package) being extracted from an archive, whereas init() is specific to package modules and is triggered specifically by the with command.

You did ask specifically about the with command, so, to directly Answer your title Question:

• Does a ModuleLoad proc() run with the "with()" command?

Direct Answer: No, ModuleLoad() is not necessarily invoked by the with command, but init() is.

• I even cannot find any documentation mentionning the init for Maple 2022.

Yes, it's hard to find. It's the 11^th bullet point in the Description section of help page ?with:

• Prior to arranging for the names exported by a package to be available at the top level, with executes either or both of two procedures that have, for historical reasons, come to be known as the system level initialization and user level initialization routines for a package, provided that they exist. A system level initialization routine is any procedure that is assigned to the name P/init, where P is the name of the package. It is called with no arguments. The user level initialization routine is a package export called init. It is called (if it exists) with no arguments after first calling the system level initialization routine. The system level initialization procedure is not called if the name P/Initialized has the value true, and this name is assigned the value true once the routine has been successfully executed. This allows you to put any once-per-session setup that needs to be done in the system level initialization, and anything that must be run each time with is called on your package can be performed in the user level initialization.

@Nikol Clearly there is a bug here. I'm not recommending that you use expand; that's too ad hoc to be generally useful. However, you're rightfully curious about why expand works. It's because it converts exp(-T) to 1/exp(T) and exp(-2*T) to 1/exp(T)^2. The expansion of the polynomial parts is irrelevant. Letting Y= exp(T) and clearing denominators allows the equation to be rewritten as a polynomial that is quadratic in two variables:

AY:= (20 + 20*T + 2*T*(T+1))*Y - 10 - (2*T + 10.0)*Y^2 = 0:

solve({AY, Y=exp(T)});
             {T = 1.411454823, Y = 4.101918631}, {T = 0., Y = 1.}

Once again, this is too ad hoc, so I'm not recommending this in general. I'm simply presenting it as a curiosity.

@rlopez The connect-the-dots phenomenon that you mentioned is often an issue with the plot command. (In Maple 2022, the new option setting adaptive= geometric, which is now the default, has largely corrected this.) But the case at hand uses implicitplot, and inspection of the plot structures reveals that that phenomenon is not the issue here. Two of the shown plots have a vertical line. (They are truly vertical lines in these cases, unlike the nearly vertical lines produced by connect-the-dots.) Looking at the plot structures for the x= -6..9 and x= -6..11 cases shows that they each contain three explicitly separate CURVES, one of which is (in each case) a vertical line composed of 201 points. In the x= -6..9 case, that vertical line is x = 1.96569296691827; in the x= -6..11 case, it's x= 2.05767071950135.

I don't know why these vertical lines are drawn. Their points are not even close to being zeros of the input expression y - x/(2-x). But it's clear that it's the high-level Maple-language plotting command that's making the error, not the plot renderer.

@fkohlhepp Option numeric is for int, not Int. Your original Question didn't show your integration commands, so I didn't know that you were already using evalf(Int(...)), which is pretty much equivalent to int(..., numeric) anyway. So, you're already using numeric integration, so switching to int(..., numeric) should "work" superficially in the sense of not giving an error, but it won't help with the computation time. So, there's no point in you using numeric.

Given that you're already using numeric integration, the slowness may be due to one of these causes:

1. Singularities: The integrand has singularities within the interval of integration (essentially, values of r where division-by-zero happens or where the integrand is not differentiable). These might still be integrable, but the algorithms to do it numerically are more complicated.
2. Loss of precision: The integrand is complicated enough that the integrator is having trouble evaluating it numerically at sufficiently high precision to get the required precision for the final integral.
3. Massively complicated: The integrand is simply so complicated that each numeric evaluation takes a long time not due to precision but simply due to the tens-of-thousands of steps required.
4. Inadvertent symbolic computations: Due to inadvertent coding of your procedure trq__s, it may be needlessly redoing symbolic computations for each numeric call. One of the most-common cases of this is a procedure which computes a derivative before evaluating it numerically.

Do you have an idea that one of those is the culprit? Showing an example plot of the integrand over the interval of integration and recording the time required to generate that plot will help.

Using some or all of these might help:

1. Put the Re(...) inside the integral.
2. Add the option epsilon= 0.5e-5 to the integration command. This will reduce the precision (precisely, increase the maximum allowed relative error) (of the final result, not of the internal computations) to 5 digits. That number can be adjusted. Even if this level of precision is unacceptable for your purpose, doing this just as a test will help to track down the source of the slowness. 
3. Add one (and only one) of these options to the integration command: method= _d01ajc or method= _d01akc. These are high-speed externally compiled numeric integrators from the Numerical Algorithms Group (NAG). 

@C_R The ' 'arctan' '(algebraic) is a structured type which is essentially used as the 2nd argument of an indets command. That indets command searches for function calls with function name arctan having exactly one argument of type algebraic (see ?type,algebraic).

Since arctan is an actual executable procedure (not just a function), specifying arctan(algebraic) would cause arctan to execute, which is potentially undesirable. (It actually wouldn't be a problem in this particular case because it would return literally arctan(algebraic) unevaluated. But why waste the execution time required for that?). The quotes prevent that execution. One pair lets the expression pass through rcurry, and the second lets it through evalindets.

To understand what the builtin evalindets does, it's useful to rewrite it in top-level Maple. When it has three arguments and no options, it's essentially the same as this:

EvalIndets:= proc(e, T::type, f)
local r:= e, k, Ind:= sort([indets(e, T)[]], 'key'= length);
    for k to nops(Ind) do
        (r,Ind):= eval([r, Ind], Ind[k]= f(Ind[k]))[]
    od;
    r
end proc
:

Now we can trace its execution for educational purposes:

`convert/atan2`:= rcurry(
    EvalIndets, ''arctan''(algebraic), arctan@op@[numer,denom]@op
):
eq:= phi = arctan(A/B) - arctan(arctan(A/B)/B):
trace(EvalIndets):
convert(eq, atan2);

You should study the output of that. It's only 11 lines.