MGL script language

Examples and Hints

Alexey Balakin, 2007-2009

Contents

  1. Basic usage
  2. Working with axes
  3. Text printing
  4. Data handling
    1. Array creation
    2. Data changing
  5. Data plotting
    1. Plots for 1D data
    2. Plots for 2D data
    3. Plots for 3D data
  6. Surface transparency
  7. Hints for MGL scripts
    1. [xyz]range commands and options
    2. Speeding up of the drawing
    3. "Compound" graphics
    4. Two axes in one plot
    5. Titles for the plot
    6. Changing of the color range
    7. Management of the point cutting
    8. Vector field visualization
    9. Several light sources
    10. Cut features
    11. Mapping visualization
    12. Nonlinear fitting
    13. PDE solving
More examples can be found at MathGL site.

Basic usage

The simplest code is

rotate 40 60 # rotate coordinates
box          # draw bounding box

Here bounding box is drawn. You can see that symbol '#' start comment (any character after it is ignored). Also you may note that arguments of commands is delimited one from another by space (or by tabulation).

The script is executed line-by-line and any settings influence only on later commands. This is the matter why I put command rotate 40 60 first (before drawing command) and this is main difference from some other plotting languages.

[contents]

Working with axes

Let me show command of axes transformation. MathGL has the following commands: subplot, inplot, aspect and rotate. The order of its calling is strictly determined. At first one changes the position of axes in image area (commands subplot and inplot). After it one may rotate plot (command rotate). Finally one may change aspects of axes (command aspect). The following code illustrates it:

subplot 2 2 0
box
text -1 1.1 1 'Just box' 'rL'
inplot0.2 0.5 0.7 1
box
text 0 1.2 1 'InPlot example'

subplot 2 2 1
rotate 60 40
aspect 1 1 1
box
text 1 1 1.5 'Rotate only' 'rR'

subplot 2 2 2
rotate 60 40
aspect 1 1 2
box
text 0 0 2 'Aspect and Rotate'

subplot 2 2 3
rotate 60 40
aspect 1 2 2
box
text 0 0 1.5 'Aspect in other direction'

Here I used command text for printing text in arbitrary position of picture. Text coordinates and size is connected with axes. However, text coordinates may be everywhere, including outside the bounding box. I shall show its features later. String arguments is denoted by symbol ' at start and at the end of the string.

MathGL library can draw not only the bounding box but also the axes, grids, labels and so on. The limits of axes and its origin (the point of intersection) are determined by command axis with numerical arguments. Ticks on axis are specified by commands xtick, ytick, ztick. If its first argument is positive then it gives the step between ticks, if it is negative then it gives the number of ticks on the axis, if it is zero then logarithmic ticks is drawn. The second argument give numbers of subticks between ticks.

Command axis with one string or without arguments draws axes. Its textual string (if present) shows in which directions the axis or axes will be drawn (by default 'xyz', command draws axes in all directions). command grid draws grid perpendicularly to specified directions. Example of axes and grid drawing is:

axis -1 -1 -1 1 1 1 # be sure that axis range is
origin 0 0 0        # set origin (point of axes intersection)

subplot 2 2 0
xtick 0.5 3 # set tick step to 0.5 and draw 3 subticks
ytick 0.5 3
box         # should be after change [xyz]tick
axis 'xy'
grid
text 0 1.3 1 'Axis and grid'

xtick -6    # set 6 ticks with step 1/5 of axis range
ytick -6    # setup commands can be in any place

subplot 2 2 1
rotate 60 40
axis        # draw all axes
xlabel 'x'
ylabel 'y'
zlabel 'z'
text 0 0 1.5 'Axis and labels'

subplot 2 2 2
rotate 60 40
xtick 0.2
ytick 0.2
ztick 0.2   # too low step of ticks
axis -1 -1 -1 1 1 1
origin -1 -1 -1
axis
grid
text 0 0 1.5 'Shift origin and add grid'
text 0 0 1.2 '(note, too many ticks)'

