## Puzzle - cut it into 2 equal parts

by:

The following puzzle prompted me to write this post: "A figure is drawn on checkered paper that needs to be cut into 2 equal parts (the cuts must pass along the sides of the squares.)" (parts are called equal if, after cutting, they can be superimposed on one another, that is, if one of them can be moved, rotated and (if need to) flip so that they completely coincide) (see the first picture below).
I could not solve it manually and wrote a procedure called  CutTwoParts  that does this automatically (of course, this procedure applies to other similar puzzles). This procedure uses my procedure  AreIsometric  published earlier  https://www.mapleprimes.com/posts/200157-Testing-Of-Two-Plane-Sets-For-Isometry  (for convenience, I have included its text here). In the procedure  CutTwoParts  the figure is specified by the coordinates of the centers of the squares of which it consists).

I advise everyone to first try to solve this puzzle manually in order to feel its non-triviality, and only then load the worksheet with the procedure for automatic solution.

For some reason, the worksheet did not load and I was only able to insert the link.

## First equilibrium condition with Maple

Maple 2020

With this application our students of science and engineering in the areas of physics will check the first condition of balance using Maple technology. Only with entering mass and angles we obtain graphs and data for a better interpretation.

First_equilibrium_condition.zip

Lenin AC

## Add Reverb to Audio

by: Maple 2020

So here's something silly but cool you can do with Maple while you're "working" from home.

• Record a few seconds of your voice on a microphone that's close to your mouth (probably using a headset). This is your dry audio.
• On your phone, record a single clap of your hands in an enclosed space, like your shower cubicle or a closet. Trim this audio to the clap, and the reverb created by your enclosed space. This is your impulse response.
• Send both sound files to whatever computer you have Maple on.
• Using AudioTools:-Convolution, convolve the dry audio with the impulse response . This your wet audio and should sound a little bit like your voice was recorded in your enclosed space.

Here's some code. I've also attached my dry audio, an impulse response recorded in my shower (yes, I stood inside my shower, closed the door, and recorded a single clap of my hands on my phone), and the resulting wet audio.

with( AudioTools ):
impulse_response := Read( "clap_sc.wav" ):
wet_audio := Normalize( Convolution( dry_audio, impulse_response ) ):
Write("wet_audio.wav", wet_audio );


A full Maple worksheet is here.

AudioSamplesForReverb.zip

## Exploring the CoVid19 outbreak

by: Maple

Hi,

Two weeks ago, I started loading data on the CoVid19 outbreak in order to understand, out of any official communication from any country, what is really going on.

From february 29 to march 9 these data come from https://bnonews.com/index.php/2020/02/the-latest-coronavirus-cases/ and from 10 march until now from https://www.worldometers.info/coronavirus/#repro.In all cases the loading is done manually (copy-paste onto a LibreOffice spreadsheet plus correction and save into a xls file)Â for I wasn't capable to find csv data (csv data do exist here https://github.com/CSSEGISandData/COVID-19, by they end febreuary 15th).
So I copied-pasted the results from the two sources above into a LibreOffice spreadsheet, adjusted the names of some countries for they appeared differently (for instance "United States" instead of "USA"), removed the unnessary commas and saved the result in a xls file.

I also used data from https://www.worldometers.info/world-population/population-by-country/ to get the populations of more than 260 countries around the world and, finally, csv data from https://ourworldindata.org/coronavirus#covid-19-tests to get synthetic histories of confirmed and death cases (I have discovered this site only yesterday evening and I think it could replace all the data I initially loaded).

The two worksheet here are aimed to exploratory and visualization only.
An other one is in progress whose goal is to infer the true death rate (also known as CFR, Case Fatality Rate).

No analysis is presented, if for no other reason than that the available data (except the numbers of deaths) are extremely dependent on the testing policies in place. But some features can be drawn from the data used here.
For instance, if you select country = "China" in file Covid19_Evolution_bis.mw, you will observe very well known behaviour which is that the "Apparent Death Rate", I defined as the ratio of the cumulated number of death at time t by the cumulatibe number of confirmed cases at the same time, is always an underestimation of the death rate one can only known once the outbreak has ended. With this in mind, changing the country in this worksheet from China to Italy seems to lead to frightening  scary interpolations... But here again, without knowing the test policy no solid conclusion can be drawn: maybe Italy tests mainly elder people with accute symptoms, thus the huge "Apparent Death Rate" Italy seems to have?

