diff --git a/vignettes/SimulatingMicrostructures.Rmd b/vignettes/SimulatingMicrostructures.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..5e6bb56e2811a7d80174f48b151d2859a50427a1
--- /dev/null
+++ b/vignettes/SimulatingMicrostructures.Rmd
@@ -0,0 +1,342 @@
+---
+title: "Simulating microstructures with gmGeostats"
+output: rmarkdown::html_vignette
+vignette: >
+ %\VignetteIndexEntry{my-vignette}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ echo = TRUE,
+ comment = "#>",
+ fig.width = 7,
+ fig.height = 6
+)
+
+# ToDo: anistropoy: what means the ratio exactly?
+```
+
+```{r setup, echo=FALSE, message=FALSE}
+library(gmGeostats)
+library(ggplot2)
+library(data.table) # only needed, because Solveig is faster with melting in data.table than with dplyr and Co. can be replaced by dplyr, if someone changes the code from data.table to dplyr later.
+```
+
+# Part 1: Granite/Gneiss
+
+Granites or Gneiss consist mainly of Quartz, Feldspars and micas, e.g. Biotite and/or Muscovite.
+Due to the brittle characteristics these types of rocks have often fractures or joint systems.
+With `gmGeostats` mineral composition, anisotrpy of certain mineral phases and features like joints and fractures can be simulated, e.g. if simulations of crystalline rocks are required for testing code or work flows.
+
+## General idea - Turning bands algorithm
+
+here comes some nice background from Raimon
+
+## Example for simulation
+
+The simulation consists of five steps:
+
+1. decide for the number of phases to be simulated. Fracture and joint systems count as a "phase", e.g. if there is one joint system, this would be one phase, if there are two joint systems this would be simulated by two additional phases.
+2. define the characteristics of each phase, i.e. the crystal elongation (or void shape) and orientation/foliation
+3. generate the variograms by using a function to define for each phase the model of the phase-variogram, the nugget, the sill, the range and the anistropy with respect to the range.
+4. use the algorithm of Turning Bands to calculate the probability for each phase to occur.
+5. Visualize the result.
+
+In the following example, we will show the implementation for 2D.
+Because the code expects 3D the third dimension is always "muted" by using zeros.
+
+### Number of phases
+
+For the example let's assume five phases, Quartz (qu), K-Feldspar (fsp), Biotite (bi), a fracture system (frac) and a joints system (joints).
+
+Generate a template matrix with only zeros, which has the dimensions `n x n`:
+```{r phases}
+n_phases = 5
+
+template = matrix(0, nrow = n_phases, ncol = n_phases)
+phases = c("qu", "fsp", "bi", "frac", "joints")
+colnames(template) <- rownames(template) <- phases
+```
+In a later step when the five variograms are stacked, these five matrices define the sill for each phase.
+Therefore, on the diagonal we replace the `0` with a `1` for the respective phase:
+
+```{r}
+s1 = template
+s1[1,1] <-1
+s2 = template
+s2[2,2] <-1
+s3 = template
+s3[3,3] <-1
+s4 = template
+s4[4,4] <-1
+s5 = template
+s5[5,5] <-1
+```
+
+The example for Quartz looks like this:
+```{r}
+s1
+```
+
+### Characteristics of phases
+
+The elongation and orientation of the elongation for each phase can be simulated by the function `anis2D.par2A`.
+The function parameter `ratio` provides the ratio of ??? (length to width?).
+E.g. `ratio = 1` is a sphere (isotrop), while `ratio = 2` would describe a mineral which is ??? (double as long as thick?).
+The function parameter `angle` gives the orientation (???, does not really, the orientation seems to be somewhat different) of this anisotropy, counted counterclockwise from the east (E).
+E.g. `angle = 45` with `ratio = 2` would result into North-South elongated crystals.
+
+For Quartz, let's assume we have mostly roundish/quadratic crystals and only a minor preferred orientation of the crystals, so that they are slightly oriented along NNE-SSW:
+```{r p-qu}
+p1_qu = anis2D.par2A(ratio = .9, angle = 45)
+```
+
+The Feldspars have a more pronounced elongated shape, also with preferred orientation into the same direction as the Quartz grains:
+
+```{r p-fsp}
+p2_fsp = anis2D.par2A(ratio = .6, angle = 45)
+```
+
+If the anistropy of the mineral phases had been generated by the function `anis2D.par2A`, then we need to calculate the cross product of the matrices to gain the anistropy matrices `A`:
+
+```{r}
+A1_qu = tcrossprod(p1_qu)
+A2_fsp = tcrossprod(p2_fsp)
+```
+
+This results into
+```{r}
+A1_qu
+```
+
+for Quartz and into
+```{r}
+A2_fsp
+```
+
+for Feldspar.
+
+There is also the option to directly generate the anisotropy matrix.
+In this case, the matrix can be written manually:
+The Biotite occurs as thin plates and has a slightly different preferred orientation then the two main minerals:
+```{r}
+A3_bi = matrix(data = c(1, .6, 0, .6, 1, 0, 0, 0, 1), ncol = 3)
+A3_bi
+```
+
+For the fracture and joint systems the anisotropy is very pronounced, because they are very thin but passing through longer distances through the rock.
+Accordingly, the anisotropy matrix should be chosen similar to this:
+```{r}
+A4_frac = matrix(data = c(1, -.999, 0, -.999, 1, 0, 0, 0, 1), ncol = 3)
+A4_frac
+```
+
+And let's assume, that the joint system is orthogonal to the fracture system:
+```{r}
+A5_joints = matrix(data = c(1, .999, 0, .999, 1, 0, 0, 0, 1), ncol = 3)
+A5_joints
+```
+
+### Generate variograms
+
+Fore each phase, variogram functions will predict the probability to occur.
+The function `setCgram` generates these variograms.
+It has four parameters:
+
+- The `type` provides the variogram model, e.g. Gaussian or Exponential.
+Gaussian models provide smooth transitions to other minerals, and may be used if bands of minerals should be modelled.
+Exponential models provide rather abrupt changes between mineral phases and should be preferred model for sharp edges.
+In this example, Quarts and Feldspar are modelled by Gaussion types, and Biotite and the fractures by Exponential type.
+- The `nugget` gives the nugget effect, but in the case of microstructures it can be omitted.
+- For the `sill`, the partial sill of the correlation function, the microstructure require simply the matrices with one `1` in the diagonal.
+We produced these matrices in the first part, where the number of phases was defined.
+- The `anisRages` use the information from the anistropy matrices generated in the previous part.
+This anisotropy range has to be downscaled a bit, therefore a factor is applied prior to the the anisotropy matrices:
+
+```{r}
+v1_qu = setCgram(type = vg.Gau, sill = s1, anisRanges = 0.015*A1_qu)
+v2_fsp = setCgram(type = vg.Gau, sill = s2, anisRanges = 0.015*A2_fsp)
+v3_bi = setCgram(type = vg.Exp, sill = s3, anisRanges = 0.015*A3_bi)
+v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
+v5_joints = setCgram(type = vg.Exp, sill = s5, anisRanges = 0.05*A5_joints)
+```
+
+Plotting one of these variograms shows this:
+```{r}
+plot(v3_bi)
+```
+
+The factor applied to the anisotropy matrix is can be changed according to target patterns, which should be achieved by the simulation.
+
+The variograms have now to be stacked:
+```{r}
+vm = v5_joints + v4_frac + v3_bi + v2_fsp + v1_qu
+```
+Please note, that the order of the stacking is important.
+The first phase will have least influence, the last phase will have the main influence.
+
+The visualisation of the stacked variograms looks as follows:
+
+```{r}
+plot(vm)
+```
+
+### Apply Turning Bands algorithm
+
+In a first step, a grid of nodes has to be generated, to determine the number and extend of the pixels of the finalized "image".
+This can be done by `expand.grid`:
+```{r}
+# make grid
+resol=200
+Xgrid = as.matrix(expand.grid(x = seq(0,1, length.out = resol),
+ y = seq(0,1, length.out = resol),
+ z = 0))
+```
+
+Pleas note: `Xgrid` needs to have the class `matrix`, and the third dimension is muted by using only zeros for the z-direction.
+
+The data are generated by:
+```{r, message=FALSE}
+set.seed(122487) # 122487
+rock1 = gsi.TurningBands(X = Xgrid,
+ vgram = vm,
+ nBands = 1000)
+colnames(rock1)[-(1:3)] = phases # apply the names
+```
+
+This gives back a large matrix which looks like this:
+```{r}
+rock = data.frame(rock1)
+head(rock)
+```
+
+For each 'row' in the resulting image, i.e. keep a y-value constant, and look along the x-axis, the five phases have a certain probability to appear.
+For example for y = 0.10025 :
+```{r, echo=F}
+rock_melt = melt(data.table(rock), id.vars = c("x", "y", "z"),
+ variable.name = "phase",
+ value.name = "relation")
+rock_melt[, phase := factor(phase, levels = phases)]
+ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase)) +
+ theme_bw() +
+ #geom_hline(yintercept = 0, color = "red") +
+ geom_line() +
+ scale_color_brewer(type = "qual", palette = "Set1")
+```
+
+For each pixel, i.e. for each x/y combination respectively each row in `rock`, the phase with maximum value is chosen:
+```{r}
+rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x))
+dim(rocksims) = c(resol, resol)
+```
+
+### Visualization
+
+The resulting image is:
+```{r, echo=F}
+# somehow ggplot doesnt work...
+# phase_info = data.table(n_phase = c(1:5), phase = factor(phases, levels = phases))
+# rocksims_melt = melt(data.table(data.frame(rocksims), keep.rownames = T), id.vars = "rn", value.name = "n_phase")
+# rocksims_melt = rocksims_melt[phase_info, on = "n_phase"]
+# # plot
+# ggplot(rocksims_melt, aes(y = rn, x = variable, fill = phase)) +
+# geom_tile() +
+# coord_equal() +
+# scale_fill_brewer(type = "qual", palette = "Set1") +
+# theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank())
+```
+```{r}
+image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
+legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
+```
+
+
+
+
+If one of the phases should be less pronounced, e.g. here in the example the joints are very frequent, then the phases can be regulated by weights.
+The weights simply increase or diminish the values of the `rock` matrix for the phase.
+For example:
+```{r}
+ph_weights = c(0,0,0,0,-1) # diminish only the joints
+rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x + ph_weights))
+dim(rocksims) = c(resol, resol)
+```
+
+Then the the line is pushed down:
+```{r, echo=F}
+rock_melt = melt(data.table(rock)[, joints := joints -1], id.vars = c("x", "y", "z"),
+ variable.name = "phase",
+ value.name = "relation")
+rock_melt[, phase := factor(phase, levels = phases)]
+ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase)) +
+ theme_bw() +
+ #geom_hline(yintercept = 0, color = "red") +
+ geom_line() +
+ scale_color_brewer(type = "qual", palette = "Set1")
+```
+
+And then the pattern looks like:
+
+```{r}
+image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
+legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
+```
+
+### Example with more realistic fractures
+
+If single the phases should be changed then only certain steps have to be retraced.
+For example, in the simulation above, the geologist would expect, that the fractures rather follow the biotite layers, instead of being orthogonal to the joints.
+
+Changing the orientation of the fractures:
+```{r}
+p4_frac = anis2D.par2A(ratio = .01, angle = 95)
+A4_frac = tcrossprod(p4_frac)
+A4_frac
+```
+
+Then recalculate the variogram for the fractures and stack all variograms again:
+```{r}
+v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
+vm = v5_joints + v4_frac + v3_bi + v2_fsp + v1_qu
+```
+
+Plot:
+```{r}
+plot(vm)
+```
+
+Applying the algorithm of the turning bands:
+```{r}
+set.seed(384952)
+rock2 = gsi.TurningBands(X = Xgrid, vgram = vm, nBands = 1000)
+colnames(rock2)[-(1:3)] = phases
+```
+
+Get the phases per pixel, including the weights for the joints:
+```{r}
+ph_weights = c(0,0,0,0,-1) # diminish only the joints
+rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
+dim(rocksims) = c(resol, resol)
+```
+
+and then visualize:
+```{r}
+image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
+legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
+```
+
+In this plot, the fracture seem to be too big.
+Again, this can be easily resolved by downweighing the fracturesa (and joints) and increase a bit the Quartz:
+```{r}
+ph_weights = c(.1,0,0,-1.2,-1) # diminish joints and fracture, increase quartz
+rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
+dim(rocksims) = c(resol, resol)
+```
+
+```{r}
+image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
+legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
+```