subplot 2 2 3
rotate 60 40
xtick -5    # decrease number of ticks
ytick -5
ztick -5
axis 'yz'
ylabel 'Y axis' 0
zlabel 'Z axis' 0
text 0 0 1.5 'Remove X axis, and'
text 0 0 1.2 'decrease number of ticks'

This example shows the importance of correct choosing of the number of ticks on axis. If tick step is too small then its text may overlap and become unreadable. This code has the example of [xyz]label command. It draws label for axis in specified direction. The text position on axis is specified by second argument of [xyz]label command. If it is positive then text is drawn near maximum of axis, if negative then near minimum of axis, if zero then at center of axis.

Next point is the using of curvilinear coordinates. In difference from other systems of plot creation, MathGL uses textual formulas for connection of old (data) and new (output) coordinates. This allows one to plot in arbitrary coordinates. The following code plots the line y=0, z=0 in Cartesian, polar, parabolic and spiral coordinates:

var x 50 -1 1   # set x to change from -1 to 1
var y 50 0.5 0.5# boundaries is the same => all values set to 0.5
new z 50        # new variable is always zero filled
axis -1 -1 -1 1 1 1
origin -1 1 -1

subplot 2 2 0
rotate 60 40
plot x y z 'r2'
axis
grid
text 0 1.3 1 'Cartesian'

subplot 2 2 1
axis 'y*sin(pi*x)' 'y*cos(pi*x)' ''
rotate 60 40
plot x y z 'r2'
axis
grid
text 0 1.3 1 'Cylindrical'

subplot 2 2 2
rotate 60 40
axis '2*y*x' 'y*y - x*x' ''
plot x y z 'r2'
axis
grid
text 0 1.3 1 'Parabolic'

subplot 2 2 3
rotate 60 40
axis 'y*sin(pi*x)' 'y*cos(pi*x)' 'x+z'
plot x y z 'r2'
axis
grid
text 0 1.3 1 'Spiral'

[contents]

Text printing

MathGL prints text by vector Hershey font. There are commands for manual specifying of text position (like text) and for its automatic selection (like [xyz]label, legend and so on). MathGL prints text always in specified position even if it lies outside the bounding box. The default size of font is specified by second argument of font command. However, the actual size of output string depends on position of axes (depends on commands subplot, inplot). Switching of font style (italic, bold, gothic, script and so on) can be done for the whole string (by command parameter) or inside the string. By default MathGL parses TeX-like commands for symbols and indexes. Example of MathGL font drawing is:

text 0 1 0 'Text can be in english и на русском'
text 0 0.6 0 'It can be \\wire{wire} and \\big{big}'
text 0 0.2 0 'One can change style in string: \\b{bold}, \\i{italic, \\b{both}}'
text 0 -0.2 0 'Easy to \\a{overline} or \\u{underline}'
text 0 -0.6 0 'Easy to change indexes ^{up} _{down} @{center}'
text 0 -1 0 'It parse TeX: \\int \\alpha \\cdot \\sqrt3{sin(\\pi x)^2 + \\gamma_{i_k}} dx'

[contents]

Data handling

