@digerdiga

1. This clearly looks like a different function, it slopes in the CO direction, not in the NO2 direction. If you made any changes to the code, upload your worksheet. rcurry is just syntactic sugar here.

2. A different tilde. You want ?elementwise. Also, look for showassumed in ?assume.

3. Could have been expand, but expanding (a+b)*x is redundant if I want the powers of x. normal is not a keyword here, it's a function (the four-argument form in ?collect). Could have been simplify instead of normal.

4. You can actually do the computation and you'll see what's going on. You're right, you can do the balancing at p=3/2. And also at p=1. We had to choose one of the roots at each step. Say we chose c1*NO2 instead of c2*NO2. Balancing at p=1 for the second time will give you 0 and (c2-c1)*NO2, i.e., it will offer you a different branch.

@vv

According to Littlewood, that was said of Jordan's writings :).

Yes, there is a small parameter NO2 and a fixed parameter CO, it can be improved but not by much.

I agree about the example itself, even just scanning through all the monomials isn't easy, but that's part of the point. A toy example can be immediately done with algcurves:-puiseux.

I added another answer here, using a different method and adding some explanation about how we decide which branch to take. Also fixed the use of convert(..., rational), which should have been convert(..., rational, exact).

@vv

I'm not saying it's wrong, only that it's not documented. The question is what exactly evalf[10](expr) does. Normally procedures evaluate the arguments first, and then the body of the procedure executes. It doesn't seem to be mentioned anywhere that that's not how evalf works. In fact, the example with evalf/NewFunction at ?evalf/details gives an exactly opposite impression: my interpretation of why evalf(NewFunction(Pi)) doesn't work would be that NewFunction(Pi) is already evaluated before evalf gets to it. But the example with series(f(1.*s)) shows that it's more complicated than that: maybe evalf just changes Digits locally, but it does something before evaluating the arguments.

I only just now realized that the way to avoid the magenta message of death is to use try/catch.

@NorwegianStudent

Well, Christopher2222's question is about creating a very specific mask, which is why that was much easier to work with. Embed(abs(fft2)^(1/4)) before applying the mask clearly shows the eight bright dots. The same after applying the mask shows a dark ring where the dots used to be. To roughly determine the positions of the dots, use Preview (in theory, you can also right-click on the preview and set Manipulator to Point Probe, but I couldn't get it to show the coordinates). You can use an image editor to get the coordinates or write an Explore to adjust the mask interactively.

You can do the same for your image. Normalize fft and look at PlotHistogram. Embed(abs(fft)^(1/2)) shows two bright dots on the vertical midline. Changing 14000<e<16000 to 500<e<700 creates a ring overlapping those points. As a result, I get something very blocky and with noisy gridlines, but with the horizontal dark stripes gone. If it's not good enough, you'll have to be more specific about which areas you want to apply a mask to.

Also, wrap everything in try/catch, because otherwise you risk losing all your work if Maple freezes trying to print a huge error message:

try Embed(Scale(abs(fft2)^(1/2), 4)) catch: lasterror end try

@vv

Sure, I'm just saying that (in 2017.3) (Digits:=30;evalf[10](e)) cannot be used for computing at Digits=30 and then truncating, but evalf[10](evalf[30](e)) can.

@Carl Love

Curiously, it does work for the case at hand:

evalf[10](evalf[30](series(f(1.*s), s = 0, 1)));
                              O(s)

So it seems as though evalf[10] sets up a context overriding the global setting of Digits, and evalf[30] sets up a nested context.

@acer

After looking at the userinfo output, I maybe probably see what you mean. evalf/Int integrates a series expansion of f(1.*s) on a subinterval. That looks like a gamble, because for some settings of Digits the 1/s term is going to disappear, for other settings there will be a spurious non-zero coefficient (no matter how high Digits is set).

And evalf/Int rationalizes the expression at some point, which doesn't help, because then there is always a 2*I/Pi-hugerational coefficient at 1/s.

Switching to another method (_Dexp) helps.

