It is relatively trivial to replicate the calculation steps described in the article. See the attached

simEq.mw

r[i] is an "indexable" entity, whereas r__i is not. The latter is an "atomic" name. You have two choices

1. use square brackets '[]' for indexable entries throughout
2. if you really need/want to use the double underscore notation in variable names, then you can use these in indexable statements, but you have to play a few tricks, such as 'cat(r__, i)', where the 'i' will be evaluated before being appended to the 'r__' and so the desired name will be produced

Both methods ar illustrated in the attached

subProp.mw

This one uses a bandstop filter with a gaussian characteristic - mainly because the use of "brick-wall" filters can(?) introduce "artefacts" akin to Gibbs phenomena for 1D signals. Frankly I wouldn't know how to recognise such artefacts in a 2D image, so I'd just like to avoid the possiblity.

You can play with the charactersitics of the "gaussian" in this worksheet by adjusting the value of "centrefreq" which defines the centre of the bandstop. I have set this equl to sqrt(15000), as Maxim has done

Other than "crop" the original image, I have tried to keep this as basic as possible, becuase there are various things in Maxim's code which I don't get - not saying they are wrong, just not sure why they exist

One *feature* which requires care is that on performing an inverseTransform, the resutant array will be of type complex, and the entries may/will contain entries with (small) complex parts. It is necessary both to use Re~() to zero out the complex parts, and then force the datatype of the array to be float[8], so that the resulting array will be recognised as a viable image

Is the result better or worse than that given by Maxim's code? - I'd say about the same! I'm not an image expert.

I would confim that commands within the ImageTools() package are absolutely brilliant at causing Maple to lock-up if you get anything mildly wrong - whihc makes writing/debugging/playingAround rather frustrating.

ImProc.mw

(PS code developed in Maple 2017 - haven't ot the time/energy to verify in earlier Maple versions)

If I fix all the parentheses errors which Preben has identified, then I can get plots for all the dependent variables, See the attached

odeProb.mw

but with a few caveats, which are

1. Units in Maple are a bit of mess, and have changed quite a lot in recent versions (ie from Maple 2015 to Maple2017)
2. American units are a bit of a mess. I'm European, so only ever really use SI units - probably makes me biased!
3. Uing palette input, or 2D-input in your worksheet will probably add further difficulties - so I use 1D-input and I type everything - just call me old-fashioned

I would recommend that you set up your own 'system' of units, and including in this system only those units you actually want/need. So, for example, if you don't want "poundals", then don't include it in your system.

The attached defines a system of units called 'CivEng'. This is a pretty "skinny" set of units - you will almost certainly need a few more which I don't know about (I'm not a civil engineer). Calculations are then performed, using this restricted system of units. The only thing I found vaguely disappointing was that the unit-product ksi*inch^2 was not automatically simplified to 'lbf'. I had to use a simplify() command. I guess you can't have everything.:-(

unitProb.mw

Just wind up the numpoinys option in the plot3d coomand an wai a few seconds, as in

restart;
f:=unapply(piecewise(x^2+y^2 <= 1 and x>=0, 4*x*y^2-x^2, undefined),x,y):
plot3d(f,0..1, -1..1, numpoints=100000);

if you convert everytning to 1D input, it works perfectly, as in the following

Download twoDoneD.mw

Aaah - the joys of using 2D input mode

the "similar" command would be

DiscreteTransforms:-FourierTransform()

but I can't check this because you have uploaded a png file of your complete worksheet with input, output, all the prettily formatted text etc.

Could I extract the relevant input data from this- ie the original Saturn piccie. Maybe, but I'm not going to, because it woud be so much eaiser if you just uploaded a png file which contained nothing but the initial data - ie the Saurn piccie

and I'm already getting bored repeating myself -

The answer to the question

I'm aware it can be done, but my question is how I go about 2d-transforming an image in Maple, i.e what would the commands be.

is as follows

1. Maple can import .png images. Unfortunately I cannot meaningfully do this because you do not supply a png image of the starting data - the images you posted are all combined into a single png file, and I do not plan on wasting my evening trying to extract a subsection of this png file, when it ought to be perfectly possible for you to come with a png file for your starting image. If ytou can't manage to do this, then you are going to get precisely nowhere.
2. So if you want to import an image into Maple, then I recommend that you check out the ImageTools:-Read() command. Strangely enough, this is designed to read images from a file into Maple (and png is one of the supported formats). So please tell me - which part of the ImageTools:-Read() command do you not understand?
3. In order to determine whether or not you have succcessfully imported the image into Maple, you can then use the ImageTools:-View() command to display the imported image - is this difficult to understand?
4. Now when I work through the examples on the ImageTools help pages, it seems that png images are stored as 2-D Arrays of hardware floats
5. Once you have successfukky imported an image into Maple, whihc is stored as a 2D array of hardware floats, I would then examine carefully the help for  DiscreteTransforms:-FourierTransform(). Guess what - this will actually generate the Fourier Transform of 1-D data, 2-D data, ..up to 5-D data, and will accept an array of hardware floats as its input

Other than an inability to read the help pages, what exactly is your problem?