Variable name is arbitrary combination of symbols (except spaces and ') started from a letter and with length less than 64. It is possible to use subarrays (like in subdata command) as command argument. For example, a(1) or a(1,:) or a(1,:,:) is second row, a(:,2) or a(:,2,:) is third column, a(:,:,0) is first slice and so on. Special variables nan=#QNAN, pi=3.1415926..., on=1, off=0, :=-1 are treated as number if they were not redefined by user. Before the first using all variables must be defined with the help of commands new, var, list, copy, read and so on. Note, that if a variable with the same name exists then it will be overwritten during definition.

Array creation

One can create a variable and fill it by several ways. Let me show it:

[contents]

Data changing

MathGL has command for data processing: differentiating, integrating, smoothing and so on. Let us consider some examples. The simplest ones are integration and differentiation. The direction in which operation will be performed is specified by textual string, which may contain symbols 'x', 'y' or 'z'. For example, call diff a 'x' will differentiate data along x direction; call integrate a 'xy' will do double integration of data along x and y direction; call diff2 c 'xyz' will apply 3d Laplace operator to data and so on. Example of this operations on 2d array a = x*y is presented in code:

new a 30 40
modify a 'x*y'
axis 0 0 0 1 1 1
subplot 2 2 0
rotate 60 40
surf a
box
text 0.7 1 1.2 'a(x,y)'

subplot 2 2 1
rotate 60 40
diff a 'x'
surf a
box
text 0.7 1 1.2 'da/dx'

subplot 2 2 2
rotate 60 40
integrate a 'xy'
surf a
box
text 0.7 1 1.2 '\int da/dx dxdy'

subplot 2 2 3
rotate 60 40
diff2 a 'y'
surf a
box
text 0.7 1 1.2 '\int {d^2}a/dxdy dx'

Data smoothing (command smooth) is more interesting and important. This command has 2 main arguments: type of smooting and its direction. Now 4 methods are supported: 0 – does nothing for delta=0 or approaches data to zero with the step delta, 1 – linear averaging by 3 points (default), 2 – linear averaging by 5 points, 3 – quadratic (parabolic) averaging by 5 points. Let me demonstrate it for 1d case:

new y0 30
modify y0 '0.4*sin(2*pi*x)+0.3*cos(3*pi*x)-0.4*sin(4*pi*x)+0.2*rnd'

copy y1 y0
smooth y1 1
copy y2 y0
smooth y2 2
copy y3 y0
smooth y3 3

plot y0 'k'; legend 'NONE'
plot y1 'r'; legend 'LINE\_3'
plot y2 'g'; legend 'LINE\_5'
plot y3 'b'; legend 'QUAD\_5'
legend
box

Finally one can create new data arrays on base of existing one: extract slice, row or column of data (subdata), summarize along some of direction(s) (sum), find distribution of data elements (hist). Note that all these commands require the definedvariables.

[contents]

Data plotting

Let me now show how to plot data. MathGL generally has 2 types of plotting commands. Simple variant require singles data array for plotting, other data (coordinates) are considered uniformly distributed in interval specified by axis or [xyz]range commands. Second variant requires data arrays for all coordinates. It allows one to plot rather complex multivalent curves and surfaces (in case of parametric dependencies). Arguments setting to default values allows one to plot data in standard form. Manual arguments setting gives possibility for fine tuning of colors, positions and view of graphics. Note that the call of drawing command adds something to picture butdoes not clear the previous plots (as it does in Matlab). Another difference from Matlab is that all setup (like transparency, lightning, axis borders and so on) must be specified before plotting commands.

Plots for 1D data

Term "1D data" means that data depend on single index (parameter) like curve in parametric form {x(i),y(i),z(i)}, i=1...n. There are 5 generally different types of data representations: simple line plot (plot), line plot with filling under it (area), stairs plot (step), bar plot (bars) and vertical lines (stem).Each type of plotting has similar interface. There are 3D version and two 2D versions. One of last requires single array. The parameters of line and marks are specified by the string argument. If the string parameter is '' then solid line with color from palette is used.

Below I shall show the features of 1D plotting on base of plot command. Let us start from sinus plot:

new y0 50
fill y0 'sin(pi*x)'
subplot 2 2 0
plot y0
box

Style of line is not specified in plot command. So MathGL uses the solid line with first color of palette (this is blue). Next subplot shows array y1 with 2 rows:

subplot 2 2 1
new y1 50 2
modify y1 'sin(pi*2*x-pi*y/2)/(1+y)'
plot y1
box

As previously I did not specify the style of lines. As a result, MathGL again uses solid line with next colors in palette (there are green and red). Now let us plot a circle on the same subplot. The circle is parametric curve x=cos(π t), y=sin(π t). I will set the color of the circle (dark yellow, 'Y') and put marks ('+') at point position:

new x 50
modify x 'cos(pi*2*x-pi)'
plot x y0 'Y+'

Note that solid line is used because I did not specify the type of line.

Drawing in 3D space is mostly the same. Let us draw spiral with default line style. Now its color is 4-th color from palette (this is cyan):

subplot 2 2 2
rotate 60 40
new z 50
modify z '2*x-1'
plot x y0 z 'bo '
box

Note that line style is empty ('') here. Usage of other 1D plotting commands looks similar:

subplot 2,2,3
rotate 60 40
bars x y0 z 'r'
box

[contents]

Plots for 2D data

Surfaces surf and other 2D plots (mesh, dens, cont, boxs, axial, tile, fall, belt, contf) are drown the same simpler as 1D one. The difference is that the string parameter specifies not line style but the color scheme of the plot. Here I draw attention on 4 most interesting color schemes. There is gray scheme where color is changed from black to white (string 'kw') or from white to black (string 'wk'). Another scheme is useful for accentuation of negative (by blue color) and positive (by red color) regions on plot (string 'BbwrR'). Last one is the popular "jet" scheme (string 'BbcyrR').

Now I shall show the example of a surface drawing. At first let us switch lightning on

light on

and draw the surface, considering coordinates x,y to be uniformly distributed in bounding box

new a0 50 40
modify a0 '0.6*sin(2*pi*x)*sin(3*pi*y)+0.4*cos(3*pi*(x*y))'
subplot 2,2,0
rotate 60 40
surf a0
box

Color scheme was not specified. So previous color scheme is used. In this case it is default color scheme ("jet") for the first plot. Next example is a sphere. The sphere is parametrically specified surface:

new x 50 40
new y 50 40
new z 50 40
modify x '0.8*sin(2*pi*x)*sin(pi*y)'
modify y '0.8*cos(2*pi*x)*sin(pi*y)'
modify z '0.8*cos(pi*y)'
subplot 2,2,1
rotate 60 40
surf x y z 'BbwrR'
box

I set color scheme to 'BbwrR' that corresponds to red top and blue bottom of the sphere.

Surfaces will be plotted for each of slice of the data if nz>1. Next example draws surfaces for data arrays with nz=3:

new a1 50 40 3
modify a1 '0.6*sin(2*pi*x)*sin(3*pi*y)+0.4*cos(3*pi*(x*y))'
modify a1 '0.6*cos(2*pi*x)*cos(3*pi*y)+0.4*sin(3*pi*(x*y))' 1
modify a1 '0.6*cos(2*pi*x)*cos(3*pi*y)+0.4*cos(3*pi*(x*y))' 2
subplot 2 2 2
rotate 60 40
alpha on
surf a1
box

Note that it may entail a confusion. However, if one will use density plot then the picture will look better:

subplot 2 2 3
rotate 60 40
dens a1
box

Note that the previous color scheme is used in last plots because there are no direct specification of the one.

Drawing of other 2D plots is analogous. The only peculiarity is the usage of flag '#' in string argument. By default this flag switches on the drawing of a grid on plot (grid or mesh for plots in plain or in volume). However, for isosurfaces (including surfaces of rotation axial) this flag switches the face drawing off. Figure becomes wired. The following code gives example of flag '#' using (compare with normal command drawing as in its description):

alpha on
light on
new a 30 20
modify a '0.6*sin(2*pi*x)*sin(3*pi*y) + 0.4*cos(3*pi*(x*y))'

subplot 2 2 0
rotate 40 60
surf a 'BbcyrR#'
box
subplot 2 2 1
rotate 40 60
dens a 'BbcyrR#'
box
subplot 2 2 2
rotate 40 60
cont a 'BbcyrR#'
subplot 2 2 3
rotate 40 60
axial a 'BbcyrR#'
box

[contents]

Plots for 3D data

Drawing procedures for 3D plot looks analogously ones for 1D and 2D plots described above. There are 3 general types of 3D plots: (i) plots on slices or on projections (dens3, cont3), (ii) isosurfaces (surf3), (iii) cloud-like plots (cloud). Plots on slice are clear enough – one specifies a slice (as its index or as value of coordinate) and MathGL draws contour lines or density plot on slice plane. Isosurface gives more information. Isosurface is 3D analog of contour line cont. It shows the region where values of data array exceed specified isosurface level. Plot becomes more informative if one adds transparency, lightning or set color scheme depending on coordinates. Generalization of isosurface is cloud-like plot. This plot draws by darker color and less transparent regions with high values of data. Contrary the regions with low values are transparent. For plotting of the phase of fields (or beams or pulses) one can use isosurface which transparency depends on other data array (see command surf3a). As example of 3D data plots let us draw Gaussian beam diffraction in space. Beam propagates along x axis:

alpha on
light on
new a 30 30 30
new b 30 30 30
modify a 'exp(-16*((z-0.5)^2+(y-0.5)^2)/(1+4*x^2))'
modify b '16*((z-0.5)^2+(y-0.5)^2)*(x)/(1+4*x^2)'
crange 0 1

subplot 2 2 0
rotate 40 60
surf3 a 'wgk'
box
subplot 2 2 1
rotate 40 60
densa a
box
subplot 2 2 2
rotate 40 60
cloud a
box
subplot 2 2 3
rotate 40 60
surf3a b a 'q'
box

[contents]

Surface transparency

MathGL library has advanced features for setting and handling the surface transparency. The simplest way to add transparency is the using of command alpha. As result all further surfaces (and isosurfaces, density plots and so on) become transparent. However, its look can be additionally improved.

First, the selected surface will be non-transparent if one sets transparento ff before the surface drawing and sets it on after the drawing.

Second, the value of transparency can be different from surface to surface. To do it just change the value of alphadef before the drawing of selected surface. If its value is close to 0 then the surface becomes more and more transparent. Contrary, if its value is close to 1 then the surface becomes practically non-transparent. This is some analog of transparent off. One may use option alpha here to change alphadef only for one graphics:

alpha on
surf a          # this plot is transpared
surf b; alpha 1 # this plot is not transpared
surf c          # and this plot it transpared again

Third feature is the changing of the way how the light goes through overlapped surfaces. The command transptype defines it. By default the usual transparency is used (transptype 0) – surfaces below is less visible than upper ones. A "glass-like" transparency (transptype 1) has a different look when the surface just decreases the background light (the surfaces are commutable in this case). A "neon-like" transparency (transptype 2) has more interesting look. In this case a surface is the light source (like lamp on dark background) and just adds some intensity to color. At this the library automatically sets black color for background and changes default line color to white.

[contents]

Hints for MGL scripts

In this section I have included some small hints and advices for the improving of the quality of plots and for the demonstration of some non-trivial features of MathGL library. In contrast to previous examples here I showed mostly the idea but not the whole drawing command.

[xyz]range commands and options

The usual way to change bounding box is the using of axis command. But if you don't know the range of data or you want to change only one dimension (for example x) then the [xyzc]range command is a solution. This command has a variant for the manual rage setting (for example, xrange 0 1) and has a variant for automatic data rage determining (for example, zrange a).

It is much more interesting to use these commands as options. In these options one should define the range of auto-filled variables. For example, command plot y; xrange 0 2 will plot data y with x-coordinate distributed in range [0, 2]. At this, the range may exceed the bounding box (so that only part of data will be plotted).

One more interesting feature is that auto-filled coordinate vary from the first argument of the option to the second one. This allows to make mirrored plots:

new y 50
modify y 'x^2'
new a 30 40
modify a 'exp(-x^2-4*(2*y-1)^2)'

subplot 2 1 0
plot y 'b'             # normal plot
plot y 'r'; xrange 0 1 # plot in some part of graphics
plot y 'g'; xrange 1 0 # the same but mirrored
box

subplot 2 1 1
rotate 50 30
light on
surf a 'r'; xrange 0 1 # one part of pulse
surf a 'b'; xrange 0 -1# second part

Speeding up of the drawing

The speeding up of the drawing is important task when the script is used often (for example, in UDAV program). There are several things that can be done.

"Compound" graphics

As I have noted above, MathGL commands (except the special one, like clf) do not erase the previous plotting but just add the new one. It allows one to draw easily "compound" plots. For example, popular Matlab command surfc can be emulated in MathGL by 2 calls:

surf a
cont a 0 7 -1    # draw contours at z = -1

Here a is 2-dimensional data for the plotting, -1 is the value of z-coordinate at which the contour should be plotted (at the bottom in this example). Analogously one can draw density plot instead of contour lines and so on.

Another nice plot is contour lines plotted directly on the surface:

light on
surf a 'BbcyrR'  # select 'jet' colormap for the surface
cont a 'y'       # and yellow color for contours

The possible difficulties arise in black&white case, when the color of the surface can be close to the color of a contour line. In that case I may suggest the following code:

light on
surf a 'kw'      # select 'gray' colormap for the surface
crange -1 0      # first draw for darker surface colors
cont a 'w       # white contours
crange 0 1       # now draw for brighter surface colors
cont a 'k'       # black contours
crange -1 1      # return color range to original state

The idea is to divide the color range on 2 parts (dark and bright) and to select the contrasting color for contour lines for each of part.

Similarly one can plot flow thread over density plot of vector field amplitude (this is another amusing plot from Matlab) and so on. The list of compound graphics can be prolonged but I hope that the general idea is clear.

Two axes in one plot

Developing the previous hint one can make a plot with 2 or more axes. The idea is that the change of settings does not influence on the already drawn graphics. So, for 2-axis plot let us set the first axis and draw everything concerning it. Then let us setup the second axis and draw things for the second axis. The corresponding code is:

# set up first axis
axis -1 -1 -1 1 1 1
origin -1 -1 -1
axis             # draw it
plot y1 'b'      # draw something in first axis<
# set up second axis
axis 0 0 0 1 1 1
origin 1 1 1
axis             # draw it
stem y2 'r'      # draw something in second axis<

Note, that the first and the second axes look better if being placed in different corners. In the code presented above the first axis is placed in the left-bottom corner, and the second one is placed in the right-top corner.

Titles for the plot

The printing of nice titles for the plot is not so trivial task in general case. The problem is that rotation and aspect change lead to different looks for titles of different subplots. So, the resulting look is not so good as it could be. The solution is simple – to print titles exactly after subplot call and before any rotation, aspect change and so on! Analogously, the title for the whole picture looks better if it is printed first (before any subplot calls).

Changing of the color range

By default (for the user comfort), the color range is set equal to z-range of the plot. However, there are different ranges. So, one can obtain an amusing plot by the change of color range manually. For example, there are plots with one-color bottom (or top) or practically bi-color picture and so on. For example, compare 2 surfaces:

subplot 2 1 0
axis -1 -1 -1 1 1 1
surf a           # usual coloring range<
subplot 2 1 1
crange 0 1
surf a           # bottom of the surface have one-colour filling

Management of the point cutting

Sometimes an experimental or numerical surface has outstanding points. Visualization of such surface will lead to the hole(s) in place of such points. The standard method of "fighting" – tochange data values – is not always good and is not so convenient. MathGL library has another method – to set cut off. As a consequence, all outstanding point will be projected on the bounding box.

Such method is good not only for outstanding points but also for the case when one need to plane the bottom or the top of the plot. Exactly such case is demonstrated in the code:

new a 20 30      # create some data
modify a '0.6*sin(2*pi*x)*sin(3*pi*y) + 0.4*cos(3*pi*(x*y))'
# set lower border above the data minimal value
axis -1 -1 0 1 1 1
cut off          # set off cutting flag
surf a           # and draw the surface

It is an interesting result, is not it? Note that the same result can be achieved by options cut:

surf a; cut off  # the same surface

Vector field visualization

Vector field visualization (especially in 3d case vect or vectc) may look tangly – there are too many overlapping lines. I may suggest 2 ways to solve this problem. The first one is to change meshnum the number of hachures for decreasing. The second way is to use the flow thread chart flow. Unfortunately, I don't know any other methods to visualize 3d vector field. If you know any, e-mail me and I shall add it to MatGL.

Several light sources

In contrast to the most of the other programs, MathGL supports several (up to 10) light sources. Moreover, the color each of them can be different: white (this is usual), yellow, red, cyan, green and so on. The use of several light sources may be interesting for the highlighting of some peculiarities of the plot or just to make an amusing picture. Note, each light source can be switched on/off individually.

Cut features

MathGL library has a feature for cutting of points in some region cut. Such an excision can be used to improve the look of the graphics. Moreover, this cutting may help to show an internal structure of an object (like isocaps plot in Matlab). For example, let us use the standard data array and show its interior:

new c 61 51 40   # create the data
crange -0.5 1
modify c '(-2*((2*x-1)^2 + (2*y-1)^2 + (2*z-1)^4 - (2*z-1)^2 - 0.1))'
cut 0 -1 -1 1 0 1.1

rotate 40 60
surf3 c -0.5 'BbcyrR'
contf3 c 'x' -1 'BbcyrR' 10
contf3 c 'y' -1 'BbcyrR' 10
contf3 c 'z' 0 'BbcyrR' 10
contf3 c 'z' 39 'BbcyrR' 10

Mapping visualization

Sometime ago I worked with mapping and have a question about its visualization. Let me remind you that mapping is some transformation rule for one set of number to another one. The 1d mapping is just an ordinary command – it takes a number and transforms it to another one. The 2d mapping (which I have used) is a pair of commands which take 2 numbers and transform them to another 2 ones. Except general plots (like surfc, surfa) there is a special plot – Arnold diagram. It shows the area which is the result of mapping of someinitial area (usually square).

I tried to make such plot in map. It shows the set of points or set of faces, which final position is the result of mapping. At this, the color gives information about their initial position and the height describes Jacobian value of the transformation. Unfortunately, it looks good only for the simplest mapping but for the real multivalent quasi-chaotic mapping it produces a confusion. So, use it if you like :).