The work has been done with Maple 2015 and some graphics can be improved if a newer version is used (for instance, as Maple 2015 doesn't allow to change the direction of tickmarks, I overcome this limitation by assigning the date to the vertical axis on some plots).
The second Explore plot could probably be improved by using newer versions or Maplets or Embeded components.

Explore data from https://bnonews.com/index.php/2020/02/the-latest-coronavirus-cases/ and https://www.worldometers.info/coronavirus/#repro
Files to use
Covid19_Evolution.mw
Covid19_Data.m.zip
Population.xls

Explore data from  https://ourworldindata.org/coronavirus#covid-19-tests
Files to use
Covid19_Evolution_bis.mw
daily-deaths-covid-19-who.xls
total-cases-covid-19-who.xls
Population.xls

I would be interested by any open collaboration with people interested by this post (it's not in my intention to write papers on the subject, my only motivation is scientific curiosity).

## Optimizing Mini Golf in MapleSim

Global Optimization , MapleSim 2019

Playing mini-golf recently, I realized that my protractor can only help me so far since it can't calculate the speed of the swing needed.  I decided a more sophisticated tool was needed and modeled a trick-shot in MapleSim.

To start, I laid out the obstacles, the ball and club, the ground, and some additional visualizations in the MapleSim environment.

When running the simulation, my first result wasn't even close to the hole (similar to when I play in real life!).

The model clearly needed to be optimized. I went to the Optimization app in MapleSim (this can be found under Add Apps or Templates  on the left hand side).

Inside the app I clicked "Load System" then selected the parameters I wanted to optimize.

For this case, I'm optimizing 's' (the speed of the club) and 'theta' (the angle of the club). For the Objective Function I added a Relative Translation Sensor to the model and attached a probe to the Vector Norm of the output.

Inside the app, I switched to the Objective Function section.  Selecting Probes, I added the new probe as the Objective Function by giving it a weight of 1.

Scrolling down to "Execute Parameter Optimization", I checked the "Use Global Optimization Toolbox" checkbox, and clicked Run Parameter Optimization.

Following a run time of 120 seconds, the app returns the graph of the objective function.

Below the plot, optimal values for the parameters are given. Plugging these back into the parameter block for the simulation we see that the ball does in fact go into the hole. Success!

Mini_golf_Global_Optimization.msim

## A riddle !

by: Maple 2018

Can you guess what P() produces, without executing it?

P:=proc(N:=infinity) local q,r,t,k,n,l,h, f;
q,r,t,k,n,l,h := 1,0,1,1,3,3,0:
while h<N do
if 4*q+r-t < n*t
then f:=if(++h mod 50=0,"\n",if(h mod 10=0," ","")); printf("%d"||f,n);
q,r,t,k,n,l := 10*q,10*(r-n*t),t,k,iquo(10*(3*q+r),t)-10*n,l
else q,r,t,k,n,l := q*k,(2*q+r)*l,t*l,k+1,iquo(q*(7*k+2)+r*l,t*l),l+2
fi
od: NULL
end:


I hope you will like it (maybe after execution).

## Feynman Diagrams - the scattering matrix in coordi...

by: Maple

Feynman Diagrams
The scattering matrix in coordinates and momentum representation

 Mathematical methods for particle physics was one of the weak spots in the Physics package. There existed a FeynmanDiagrams command, but its capabilities were too minimal. People working in the area asked for more functionality. These diagrams are the cornerstone of calculations in particle physics (collisions involving from the electron to the Higgs boson), for example at the CERN. As an introduction for people curious, not working in the area, see "Why Feynman Diagrams are so important".
 This post is thus about a new development in Physics: a full rewriting of the FeynmanDiagrams command, now including a myriad of new capabilities (mainly a. b. and c. in the Introduction), reversing the previous status of things entirely. This is work in collaboration with Davide Polvara from Durham University, Centre for Particle Theory.
 The complexity of this material is high, so the introduction to the presentation below is as brief as it can get, emphasizing the examples instead. This material is reproducible in Maple 2019.2 after installing the Physics Updates, v.598 or higher.
 At the end they are attached the worksheet corresponding to this presentation and a PDF version of it, as well as the new FeynmanDiagrams help page with all the explanatory details.

Introduction

 A scattering matrix  relates the initial and final states,  and , of an interacting system. In an 4-dimensional spacetime with coordinates ,  can be written as:

 where  is the imaginary unit  and  is the interaction Lagrangian, written in terms of quantum fields  depending on the spacetime coordinates  . The T symbol means time-ordered. For the terminology used in this page, see for instance chapter IV, "The Scattering Matrix", of ref.[1] Bogoliubov, N.N., and Shirkov, D.V. Quantum Fields.
 This exponential can be expanded as

 where

 and  is the time-ordered product of n interaction Lagrangians evaluated at different points. The S matrix formulation is at the core of perturbative approaches in relativistic Quantum Field Theory.
 In connection, the FeynmanDiagrams  command has been rewritten entirely for Maple 2020. In brief, the new functionality includes computing:
 a. The expansion  in coordinates representation up to arbitrary order (the limitation is now only your hardware)
 b. The S-matrix element  in momentum representation up to arbitrary order for given number of loops and initial and final particles (the contents of the  and  states); optionally, also the transition probability density, constructed using the square of the scattering matrix element , as shown in formula (13) of sec. 21.1 of ref.[1].
 c. The Feynman diagrams (drawings) related to the different terms of the expansion of S or of its matrix elements .
 Interaction Lagrangians involving derivatives of fields, typically appearing in non-Abelian gauge theories, are also handled, and several options are provided enabling restricting the outcome in different ways, regarding the incoming and outgoing particles, the number of loops, vertices or external legs, the propagators and normal products, or whether to compute tadpoles and 1-particle reducible terms.

Examples

For illustration purposes set three coordinate systems , and set  to represent a quantum operator

 >
 >
 (1.1)

Let  be the interaction Lagrangian

 >
 (1.2)

The expansion of  in coordinates representation, computed by default up to order = 3 (you can change that using the option order = n), by definition containing all possible configurations of external legs, displaying the related Feynman Diagrams, is given by

 >

 (1.3)

The expansion of   in coordinates representation to a specific order shows in a compact way the topology of the underlying Feynman diagrams. Each integral is represented with a new command, FeynmanIntegral , that works both in coordinates and momentum representation. To each term of the integrands corresponds a diagram, and the correspondence is always clear from the symmetry factors.

In a typical situation, one wants to compute a specific term, or scattering process, instead of the S matrix up to some order with all possible configurations of external legs. For example, to compute only the terms of this result that correspond to diagrams with 1 loop use numberofloops = 1 (for tree-level, use numerofloops = 0)

 >

 (1.4)

In the result above there are two terms, with 4 and 6 external legs respectively.

A scattering process with matrix element  in momentum representation, corresponding to the term with 4 external legs (symmetry factor = 72), could be any process where the total number of incoming + outgoing parties is equal to 4. For example, one with 2 incoming and 2 outgoing particles. The transition probability for that process is given by

 >

 (1.5)

When computing in momentum representation, only the topology of the corresponding Feynman diagrams is shown (i.e. the diagrams associated to the corresponding Feynman integral in coordinates representation).

The transition matrix element  is related to the transition probability density  (formula (13) of sec. 21.1 of ref.[1]) by

where  represent the particle densities of each of the s particles in the initial state , the  (Dirac) is the expected singular factor due to the conservation of the energy-momentum and the amplitude is related to  via

To directly get the probability density  instead ofuse the option output = probabilitydensity

 >
 (1.6)

In practice, the most common computations involve processes with 2 or 4 external legs. To restrict the expansion of the scattering matrix in coordinates representation to that kind of processes use the numberofexternallegs option. For example, the following computes the expansion of  up to order = 3, restricting the outcome to the terms corresponding to diagrams with only 2 external legs

 >

 (1.7)

This result shows two Feynman integrals, with respectively 2 and 3 loops, the second integral with two terms. The transition probability density in momentum representation for a process related to the first integral (1 term with symmetry factor = 96) is then

 >

 (1.8)

In the above, for readability, the contracted spacetime indices in the square of momenta entering the amplitude F (as denominators of propagators) are implicit. To make those indices explicit, use the option putindicesinsquareofmomentum

 >
 (1.9)

This computation can also be performed to higher orders. For example, with 3 loops, in coordinates and momentum representations, corresponding to the other two terms and diagrams in (1.7)

 >
 (1.10)

A corresponding S-matrix element in momentum representation:

 >
 (1.11)

Consider the interaction Lagrangian of Quantum Electrodynamics (QED). To formulate this problem on the worksheet, start defining the vector field .

 >
 (1.12)

Set lowercase Latin letters from i to s to represent spinor indices (you can change this setting according to your preference, see Setup ), also the (anticommutative) spinor field will be represented by , so set  as an anticommutativeprefix, and set  and  as quantum operators

 >
 (1.13)

The matrix indices of the Dirac matrices  are written explicitly and use conjugate  to represent the Dirac conjugate

 >
 (1.14)

Compute , only the terms with 4 external legs, and display the diagrams: all the corresponding graphs have no loops

 >

 (1.15)

The same computation but with only 2 external legs results in the diagrams with 1 loop that correspond to the self-energy of the electron and the photon (page 218 of ref.[1])

 >

 (1.16)

where the diagram with two spinor legs is the electron self-energy. To restrict the output furthermore, for example getting only the self-energy of the photon, you can specify the normal products you want:

 >
 (1.17)

The corresponding S-matrix elements in momentum representation

 >

 (1.18)

In this result we see  spinor (see ref.[2]), and the propagator of the field  with a mass . To indicate that this field is massless use

 >
 (1.19)

Now the propagator for  is the one of a massless vector field:

 >
 (1.20)

The self-energy of the photon:

 >
 (1.21)

where  is the corresponding polarization vector.

When working with non-Abelian gauge fields, the interaction Lagrangian involves derivatives. FeynmanDiagrams  can handle that kind of interaction in momentum representation. Consider for instance a Yang-Mills theory with a massless field  where a is a SU2 index (see eq.(12) of sec. 19.4 of ref.[1]). The interaction Lagrangian can be entered as follows

 >
 (1.22)
 >
 >
 (1.23)
 >
 (1.24)

The transition probability density at tree-level for a process with two incoming and two outgoing B particles is given by

 >
 (1.25)

 (1.26)

To simplify the repeated indices, us the option simplifytensorindices. To check the indices entering a result like this one use Check ; there are no free indices, and regarding the repeated indices:

 >
 (1.27)

This process can be computed with 1 or more loops, in which case the number of terms increases significantly. As another interesting non-Abelian model, consider the interaction Lagrangian of the electro-weak part of the Standard Model

 >
 (1.28)
 >
 (1.29)
 >
 (1.30)
 >
 (1.31)
 >
 (1.32)
 >
 (1.33)
 >
 (1.34)

This interaction Lagrangian contains six different terms. The S-matrix element for the tree-level process with two incoming and two outgoing W particles is shown in the help page for FeynmanDiagrams .

 >
 References [1] Bogoliubov, N.N., and Shirkov, D.V. Quantum Fields. Benjamin Cummings, 1982. [2] Weinberg, S., The Quantum Theory Of Fields. Cambridge University Press, 2005.

FeynmanDiagrams_and_the_Scattering_Matrix.PDF

FeynmanDiagrams_and_the_Scattering_Matrix.mw

FeynmanDiagrams_-_help_page.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

## Happy New Year!

by: Maple
with(plots):
 >

## Putnam Mathematical Competition 2019

by: Maple 2018

Maple can easily solve the B4 problem of the Putnam Mathematical Competition 2019  link

B4.  Let F be the set of functions f(x,y) that are twice continuously differentiable for x≥1, y≥1 and that satisfy the following two equations:

For each f2F, let

Determine m(f), and show that it is independent of the choice of f.

 > # Solution
 > pdsolve({ x*diff(f(x,y),x)+y*diff(f(x,y),y) = x*y*ln(x*y), x^2*diff(f(x,y),x,x)+y^2*diff(f(x,y),y,y) = x*y });
 (1)
 > f:=unapply(rhs(%[]), x,y);
 (2)
 > h := f(s+1, s+1) - f(s+1, s) - f(s, s+1) + f(s, s);
 (3)
 > minimize(h, s=1..infinity);
 (4)
 > answer = simplify(%);
 (5)
 >

## Triangulation, hatching and texturing of simple...

by: Maple

Hi,

As an amusement,  I decided several months ago to develop some procedures to fill a simple polygon* by hatches or simple textures.

* A simple polygon is a polygon  whose sides either do not intersect or either have a common vertex.

This work started with the very simple observation that if I was capable to hatch or texture a triangle, I will be too to hatch or texture any simple polygon once triangulated.
I also did some work to extend this work to non-simple polygons but there remains some issues to fix (which explains while it is not deliverd here).

The main ideat I used for hatching and texturing is based upon the description of each triangles by a set of 3 inequalities that any interior point must verify.
A hatch of this triangle is thius a segment whose each point is interior.
The closely related idea is used for texturing. Given a simple polygon, periodically replicated to form the texture, the set of points of each replicate that are interior to a given triangle must verify a set of inequalities (the 3 that describe the triangle, plus N if the pattern of the texture is a simple polygon with N sides).

Unfortunately I never finalise this work.
Recently @Christian Wolinski asked a question about texturing that reminded me this "ancient" work of mine.
So I decided to post it as it is, programatically imperfect, lengthy to run, and most of all french-written for a large part.
I guess it is a quite unorthodox way to proceed but some here could be interested by this work to take it back and improve it.

The module named "trianguler" (= triangulate) is a translation into Maple of Frederic Legrand's Python code (full reference given in the worksheet).
I added my own procedure "hachurer" (= hatching) to this module.
The texturing part is not included in this module for it was still in development.

A lot of improvements can be done that I could list, but I prefer not to be too intrusive in your evaluation of this work. So make your own idea about it and do not hesitate to ask me any informations you might need (in particular translation questions).

PS: this work has been done with Maple 2015.2

 > restart:

Reference: http://www.f-legrand.fr/scidoc/docmml/graphie/geometrie/polygone/polygone.html
(in french)
reference herein : M. de Berg, O. Cheong, M. van Kreveld, M. Overmars,
Computational geometry,  (Springer, 2010)

Direct translation of the Frederic Legrand's Python code

Meaning of the different french terms

voisin_sommet  (n, i, di)
let L the list [1, ..., n] where n is the number of vertices
This procedure returns the index of the neighbour (voisin) of the vertex (sommet) i when L is rotated by di

equation_droite  (P0, P1, M)
Let P0 and P1 two vertices and M an arbitrary point.
Let (P0, P1) the vector originated at P0 and ending at P1 (idem for (P0, M)) and  the unitary vector in the Z direction.
This procedure returns (P0, P1) o (P0, M) .

point_dans_triangle  (triangle, M) P1, P2]
This procedure returns "true" if point M is within (strictly) the  triangle "triangle" and "false" if not.

sommet_distance_maximale  (polygone, P0, P1, P2, indices)
Given a polygon (polygone) threes vertices P0, P1 and P2 and a list of indices , this procedure returns
the vertex of the polygon "polygone" which verifies: 1/ this vertex is strictly within
the triangle [P0, P1, P2] and 2/ it is the farthest from side [P1, P2] amid all the vertices that verifies point 1/.
If there is no such vertex the procedure returns NULL.

sommet_gauche (polygone)
This procedure returns the index of the leftmost ("gauche" means "left") vertex in polygon "polygone".
If more than one vertices have the same minimum abscissa value then only the first one is returned.

nouveau_polygone(polygone,i_debut,i_fin)
This procedure creates a new polygon from index i_debut (could be translated by i_first) to i_end (i_last)

trianguler_polygone_recursif(polygone)
This procedure recursively divides a polygon in two parts A and B from its leftmost vertex.
If A (B) is a triangle the list "liste_triangles" (mening "list of triangles") is augmented by A (B);
if not the procedure executes recursively on A and B.

trianguler_polygone(polygone)
This procedure triangulates the polygon "polygon"

hachurer(shapes, hatch_angle, hatch_number, hatch_color)
This procedure generates stes of hatches of different angles, colors and numbers

Limitations:
1/ "polygone" is a simply connected polygon
2/  two different sides S and S', either do not intersect or either have a common vertex

 > trianguler := module() export voisin_sommet, equation_droite, interieur_forme, point_dans_triangle, sommet_distance_maximale,        sommet_gauche, nouveau_polygone, trianguler_polygone_recursif, trianguler_polygone, hachurer: #------------------------------------------------------------------- voisin_sommet := (n, i, di) -> ListTools:-Rotate([$1..n], di)[i]: #------------------------------------------------------------------- equation_droite := proc(P0, P1, M) (P1[1]-P0[1])*(M[2]-P0[2]) - (P1[2]-P0[2])*(M[1]-P0[1]) end proc: #------------------------------------------------------------------- interieur_forme := proc(forme, M) local N: N := numelems(forme); { seq( equation_droite(forme[n], forme[piecewise(n=N, 1, n+1)], M) >= 0, n=1..N) } end proc: #------------------------------------------------------------------- point_dans_triangle := proc(triangle, M) and( is( equation_droite(triangle[1], triangle[2], M) > 0 ), is( equation_droite(triangle[2], triangle[3], M) > 0 ), is( equation_droite(triangle[3], triangle[1], M) > 0 ) ) end proc: #------------------------------------------------------------------- sommet_distance_maximale := proc(polygone, P0, P1, P2, indices) local n, distance, j, i, M, d; n := numelems(polygone): distance := 0: j := NULL: for i from 1 to n do if not(member(i, indices)) then M := polygone[i]; if point_dans_triangle([P0, P1, P2], M) then d := abs(equation_droite(P1, P2, M)): if d > distance then distance := d: j := i end if: end if: end if: end do: return j: end proc: #------------------------------------------------------------------- sommet_gauche := polygone -> sort(polygone, key=(x->x[1]), output=permutation)[1]: #------------------------------------------------------------------- nouveau_polygone := proc(polygone, i_debut, i_fin) local n, p, i: n := numelems(polygone): p := NULL: i := i_debut: while i <> i_fin do p := p, polygone[i]: i := voisin_sommet(n, i, 1) end do: p := [p, polygone[i_fin]] end proc: #------------------------------------------------------------------- trianguler_polygone_recursif := proc(polygone) local n, j0, j1, j2, P0, P1, P2, j, polygone_1, polygone_2: global liste_triangles: n := numelems(polygone): j0 := sommet_gauche(polygone): j1 := voisin_sommet(n, j0, +1): j2 := voisin_sommet(n, j0, -1): P0 := polygone[j0]: P1 := polygone[j1]: P2 := polygone[j2]: j := sommet_distance_maximale(polygone, P0, P1, P2, [j0, j1, j2]): if not(j::posint) then liste_triangles := liste_triangles, [P0, P1, P2]: polygone_1 := nouveau_polygone(polygone,j1,j2): if numelems(polygone_1) = 3 then liste_triangles := liste_triangles, polygone_1: else thisproc(polygone_1) end if: else polygone_1 := nouveau_polygone(polygone, j0, j ): polygone_2 := nouveau_polygone(polygone, j , j0): if numelems(polygone_1) = 3 then liste_triangles := liste_triangles, polygone_1: else thisproc(polygone_1) end if: if numelems(polygone_2) = 3 then liste_triangles := liste_triangles, polygone_2: else thisproc(polygone_2) end if: end if: return [liste_triangles]: end proc: #------------------------------------------------------------------- trianguler_polygone := proc(polygone) trianguler_polygone_recursif(polygone): return liste_triangles: end proc: #------------------------------------------------------------------- hachurer := proc(shapes, hatch_angle::list, hatch_number::list, hatch_color::list) local A, La, Lp; local N, P, _sides, L_sides, Xshape, ch, rel, p_rel, n, sol, p_range: local AllHatches, window, p, _segment: local NT, ka, N_Hatches, p_range_t, nt, shape, p_hatches, WhatHatches: #----------------------------------------------------------------- # internal functions: # # La(x, y, alpha, p) is the implicit equation of a straight line of angle alpha relatively # to the horizontal axis and intercept p # # Lp(x, y, P) is the implicit equation of a straight line passing through points P[1] and P[2] # # interieur_triangle(triangle, M) La := (x, y, alpha, p) -> cos(alpha)*x - sin(alpha)*y + p; Lp := proc(x, y, P::list) (x-P[1][1])*(P[2][2]-P[1][2]) - (y-P[1][2])*(P[2][1]-P[1][1] = 0) end proc; p_range := [+infinity, -infinity]: NT := numelems(shapes): AllHatches := NULL: for ka from 1 to numelems(hatch_angle) do A := hatch_angle[ka]: N_Hatches := hatch_number[ka]: p_range_t := NULL: _sides := []: L_sides := []: rel := []: for nt from 1 to NT do shape := shapes[nt]: # _sides : two points description of the sides of the shape # L_sides : implicit equations of the straight lines that support the sides N := [1, 2, 3]; P := [2, 3, 1]; _sides := [ _sides[] , [ seq([shape[n], shape[P[n]]], n in N) ] ]; L_sides := [ L_sides[], Lp~(x, y, _sides[-1]) ]; # Inequalities that define the interior of the shape rel := [ rel[], trianguler:-interieur_forme(shape, [x, y]) ]; # Given the orientation of the hatches we search here the extreme values of # the intercept p for wich a straight line of equation La(x, y, alpha, p) # cuts the shape. p_rel := NULL: for n from 1 to numelems(L_sides[-1]) do sol := solve({La(x, y, A, q), lhs(L_sides[-1][n])} union rel[-1], [x, y]); p_rel := p_rel, if(sol <> [], [rhs(op(1, %)), rhs(op(3, %))], [+infinity, -infinity]); end do: p_range_t := p_range_t, evalf(min(op~(1, [p_rel]))..max(op~(2, [p_rel]))); p_range := evalf(min(op~(1, [p_rel]), op(1, p_range))..max(op~(2, [p_rel]), op(2, p_range))); end do: # end of the loop over triangles p_range_t := [p_range_t]: p_hatches := [seq(p_range, (op(2, p_range)-op(1, p_range))/N_Hatches)]: # Building of the hatches # # This construction is far from being optimal. # Here again the main goal was to obtain the hatches with a minimum effort # if algorithmic development. window := min(op~(1..shape))..max(op~(1..shape)): WhatHatches := map(v -> map(u -> if verify(u, v, 'interval'('closed') ) then u end if, p_hatches), p_range_t): for nt from 1 to NT do for p in WhatHatches[nt] do _segment := []: for n from 1 to numelems(L_sides[nt]) do _segment := _segment, evalf( solve({La(x, y, A, p), lhs(L_sides[nt][n])} union rel[nt], [x, y]) ); end do; map(u -> u[], [_segment]); AllHatches := AllHatches, plot(map(u -> rhs~(u), %), color=hatch_color[ka]): end do: end do; end do: # end of loop over hatch angles plots:-display( PLOT(POLYGONS(polygone, COLOR(RGB, 1$3))),   AllHatches,   scaling=constrained ) end proc: end module:

Legrand's example (see reference above)

 > global liste_triangles: liste_triangles := NULL:
 > polygone := [[0,0],[0.5,-1],[1.5,-0.2],[2,-0.5],[2,0],[1.5,1],[0.3,0],[0.5,1]]: trianguler:-trianguler_polygone(polygone): PLOT(seq(POLYGONS(u, COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12)), u in liste_triangles), VIEW(0..2, -2..2))
 > trianguler:-hachurer([liste_triangles], [-Pi/4, Pi/4], [40, 40], [red, blue])
 > F := (P, a, b) -> map(p -> [p[1]+a, p[2]+b], P):
 > MOTIF  := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]]; motifs := [ seq(seq(F(MOTIF, 0+i*0.075, 0+j*0.075), i=0..26), j=-14..13) ]: plots:-display(   plot([polygone[], polygone[1]], color=red, filled),   map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs) ): texture    := NULL: rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):    for ref in liste_triangles do   ref;   #   # the three lines below are used to define REF counter clockwise   #   g           := trianguler:-sommet_gauche(ref):   bas         := sort(op~(2, ref), output=permutation);   REF         := ref[[g, op(map(u -> if u<>g then u end if, bas))]];   rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);   texture_ref := map(u -> plots:-inequal(rel_ref union u, x=0..2, y=-1..1, color=blue, 'nolines'), rel_motifs):   texture     := texture, texture_ref: end do: plots:-display(   plot([polygone[], polygone[1]], color=red, scaling=constrained),   texture )
 > MOTIF  := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]]; motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.05, 0.1)+i*0.1, 0+j*0.05), i=-0.2..20), j=-20..20) ]: plots:-display(   plot([polygone[], polygone[1]], color=red, filled),   map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs) ): texture    := NULL: rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):    for ref in liste_triangles do   ref;   g := trianguler:-sommet_gauche(ref):   bas := sort(op~(2, ref), output=permutation);   REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];   rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);   texture_ref := map(u -> plots:-inequal(rel_ref union u, x=0..2, y=-1..1, color=blue, 'nolines'), rel_motifs):   texture     := texture, texture_ref: end do: plots:-display(   plot([polygone[], polygone[1]], color=red, scaling=constrained),   texture )
 > MOTIF  := [[0, 0], [0.4, 0], [0.4, 0.14], [0, 0.14]]: motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.4, 0.2)+i*0.4, 0+j*0.14), i=-1..4), j=-8..7) ]: plots:-display(   plot([polygone[], polygone[1]], color=red, filled),   map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs) ): palettes := ColorTools:-PaletteNames(): ColorTools:-GetPalette("HTML"): couleurs := [SandyBrown, PeachPuff, Peru, Linen, Bisque, Burlywood, Tan, AntiqueWhite,      NavajoWhite, BlanchedAlmond, PapayaWhip, Moccasin, Wheat]: nc   := numelems(couleurs): roll := rand(1..nc): motifs_nb      := numelems(motifs): motifs_couleur := [ seq(cat("HTML ", couleurs[roll()]), n=1..motifs_nb) ]: texture    := NULL: rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):    for ref in liste_triangles do   ref;   g := trianguler:-sommet_gauche(ref):   bas := sort(op~(2, ref), output=permutation);   REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];   rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);   texture_ref := map(n -> plots:-inequal(rel_ref union rel_motifs[n], x=0..2, y=-1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]): texture := texture, texture_ref: end do: plots:-display( plot([polygone[], polygone[1]], color=red, scaling=constrained), texture )  > MOTIF := [ seq(0.1*~[cos(Pi/6+Pi/3*i), sin(Pi/6+Pi/3*i)], i=0..5) ]: motifs := [ seq(seq(F(MOTIF, i*0.2*cos(Pi/6)+piecewise(j::odd, 0, 0.08), j*0.3*sin(Pi/6)), i=0..12), j=-6..6) ]: plots:-display( plot([polygone[], polygone[1]], color=red, filled), map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs) ): motifs_nb := numelems(motifs): motifs_couleur := [ seq(if(n::even, yellow, brown) , n=1..motifs_nb) ]: texture := NULL: rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs): for ref in liste_triangles do ref; g := trianguler:-sommet_gauche(ref): bas := sort(op~(2, ref), output=permutation); REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]]; rel_ref := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref); texture_ref := map(n -> plots:-inequal(rel_ref union rel_motifs[n], x=0..2, y=-1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):   texture     := texture, texture_ref: end do: plots:-display(   plot([polygone[], polygone[1]], color=red, scaling=constrained),   texture )
 >