By the way, I was surprised by this:

Digits := 30:

evalf[10](series(f(1.*s), s = 0, 1));
                                     (1.*10^(-10)*I)/s+O(s)

evalf[10]~([series(f(1.*s), s = 0, 1)]);
                                           [O(s)]

Didn't know that evalf has non-standard evaluation rules too.

P.S. Maple evaluates the integral of s^p*cos(alpha*s)*gr(beta, s), but the result seems to be wrong for every p where the integral exists (-2<p<1/2):

int(s^p*cos(alpha*s)*gr(beta, s), s = 0 .. infinity);
        -(1/8)*GAMMA(1+(1/2)*p)*(8*2^(-3+p)*beta^(-3-p)*hypergeom([(1/2)*p, 1+(1/2)*p], [1/2],
        alpha^2/beta^2)*cos((1/2)*Pi*p)*GAMMA((1/2)*p)*beta*alpha^2*GAMMA(-(1/2)*p-1/2)-8*2^
        (-3+p)*beta^(-3-p)*hypergeom([(1/2)*p, 1+(1/2)*p], [1/2], -alpha^2/beta^2)*GAMMA((1/2)*
        p)*beta*alpha^2*GAMMA(-(1/2)*p-1/2)+I*2^p*(alpha^2)^(-(1/2)*p)*Pi^(3/2)*hypergeom([1+
        (1/2)*p, 3/2+(1/2)*p], [2], beta^2/alpha^2))/(beta*alpha^2*GAMMA(-(1/2)*p-1/2)*Pi)

@vv

I took this formula and coded it in Maple.

@Christopher2222

Embed(fft2^(1/4)) gives a rough idea. Before applying the mask, the circle of the bright dots is clearly visible.

In fact, since we're not doing any complicated normalizations, the code can be simplified to just this:

# careful, this time fft2 is complex-valued
fft2 := CircularShift(CircularShift(fft, floor(h/2), floor(w/2))*mask, ceil(h/2), ceil(w/2));
img2 := Create(Re~(SignalProcessing:-InverseFFT(fft2)), fit);

@digerdiga

What I'm saying is that (x^p)::anything works in Maple 2017.3, but you probably don't want to use it in your code. You want x^(p::anything).

In 2D input, this should be written as x^p::anything, that is, ::anything should go in the superscript.

@digerdiga

Essentially you're correct, the situation is like this: if p>3/2, then the dominant term is NO2^5. Then

terms := (simplify@{op}@expand@P)(NO2, 600*NO2+alpha*NO2^p) assuming p > 0:

applyrule([NO2^5 = NO2^5, NO2^p::anything = 0], terms);

[solve](add(%), alpha);
                               []

If p<3/2, then the dominant term is NO2^(2p+2). Now

applyrule([NO2^(2*p+2) = NO2^(2*p+2), NO2^p::anything = 0], terms);

[solve](add(%), alpha);
                             [0, 0]

It just means that for p other than 3/2, an expansion of the form 600*NO2+alpha*NO2^p cannot be constructed.

Regarding applyrule, there are more details in ?patmatch and ?type. In 2017.3 this actually works, but it treats x^p as a single pattern variable:

applyrule([l::list = l, (x^p)::anything = [x, p, x^p]], a^2); # l::list to avoid infinite recursion
                          [x, p, a^2]

So x and p are still free variables, but x^p as a whole is matched to a^2.

@digerdiga

Sure it does. It's exactly as you say, the lowest order is 4, and it has to vanish:

coeff(P(NO2, A*NO2), NO2, 4);

solve(%, A);
             600, 600, -228000000000000000/(4071*CO+35020000000000000),
             -57000000000000000/(4071*CO-17795000000000000)

I used applyrule because coeff works only with integer powers. If you look at the result of simplify@{op}@expand, all elements have the form k*NO2^p, there is no free term. I want to replace each term with p. applyrule(x^p::anything=p, 3*x^2) would give 3*2=6, not 2.