[contents]

Nonlinear fitting

Nonlinear fitting is rather simple. All that you need is the data to fit, the approxiamtion formula and the list of coefficients to fit (better with its initial guess values). Let me demonstrate it on the following simple example. First, let us use sin function with some random noise:

new rnd 100 #data to be fitted
modify rnd '0.4*rnd+0.1+sin(4*pi*x)'
new in 100 #data with ideal function (without noise)
modify in '0.3+sin(4*pi*x)'

and plot it

yrange -2 2
plot rnd '. '
plot in 'b'
box
text 0 2.2 'initial: y = 0.3+sin(2\pi x)' 'C:b' -1

The next step is the fitting itself. For that let me specify an initial values ini for coefficients 'abc' and do the fitting for approximation formula 'a+b*sin(c*x)'

list ini 1 1 3
fit res rnd 'a+b*sin(c*x)' 'abc' ini
Now display it
plot res 'r'
text -1 -1.3 'fitted:' 'L:r' -1
putsfit 0 -1.8 'y = ' 'C:r'

NOTE! the fitting results may have strong dependence on initial values for coefficients due to algorithm features. The problem is that in general case there are several local "optimums" for coefficients and the program returns only first found one! There are no guaranties that it will be the best. Try for example to set list ini 0 0 0 in the code above.