Let's break your requirement down into discrete steps

I'm working with Fourier Transform. I'd like to transform a noisy image to the frequency domain, create a grayscale version of the spectrum, mask the dots or lines, threshold it, multiply the binary mask image with the magnitude image and then transform back to the spatial domain. Is this possible using Maple?

1. "I'd like to transform a noisy image to the frequency domain": First of all, does your image have a 2-D (ie "grayscale") interpretation or is it 3D (ie a "color" image). In other words what format is it? Assumiing that it is "grayscale", then the DiscreteTransforms[FourierTransform] command will handle 2-D Fourier transforms, so you ought to be able to produce a 2D frequency domain representation.
2. "create a grayscale version of the spectrum ": Well starting with grayscale-type input, I think you will get this more or les by default from the 2-D Fourier transorm command. It may be necessary to do a little adjustment, perhaps using the ImageTools:-Create() command, and using the output of the 2D FT command as an initialiser.
3. "mask the dots or lines" Frankly not sure what this means: at this stage in the process you have a 2D grayscale image, so I'm not sure what "dots" and "lines" you are talking about
4. "threshold it," The ImageTools package has a Threshold() command, so asuming that you know the thresholding criteria, then you ought to be able to use this to "threshold" an "grayscale" image. I'm assuming that the result of this operation produces the binary mask image you refer to later
5. "multiply the binary mask image with the magnitude image". The ImageTools package also contains a Mask() command, for precisely this purpose
6. "transform back to the spatial domain" The DiscreteTransforms package will allow you to perform the 2D inverse transform, to get back to the spatial domain.

A few observations

1. Image processing in this way can get quite time-consuming, and memory hungry, so whilst developing the code to perform the above process, I'd work with relatively small images.
2. The FT/IFT stages will be speeded up if the FFT algorithm can be used, which requires that pixel counts be a power of two in both dimensions. If this isn't possible then the DFT/IDFT algorithms will be reasonably quick provided that both pixel counts are "highly-composite", preferably with a smal number of low value prime factors (say 2,3,5).
3. The de-noising process you describe is pretty basic - really it is a thresholding operation. Better de-noising, would probably be achieved by applying a 2D low pass filtering operation to the frequency domain image. The diffciulty here would be figuring out which low-pass filter function is optimal, in term of removing the "maximum" amount of noise without introducing too much "blurring" of the desired data

Good luck!

In 2017.3, if you change the Typesetting option from "Extended" to "Maple Standard", then the problem "goes away".

There are several possible calling sequences for the plot3d() command, and it is important to understand all of them.

Some enlightenment can be obtained by stalling the evaluation of the plot3d() command whilst permitting the evaluation of the arguments - because Maple will always "evaluate" the arguments first.

This can be achieved by putting "uneval" quotes around the function -as in writing 'plot3d'() rather than plot3d()

So read the attached *carefully*, to see if it helps

plot3dStuff.mw

Until you have grasped all of the above, there is little point in considering the possibilities of postponing evaluation of the arguments. When you make remarks such as

Knowing that quotes can sometimes solve this sort of problem I then tried:.....

This simply demonstrates that you have no idea what "quotes" are for, why they work, when they work, or anything else. You appear to believe that this is some random thing which "might work", when everything else has failed.

It is actually a better idea to understand what uneval quotes do, when they do it, why they do it, etc

  restart;
  Pe := 10:
  x := 0:
  y := 8/10:
  z := 1/18:
  f := exp(-(xip^2+yip^2)*(1/2))*exp(-((2*Pe*tp+x-xip)^2+(y-yip)^2+z^2)/(4*(t-tp)))/((4*Pi^(3/2))*(t-tp)^(3/2)):
  F := int( f,
            yip = -infinity .. infinity,
            xip = -infinity .. infinity,
            tp = 0 .. t
          );
  plot(F, t=0..0.5);
  objF:=unapply(F,t):
  Optimization:-NLPSolve( objF,
                          0..1/2,
                          maximize=true
                        );

Try the following code

restart;
diagonalise:= proc( A )
                                  uses LinearAlgebra;
                                  local eVs:= Eigenvectors(A),
                                           Tinv:= MatrixInverse(eVs[2]),
                                           Dg:= DiagonalMatrix(eVs[1]);
                                  return A, Dg, eVs[2], Tinv;
                        end proc:
diagonalise(<<2,-1>|<3,-2>>);
diagonalise(<<-2,0,24,3>|<0,-2,-6,0>|<0,0,2,0>|<0,0,0,1>>);

I *think* you are being bamboozled by the fact that Maple treats 'I' as the square root of -1, and processes things accordingly - after all the "normal" way to express/display a complex number would be a+I*b, not a+b*I, even if you wrote the latter.

The only way I can find to get around this, is to redefine the symbol used for the imaginaryUnit (ie the square root of -1) to anything other than 'I', and then apply your requirement.

So the following, will provide what you want from Maple18 to Maple 2017: I can't do earlier Maple versions because I simply don't have them loaded

restart;
interface(imaginaryunit=whatever);
print(lambda*I);