diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro
index a5b4f33089c0e6dd5fa7116243a88d4b4a06e573..91d96e44d411b2c77b8caf1c23d1462a2e4ed5b3 100644
--- a/Electrostatics/microstrip.pro
+++ b/Electrostatics/microstrip.pro
@@ -36,12 +36,12 @@
    boundary condition (zero flux of the displacement field, i.e. n.d = 0) is
    imposed on the left boundary of the domain to account for the symmetry of the
    problem, as well as on the top and right boundaries that truncate the
-   modelling domain. */
+   simulation domain. */
 
 Group {
   /* One starts by giving explicit meaningful names to the Physical regions
-     defined in the "microstrip.msh" mesh file. We only use 2 volume regions and
-     2 surface regions in this model. */
+     defined in the "microstrip.msh" mesh file. This model comprises only 
+     2 volume regions and 2 surface regions. */
 
   Air = Region[101];
   Diel1 = Region[111];
@@ -61,6 +61,10 @@ Group {
 
      Since there are no non-homogeneous Neumann conditions in this particular
      example, Sur_Neu_Ele is defined as empty.
+     
+     Note that volume elements are those that correspond to the higher dimension 
+     of the model at hand (2D elements here), surface elements correspond to the 
+     higher dimension of the model minus one (1D elements here).
      */
 
   Vol_Ele = Region[ {Air, Diel1} ];
@@ -78,8 +82,8 @@ Function {
 }
 
 Constraint {
-  /* As for material laws, the Dirichlet boundary condition is defined
-     piecewise.  The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
+  /* The Dirichlet boundary condition is also defined piecewise.  
+     The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
      below. */
 
   { Name Dirichlet_Ele; Type Assign;
@@ -97,13 +101,13 @@ Group{
      Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or
      surface regions, resp., domains of definition Dom_xxx regions may contain
      both volume and surface regions. Hence the use of the prefixes Vol_, Sur_
-     and Dom_ to avoid confusions.*/
+     and Dom_ to avoid confusion.*/
 
   Dom_Hgrad_v_Ele =  Region[ {Vol_Ele, Sur_Neu_Ele} ];
 }
 
 FunctionSpace {
-  /* The function space in which we shall pick the electric scalar potential "v"
+  /* The function space in which we pick the electric scalar potential "v"
      solution is defined by
      - a domain of definition (the "Support": "Dom_Hgrad_v_Ele")
      - a type ("Form0" means scalar field)
@@ -127,7 +131,7 @@ FunctionSpace {
       { Name sn; NameOfCoef vn; Function BF_Node;
         Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; }
       // using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an
-      // optimization, which allows GetDP to not explicitly build the list of
+      // optimization, which avoids explicitly building the list of 
       // all the nodes
     }
     Constraint {
@@ -138,7 +142,7 @@ FunctionSpace {
 
 Jacobian {
   /* Jacobians are used to specify the mapping between elements in the mesh and
-     the reference elements (defined in standardized unit cells) over which
+     the reference elements (defined in standardized/reference unit cells) over which
      integration is performed. "Vol" represents the classical 1-to-1 mapping
      between identical spatial dimensions, i.e. in this case a reference
      triangle/quadrangle onto triangles/quadrangles in the z=0 plane (2D <->
@@ -173,7 +177,7 @@ Integration {
 }
 
 Formulation {
-  /* The syntax of the Formulation section is a harder nut to crack. So let's
+  /* The syntax of the Formulation section is a harder nut to crack. So let us
      deal with it carefully.
 
      A GetDP Formulation encodes a so-called weak formulation of the original
@@ -201,11 +205,11 @@ Formulation {
      ones used to interpolate the unknown potential v.
 
      The "Integral" statement in the Formulation is a symbolic representation of
-     this weak formulation. It has got 4 semicolon separated arguments:
+     this weak formulation. It has 4 semicolon separated arguments:
      * the density [.,.] to be integrated (note the square brackets instead of
-       the parentheses), with the test-functions (always) after the comma
-     * the domain of integration,
-     * the Jacobian of the transformation reference element <-> real element,
+       the parentheses), with the test-functions (always) after the comma;
+     * the domain of integration;
+     * the Jacobian of the transformation reference element <-> real element;
      * the integration method to be used.
 
      In the density, braces around a quantity (such as "{v}") indicate that this
@@ -223,33 +227,33 @@ Formulation {
      be something like this [ ... , basis_functions_of {d v} ]. However, since
      the second term is always devoted to test functions, the operator
      "basis_functions_of" would always be there.  It can therefore be made
-     implicit, and, according to the GetDP syntax, it is omitted. So, one writes
-     simply [ ... , {d v} ].
+     implicit and, according to the GetDP syntax, omitted. So, one simply writes
+     [ ... , {d v} ].
 
-     The first term, on the other hand, can contain a much wider variety of
+     On the other hand, the first term can contain a much wider variety of
      expressions than the second one. In our case it should be expressed in
      terms of the FE expansion of "v" at the present system solution, i.e. when
      the coefficients vn_k in the expansion of "v = Sum_k vn_k sn_k" are
-     unknown. This is indicated by prefixing the braces with "Dof", which leads
-     to the following density:
+     unknown. This is indicated by prefixing the braces with "Dof" (degrees of 
+     freedom), which leads to the following density:
 
      [ epsilon[] * Dof{d v} , {d v} ],
 
      a so-called bilinear term that will contribute to the stiffness matrix of
      the electrostatic problem at hand.
 
-     Another option, which would not work here, is to evaluate the first
-     argument with the last available already computed solution, i.e. simply
-     perform the interpolation with known coefficients vn_k. For this the Dof
-     prefix would be omitted and we would have:
+     Another option, which does not work here, is to evaluate the first
+     argument with the last available computed solution, i.e. simply
+     perform the interpolation with known coefficients vn_k. The Dof
+     prefix is then omitted and we have:
 
      [ epsilon[] * {d v} , {d v} ],
 
-     a so-called linear term that will contribute to the right-hand side of the
+     a so-called linear term that contributes to the right-hand side of the
      linear system.
 
-     Both choices are commonly used in different contexts, and we shall come
-     back on this often in subsequent tutorials. */
+     Both choices are commonly used in different contexts, and we shall often 
+     come back in subsequent tutorials. */
 
   { Name Electrostatics_v; Type FemEquation;
     Quantity {
@@ -259,20 +263,21 @@ Formulation {
       Integral { [ epsilon[] * Dof{d v} , {d v} ];
 	    In Vol_Ele; Jacobian Vol; Integration Int; }
 
-   /* Additional Integral terms could be added here. For example, the
-      following term would account for non-homogeneous Neumann boundary
+   /* Additional Integral terms can be added here. For example, the
+      following term may account for non-homogeneous Neumann boundary
       conditions, provided that the function nd[] is defined:
 
       Integral { [ nd[] , {v} ];
 	    In Sur_Neu_Ele; Jacobian Sur; Integration Int; }
 
-      All the terms are added, and an implicit "= 0" is considered at the end. */
+      All the terms in the Equation environment are added, 
+      and an implicit "= 0" is considered at the end. */
     }
   }
 }
 
-/* What to do with a weak formulation is specified in the Resolution: here we
-   simply generale a linear system, solve it and save the solution (.res file)
+/* In the Resolution environment we specify what to do with a weak formulation: 
+   here we simply generate a linear system, solve it and save the solution (.res file)
    to disk. */
 
 Resolution {
@@ -291,18 +296,18 @@ Resolution {
    The first part defines, in terms of the Formulation, which itself refers to
    the FunctionSpace, a number of quantities that can be evaluated at the
    postprocessing stage. The three quantities defined here are:
-   - the scalar electric potential,
-   - the electric field,
+   - the scalar electric potential;
+   - the electric field;
    - the electric displacement.
 
    The second part consists in defining groups of post-processing operations,
    which can be invoked separately. The first group is invoked by default when
    Gmsh is run interactively. Each Operation specifies
-   - a quantity to be eveluated,
-   - the region on which the evaluation is done,
+   - a quantity to evaluate;
+   - the region on which the evaluation is done;
    - the name of the output file.
    The generated post-processing files are automatically displayed by Gmsh if
-   the "Merge result automatically" option is enabled (which is the default). */
+   the "Merge result automatically" option is enabled (the default). */
 
 PostProcessing {
   { Name EleSta_v; NameOfFormulation Electrostatics_v;