[contents]

PDE solving

Solving of Partial Differential Equations (PDE, including beam tracing) and ray tracing (or finding particle trajectory) are more or less common task. So, MathGL have several commands for that. There are ray for ray tracing, pde for PDE solving, qo2d for beam tracing in 2D case. These functions take "Hamiltonian" or equations as string values.

The ray tracing can be done by ray command. Really ray tracing equation is Hamiltonian equation for 3D space. So, the function can be also used for finding a particle trajectory (i.e. solve Hamiltonian ODE) for 1D, 2D or 3D cases. The function have a set of arguments. First of all, it is Hamiltonian which defined the media (or the equation) you are planning to use. The Hamiltonian is defined by string which may depend on coordinates 'x', 'y', 'z', time 't' (for particle dynamics) and momentums 'p'=px, 'q'=py, 'v'=pz. Next, you have to define the initial conditions for coordinates and momentums at 't'=0 and set the integrations step (default is 0.1) and its duration (default is 10). The Runge-Kutta method of 4-th order is used for integration.

define $1 'p^2+q^2-x-1+i*0.5*(y+x)*(y>-x)'
ray r $1 -0.7 -1 0 0 0.5 0 0.02 2
plot r(0) r(1) 'k'

This example calculate the reflection from linear layer (media with Hamiltonian 'p^2+q^2-x-1'=p_x^2+p_y^2-x-1). This is parabolic curve. The resulting array have 7 columns which contain data for {x,y,z,p,q,v,t}.

The solution of PDE is a bit more complicated. As previous you have to specify the equation as pseudo-differential operator H(x, ∇) which is called sometime as "Hamiltonian" (for example, in beam tracing). As previously, it is defined by string which may depend on coordinates 'x', 'y', 'z' (but not time!), momentums 'p'=(d/dx)/i k0, 'q'=(d/dy)/i k0 and field amplitude 'u'=|u|. The evolutionary coordinate is 'z' in all cases. So that, the equation look like du/dz = ik0 H(x,y, p, q, |u|)[u]. Dependence on field amplitude 'u'=|u| allows one to solve nonlinear problems too. For example, for nonlinear Shrodinger equation you may set ham='p^2 + q^2 - u^2'. Also you may specify imaginary part for wave absorption, like ham = 'p^2 + i*x*(x>0)', but only if dependence on variable 'i' is linear (i.e. H = Hre+i*Him).

Next step is specifing the initial conditions at 'z'=Min.z. The function need 2 arrays for real and for imaginary part. Note, that coordinates x,y,z are supposed to be in specified range [Min, Max]. So, the data arrays should have corresponding scales. Finally, you may set the integration step and paramter k0=k0. Also keep in mind, that internally the 2 times large box is used (for suppressing numerical reflection from boundaries) and the equation should well defined even in this extended range.

Final comment is concerning the possible form of pseudo-differential operator H. At this moment, simplified form of operator H is supported – all "mixed" terms (like 'x*p'->x*d/dx) are excluded. For example, in 2D case this operator is effectively H = f(p,z) + g(x,z,u). However commutable combinations (like 'x*q'->x*d/dy) are allowed for 3D case.

So, for example let solve the equation for beam deflected from linear layer and absorbed later. The operator will have the form 'p^2+q^2-x-1+i*0.5*(z+x)*(z>-x)' that correspond to equation ik0z u + Δ u + x ⋅ u + i (x+z)/2 ⋅ u = 0. This is typical equation for Electron Cyclotron (EC) absorption in magnetized plasmas. For initial conditions let me select the beam with plane phase front 'exp(-48*(x+0.7)^2)'. The corresponding code looks like this (see section PDE sample):

new re 128
new im 128
fill re 'exp(-48*(x+0.7)^2)'
pde a 'p^2+q^2-x-1+i*0.5*(z+x)*(z>-x)' re im 0.01 30
transpose a
crange 0 1
dens a 'wyrRk'

The last example is example of beam tracing. Beam tracing equation is special kind of PDE equation written in coordinates accompanied to a ray. Generally this is the same parameters and limitation as for PDE solving but the coordinates are defined by the ray and by parameter of grid width w in direction transverse the ray. So, you don't need to specify the range of coordinates. BUT there is limitation. The accompanied coordinates are well defined only for smooth enough rays, i.e. then the ray curvature K (which is defined as 1/K2 = (|d2r/dt2|2 |dr/dt|2 - (d2r/dt2, dr/dt)2)/|dr/dt|6) is much large then the grid width: K>>w. So, you may receive incorrect results if this condition will be broken.

You may use following code for obtaining the same solution as in previous example:

new re 128
new im 128
new xx
new yy
fill re 'exp(-48*x^2)'
qo2d a $1 re im r 1 30 xx yy
crange 0 1
dens xx yy a 'wyrRk'

[contents]