diff --git a/Formulas_Reference/Formelsammlung_Beispielrechnung.pdf b/Formulas_Reference/Formelsammlung_Beispielrechnung.pdf
index 1dcd38be2abf8949cac5ecdb7720bb0c7aeba9f1..19a7192e1dd76c8035108a6141d9350d7d7ca254 100644
Binary files a/Formulas_Reference/Formelsammlung_Beispielrechnung.pdf and b/Formulas_Reference/Formelsammlung_Beispielrechnung.pdf differ
diff --git a/Formulas_Reference/Formelsammlung_Beispielrechnung.tex b/Formulas_Reference/Formelsammlung_Beispielrechnung.tex
index 4363986468ec945cf6d64d202f5d0f8e78584ef4..d62e41beea98b7e097851c717faef41bb73c803f 100644
--- a/Formulas_Reference/Formelsammlung_Beispielrechnung.tex
+++ b/Formulas_Reference/Formelsammlung_Beispielrechnung.tex
@@ -92,7 +92,8 @@ Bemerkung: Sämtliche Barwerte werden rekursiv bestimmt, daher werden alle Forme
  $Beg$ & Versicherungsbeginn (Datum, TODO)\\[0.5em]
 
  $q_{x+t}$ & einjährige Sterbewahrscheinlichkeit der versicherten Person (aus $YOB$ und $x$ bestimmt) & \texttt{contract\$params\$transitionProbabilities\$q}\\
- $p_{x+t}$ & einjährige Überlebenswahrscheinlichkeit, $p_{x+t} = 1-q_{x+t}$ & \texttt{contract\$params\$transitionProbabilities\$p} \\
+ $i_{x+t}$ & einjährige Erkrankungswahrscheinlichkeit (z.B. Dread Disease, Krebs, etc.) & \texttt{contract\$params\$transitionProbabilities\$i} \\
+ $p_{x+t}$ & einjährige Überlebenswahrscheinlichkeit als Aktiver, $p_{x+t} = 1-q_{x+t}-i_{x+1}$ & \texttt{contract\$params\$transitionProbabilities\$p} \\
  $\omega$  & Höchstalter gemäß der benutzten Sterbetafel  & \texttt{getOmega(tarif\$mortalityTable)} \\[0.5em]
  
 
@@ -361,30 +362,39 @@ Jede Kostenart ($\alpha$/Zillmer/$\beta$/$\gamma$/$\tilde{\gamma}$) und Bemessun
 
 \pagebreak
 
-\section{Cashflows}
+\section{Normierte Cashflows (auf Einheitswerte)}
 
-\begin{longtable}{lp{6.5cm}v{2.5cm}v{2.5cm}v{2.5cm}}
-& Beschreibung & LR & ALV & ELV \\\hline
-$pr_t$  \dots          & Prämienzahlungen (vorschüssig) zu $t$ & $\delta(t<m)$ & $\delta(t<m)$ & $\delta(t<m)$ \\
-$PS$                   & Prämiensumme, $PS=\sum_{t=0}^n pr_t$\\[1em]
+% \begin{sidewaystable}
+\begin{center}
+ 
 
-$\ddot{e}_t$  \dots    & Erlebenszahlungen vorschüssig zu $t$ & $\delta(l+g\le t< n)$ & 0 & $\delta(t=n)$\\
-$e_t$  \dots           & Erlebenszahlungen nachschüssig zu $t+1$& $\delta(l+g\le t< n)$ & 0 & $\delta(t=n)$\\
-$\ddot{e}_t^{*}$ \dots & garantierte Zahlungen vorschüssig zu $t$ & $\delta(l\le t< l+g)$ & 0 & 0 \\
-$e_t^{*}$  \dots       & garant. Zahlungen nachschüssig zu $t+1$ & $\delta(l\le t< l+g)$ & 0 & 0 \\[1em]
+% \rotatebox{90}{
+\begin{tabular}{lp{6.5cm}v{1.6cm}v{1.5cm}v{1.5cm}v{1cm}v{1cm}}
+& Beschreibung & LR & ALV & ELV & T-F & DD\\\hline\\[0.01em]
+$pr_t$  \dots          & Prämienzahlungen (vorschüssig) zu $t$ & $t<m$ & $t<m$ & $t<m$ & $t<m$ & $t<m$ \\
+$PS$                   & Prämiensumme, $PS=\sum_{t=0}^n pr_t$\\[0.5em]\hline\\[0em]
 
-$a_t$  \dots           & Ablebenszahlung zu $t+1$ & 0 & $\delta(l\le t < n)$ & 0 \\
-$a_t^{(RG)}$           & Ablebenszahlungen für PRG zu $t+1$ (Ableben im Jahr $t$) & $\min(t+1,m,f)$ & 0 & $\min(t+1,m,f)$ \\
+$\ddot{e}_t$  \dots    & Erlebenszahlungen vorschüssig zu $t$ & $(l+g)...n$ & 0 & $t=n$ & 0 & 0 \\
+$e_t$  \dots           & Erlebenszahlungen nachschüssig zu $t+1$& $(l+g)...n$ & 0 & $t=n$ & 0 & 0 \\
+$\ddot{e}_t^{*}$ \dots & garantierte Zahlungen vorschüssig zu $t$ & $(l+g)...n$ & 0 & 0  & $t=n$ & 0 \\
+$e_t^{*}$  \dots       & garant. Zahlungen nachschüssig zu $t+1$ & $(l+g)...n$ & 0 & 0  & 0 & 0 \\[0.5em]\hline\\[0em]
 
-\end{longtable}
+$a_t$  \dots           & Ablebenszahlung zu $t+1$ & 0 & $l...n$ & 0 & 0 & 0 \\
+$d_t$  \dots           & Zahlung zu $t+1$ bei Eintritt der Erkrankung zwischen $t$ und $t+1$ & 0 & 0 & 0 & 0& $l...n$ \\[0.5em]\hline\\[0em]
+$a_t^{(RG)}$           & Ablebenszahlungen für PRG zu $t+1$ (Ableben im Jahr $t$) & $\min(t+1, m, f)$ & 0 & $\min(t+1,m,f)$ & 0 & 0 \\[0.5em]\hline
+
+\end{tabular}
+% }
+\end{center}
+% \end{sidewaystable}
 
 Die Cash-Flows können auch in Matrixform dargestellt werden:
 \begin{align*}
 %
  \overrightarrow{CF}^L_t &= \left(
  \begin{matrix} % TODO: Add an indication that the first row is in advance, the second in arrears
-pr_t & \ddot{e}_t^{*} & \ddot{e}_t & 0 & 0 \\
-pr^{(nachsch)}_t & e_t^{*} & e_t & a_t & a_t^{(RG)}\\
+pr_t & \ddot{e}_t^{*} & \ddot{e}_t & 0 & 0 & 0 \\
+pr^{(nachsch)}_t & e_t^{*} & e_t & a_t & d_t & a_t^{(RG)}\\
 \end{matrix}
  \right)
 % 
@@ -522,10 +532,16 @@ A^{(RG+)}_\xn(t) &= A^{(RG)}_\xn(t) - A^{(RG-)}_\xn(t)
 \intertext{Dieser Wert ist rekursiv nicht darstellbar.}
 \end{align*}
 
+\subsection{Disease-Barwert}
+\begin{align*} 
+A^d_\xn(t) &= i_{x+t} \cdot v \cdot d_t + p_{x+t} \cdot v \cdot A^d_\xn(t+1)\\
+\end{align*}
+
 \subsection{Leistungsbarwert}
 
 \begin{align*}
-BW^L_\xn(t) &= E_\xn(t) + A_\xn(t) + (1+\rho^{RG}) \cdot A^{(RG)}_\xn(t)\cdot BP_\xn
+BW^L_\xn(t) &= E_\xn(t) + A_\xn(t) + A^d_\xn(t) &\text{ohne Prämienrückgewähr}\\
+BW^{L,RG}_\xn(t) &= BW^L_\xn(t) + (1+\rho^{RG}) \cdot A^{(RG)}_\xn(t)\cdot BP_\xn &\text{inkl. Prämienrückgewähr}
 \end{align*}
 
 
@@ -556,7 +572,8 @@ allen Termen der Subscript $\xn$ unterlassen):
 \begin{align*}
   \overrightarrow{BW}^L(t) &= \left(
  \begin{matrix}
-    P(t), & E^{Gar}(t), & E(t), & A(t), & A^{(RG)}(t)
+    P(t), & %E^{Gar}(t), & E(t), & A(t), 
+    BW^L(t), & A^{(RG)}(t)
  \end{matrix}
  \right)
 % 
@@ -585,7 +602,9 @@ VK^{(VS)}(t) & VK^{(PS)}(t) & -\\
 
 \subsection{Nettoprämie:}
 \begin{align*}
-NP_\xn &= \frac{E_\xn(0) + A_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\xn(0) \cdot BP_\xn}{P_\xn(0)} \cdot \left(1+\rho\right)
+NP_\xn &= \frac{%E_\xn(0) + A_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\xn(0) \cdot BP_\xn}
+BW^{L}_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\xn(0) \cdot BP_\xn}
+{P_\xn(0)} \cdot \left(1+\rho\right)
 \end{align*}
 
 
@@ -594,7 +613,8 @@ NP_\xn &= \frac{E_\xn(0) + A_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\x
 \begin{align*}
 ZP_\xn &= 
 \frac{%
-  NP_\xn\cdot P_\xn(0) + 
+%   NP_\xn\cdot P_\xn(0) +
+  \overbrace{\left(BW^{L}_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\xn(0) \cdot BP_\xn\right) \cdot (1+\rho)}^{=NP_\xn\cdot P_\xn(0)} +
   ZK^{(VS)}_\xn(0) + 
   ZK^{(PS)}_\xn(0) \cdot BP_\xn \cdot PS + 
   ZK^{(BP)}_\xn(0)\cdot BP_\xn
@@ -609,7 +629,8 @@ Varianten:
 \begin{align*}
 ZP_\xn &= 
 \left[%
-  NP_\xn\cdot P_\xn(0) + 
+%   NP_\xn\cdot P_\xn(0) + 
+  BW^{L,RG}_\xn(0) \cdot (1+\rho) +
   \left(ZK^{(VS)}_\xn(0) + IK^{(VS)}_\xn(0) + VK^{(VS)}_\xn(0)\right) +  \right.\\
   &\qquad\left.\left(ZK^{(PS)}_\xn(0) + IK^{(PS)}_\xn(0) + VK^{(PS)}_\xn(0)\right) \cdot BP_\xn \cdot PS + \right.\\
   &\qquad\left.\left(ZK^{(BP)}_\xn(0) + IK^{(BP)}_\xn(0) + VK^{(BP)}_\xn(0)\right) \cdot BP_\xn
@@ -621,8 +642,8 @@ ZP_\xn &=
  \item % REFERENCE: ME
     Prämienrückgewähr proportional zu Zillmerprämie (für Berechnung der Zillmerprämie):
 \begin{align*}
-ZP_\xn &= \frac{E_\xn(0) + A_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\xn(0) \cdot ZP_\xn}{P_\xn(0)} \cdot \left(1+\rho\right)\\
-ZP_\xn &= \frac{E_\xn(0) + A_\xn(0) +  \cdot ZP_\xn}{P_\xn(0) - \left(1+\rho^{RG}\right)\cdot A^{(RG)}_\xn(0)\cdot \left(1+\rho\right)} \cdot \left(1+\rho\right)
+ZP_\xn &= \frac{BW^L_\xn(0) + \left(1+\rho^{RG}\right) \cdot A^{(RG)}_\xn(0) \cdot ZP_\xn}{P_\xn(0)} \cdot \left(1+\rho\right)\\
+ZP_\xn &= \frac{BW^L_\xn(0)}{P_\xn(0) - \left(1+\rho^{RG}\right)\cdot A^{(RG)}_\xn(0)\cdot \left(1+\rho\right)} \cdot \left(1+\rho\right)
 \end{align*}
  
 \end{itemize}
@@ -632,7 +653,7 @@ ZP_\xn &= \frac{E_\xn(0) + A_\xn(0) +  \cdot ZP_\xn}{P_\xn(0) - \left(1+\rho^{RG
 \begin{align*}
 BP_\xn &= \frac%
 { 
-   \left(E_\xn(0) + A_\xn(0)\right)\cdot\left(1+\rho\right) + 
+   BW^L_\xn(0)\cdot\left(1+\rho\right) + 
    \left( AK^{(VS)}_\xn(0) + IK^{(VS)}_\xn(0) + VK^{(VS)}_\xn(0) \right)
 }%
 {
@@ -654,13 +675,15 @@ Für die Berechnung der Prämien können die Koeffizienten der jeweiligen Barwer
 \begin{sidewaystable}
 \centering
 
-\newenvironment{benarray}{\big(\begin{array}{*4{C{3.4em}}C{12em}}}{\end{array}\big)}
+% \newenvironment{benarray}{\big(\begin{array}{*4{C{3.4em}}C{12em}}}{\end{array}\big)}
+\newenvironment{benarray}{\big(\begin{array}{C{4em}C{5em}C{12em}}}{\end{array}\big)}
 \newenvironment{costarray}{\left(\begin{array}{*3{C{4.7em}}}}{\end{array}\right)}
 % TODO: laufende alpha-Kosten
 \begin{tabular}{||ll|c|c||}\hline\hline
  & & Leistungen & Kosten \\ \hline\hline
  Terme & &  
-  $\begin{benarray}P_\xn(t) & E^{Gar}_\xn(t) & E_\xn(t) & A_\xn(t)  & A_\xn^{(RG)}(t)\end{benarray}$
+%   $\begin{benarray}P_\xn(t) & E^{Gar}_\xn(t) & E_\xn(t) & A_\xn(t)  & A_\xn^{(RG)}(t)\end{benarray}$
+  $\begin{benarray}P_\xn(t) & BW^L_\xn(t)  & A_\xn^{(RG)}(t)\end{benarray}$
  &
  $\begin{costarray}
 AK^{(VS)}_\xn(t) & AK^{(PS)}_\xn(t)  & AK^{(BP)}_\xn(t) \\
@@ -673,18 +696,18 @@ VK^{frei}_\xn(t) & - & -\\
 
 
 Nettoprämie & Zähler & 
-  $\begin{benarray}0 & 1+\rho & 1+\rho & 1+\rho & \left(1+\rho^{RG}\right)\cdot BP_\xn \cdot \left(1+\rho\right)\end{benarray}$
+  $\begin{benarray}0 & 1+\rho & \left(1+\rho^{RG}\right)\cdot BP_\xn \cdot \left(1+\rho\right)\end{benarray}$
  & -
 \\
 
  & Nenner & 
-  $\begin{benarray} 1 & 0 & 0 & 0 & 0\end{benarray}$
+  $\begin{benarray} 1 & 0 & 0 \end{benarray}$
  & - 
 \\\hline
 
 
 Zillmerprämie & Zähler & 
-  $\begin{benarray}0 & 1+\rho & 1+\rho & 1+\rho & \left(1+\rho^{RG}\right)\cdot BP_\xn \cdot \left(1+\rho\right)\end{benarray}$
+  $\begin{benarray}0 & 1+\rho & \left(1+\rho^{RG}\right)\cdot BP_\xn \cdot \left(1+\rho\right)\end{benarray}$
  &
  $\begin{costarray}
 0 & 0 & 0 \\
@@ -696,13 +719,13 @@ Zillmerprämie & Zähler &
 \\
 
  & Nenner & 
-  $\begin{benarray} 1 & 0 & 0 & 0 & 0\end{benarray}$
+  $\begin{benarray} 1 & 0 & 0\end{benarray}$
  &
 -
 \\\hline
 
 Bruttoprämie & Zähler & 
-   $\begin{benarray}0 & 1+\rho & 1+\rho & 1+\rho & 0\end{benarray}$
+   $\begin{benarray}0 & 1+\rho & 0\end{benarray}$
  &
  $\begin{costarray}
 1 & 0 & 0 \\
@@ -714,7 +737,7 @@ Bruttoprämie & Zähler &
 \\
 
  & Nenner & 
-  $\begin{benarray}1 & 0 & 0 & 0 & -(1+\rho)\cdot(1+\rho^{RG}) \end{benarray}$
+  $\begin{benarray}1 & 0 & -(1+\rho)\cdot(1+\rho^{RG}) \end{benarray}$
  &
  $\begin{costarray}
 0 & -PS & -1 \\
diff --git a/R/HelperFunctions.R b/R/HelperFunctions.R
index 75d525b75f90e0b0954cf23fdca66510c2301507..fc69ee54208177da0011dc3a109abfd7aa4217ac 100644
--- a/R/HelperFunctions.R
+++ b/R/HelperFunctions.R
@@ -1,8 +1,10 @@
-calculatePVSurvival = function(q, advance, arrears=c(0), ..., m=1, mCorrection = list(alpha=1, beta=0), v=1) {
+# Caution: px is not neccessarily 1-qx, because we might also have dread diseases so that px=1-qx-ix! However, the ix is not used for the survival present value
+calculatePVSurvival = function(px=1-qx, qx=1-px, advance, arrears=c(0), ..., m=1, mCorrection = list(alpha=1, beta=0), v=1) {
   # assuming advance and arrears have the same dimensions...
   init = advance[1]*0;
-  l = max(length(q), length(advance), length(arrears));
-  q = pad0(q, l, value=1);
+  l = max(length(qx), length(advance), length(arrears));
+  q = pad0(qx, l, value=1);
+  p = pad0(px, l, value=0);
   advance = pad0(advance, l, value=init);
   arrears = pad0(arrears, l, value=init);
 
@@ -11,11 +13,10 @@ calculatePVSurvival = function(q, advance, arrears=c(0), ..., m=1, mCorrection =
   res = rep(0, l+1);
   for (i in l:1) {
     # coefficients for the payemtns(including corrections for payments during the year (using the alpha(m) and beta(m)):
-    p = (1-q[i]);
-    advcoeff = mCorrection$alpha - mCorrection$beta*(1-p*v);
-    arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m)*(1-p*v);
+    advcoeff = mCorrection$alpha - mCorrection$beta*(1-p[i]*v);
+    arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m)*(1-p[i]*v);
     # The actual recursion:
-    res[i] = advance[i]*advcoeff + arrears[i]*arrcoeff + v*(1-q[i])*res[i+1];
+    res[i] = advance[i]*advcoeff + arrears[i]*arrcoeff + v*p[i]*res[i+1];
   }
   res[1:l]
 }
@@ -23,9 +24,9 @@ calculatePVSurvival = function(q, advance, arrears=c(0), ..., m=1, mCorrection =
 
 
 # TODO: So far, we are assuming, the costs array has sufficient time steps and does not need to be padded!
-calculatePVCosts = function(q, costs, ..., v=1) {
-  l = max(length(q), dim(costs)[1]);
-  q = pad0(q, l, value=1);
+calculatePVCosts = function(px=1-qx, qx=1-px, costs, ..., v=1) {
+  l = max(length(qx), dim(costs)[1]);
+  p = pad0(px, l, value=0);
   costs = costs[1:l,,];
 
   # Take the array structure from the cash flow array and initialize it with 0
@@ -34,27 +35,48 @@ calculatePVCosts = function(q, costs, ..., v=1) {
   # Backward recursion starting from the last time:
   for (i in l:1) {
     # cat("values at iteration ", i, ": ", v, q[i], costs[i,,], prev);
-    res[i,,] = costs[i,,] + v*(1-q[i])*prev;
+    res[i,,] = costs[i,,] + v*p[i]*prev;
     prev=res[i,,];
   }
   res
 }
 
-calculatePVDeath = function(q, benefits, ..., v=1) {
+calculatePVDeath = function(px, qx, benefits, ..., v=1) {
+  init = benefits[1]*0; # Preserve the possible array structure of the benefits -> vectorized calculations possible!
+  l = max(length(qx), length(benefits));
+  q = pad0(qx, l, value=1);
+  p = pad0(px, l, value=0);
+  benefits = pad0(benefits, l, value=init);
+
+  # TODO: Make this work for matrices (i.e. currently benefits are assumed to be one-dimensional vectors)
+  # TODO: Replace loop by better way (using Reduce?)
+  res = rep(init, l+1);
+  for (i in l:1) {
+    # Caution: p_x is not neccessarily 1-q_x, because we might also have dread diseases, so that px=1-qx-ix!
+    res[i] = v*q[i]*benefits[i] + v*p[i]*res[i+1];
+  }
+  res[1:l]
+}
+
+calculatePVDisease = function(px=1-qx-ix, qx=1-ix-px, ix=1-px-qx, benefits, ..., v=1) {
   init = benefits[1]*0;
-  l = max(length(q), length(benefits));
-  q = pad0(q, l, value=1);
+  l = min(length(ix), length(qx), length(benefits));
+  qx = pad0(qx, l, value=1);
+  ix = pad0(ix, l, value=0);
+  px = pad0(px, l, value=0);
   benefits = pad0(benefits, l, value=init);
 
   # TODO: Make this work for matrices (i.e. currently benefits are assumed to be one-dimensional vectors)
   # TODO: Replace loop by better way (using Reduce?)
   res = rep(init, l+1);
   for (i in l:1) {
-    res[i] = v*q[i]*benefits[i] + v*(1-q[i])*res[i+1];
+    res[i] = v*ix[i]*benefits[i] + v*px[i]*res[i+1];
   }
   res[1:l]
 }
 
+
+
 getSavingsPremium = function(reserves, v=1) {
   pad0(reserves[-1], length(reserves))*v - reserves
 }
diff --git a/R/InsuranceTarif.R b/R/InsuranceTarif.R
index 9b0442b75de00c8f2fc35787db3cb7f140c3cba2..92dfc665d05c31b3ccec50340bd11ca60272e632 100644
--- a/R/InsuranceTarif.R
+++ b/R/InsuranceTarif.R
@@ -3,7 +3,7 @@ library(lifecontingencies)
 library(objectProperties)
 library(foreach)
 
-TariffTypeEnum = setSingleEnum("TariffType", levels = c("annuity", "wholelife", "endowment", "pureendowment", "terme-fix"))
+TariffTypeEnum = setSingleEnum("TariffType", levels = c("annuity", "wholelife", "endowment", "pureendowment", "terme-fix", "dread-disease"))
 PaymentTimeEnum = setSingleEnum("PaymentTime", levels = c("in advance", "in arrears"))
 #PaymentCountEnum = setSingleEnum(PaymentCount, levels = c(1,2,3))
 
@@ -39,6 +39,7 @@ InsuranceTarif = R6Class(
 
     states = c("alive", "dead"),
     mortalityTable = NULL,
+    invalidityTable = NULL,
     i = 0, # guaranteed interest rate
     v = 1, # discount factor
     tariffType = TariffTypeEnum("wholelife"), # possible values: annuity, wholelife, endowment, pureendowment, terme-fix
@@ -74,7 +75,7 @@ InsuranceTarif = R6Class(
       ),
 
 
-    initialize = function(name = NULL, mortalityTable = NULL, i = NULL, type = "wholelife", ..., features = list(), premiumPeriod = NULL, premiumFrequencyOrder = 0, benefitFrequencyOrder = 0, costs, surrenderValueCalculation) {
+    initialize = function(name = NULL, mortalityTable = NULL, i = NULL, type = "wholelife", ..., invalidityTable=NULL, features = list(), premiumPeriod = NULL, premiumFrequencyOrder = 0, benefitFrequencyOrder = 0, costs, surrenderValueCalculation) {
       if (!missing(name))           self$name = name;
       if (!missing(mortalityTable)) self$mortalityTable = mortalityTable;
       if (!missing(i))              self$i = i;
@@ -86,11 +87,18 @@ InsuranceTarif = R6Class(
       if (!missing(premiumPeriod))  self$defaultPremiumPeriod = premiumPeriod;
       if (!missing(features))       self$features = c(features, self$features);
       if (!missing(surrenderValueCalculation)) self$surrenderValueCalculation = surrenderValueCalculation;
+      if (!missing(invalidityTable)) self$invalidityTable = invalidityTable;
 
       self$v = 1/(1+self$i);
 
       cat(paste0("Initializing Insurance Tarif ", self$name, "...\n"));
     },
+
+    # Merge a possibly passed loadings override with this tariff's default loadings:
+    getLoadings = function(..., loadings=list()) {
+      c(loadings, self$loadings)
+    },
+
     getAges = function(age, ..., YOB = 2000) {
       ages = ages(self$mortalityTable, YOB = YOB);
       if (age > 0) {
@@ -104,7 +112,16 @@ InsuranceTarif = R6Class(
       if (age > 0) {
         q    = q[-age:-1];
       }
-      df = data.frame(age=ages, q=q, p=1-q, row.names = ages-age)
+      if (!is.null(self$invalidityTable)) {
+        i = deathProbabilities(self$invalidityTable, YOB=YOB);
+        if (age > 0) {
+          i    = i[-age:-1];
+        }
+      } else {
+        i = rep(0, length(q));
+      }
+      i = pad0(i, length(q));
+      df = data.frame(age=ages, q=q, i=i, p=1-q-i, row.names = ages-age)
       df
     },
     getBasicCashFlows = function(age, ..., guaranteed = 0, policyPeriod = inf, deferral = 0, maxAge = getOmega(self$mortalityTable)) {
@@ -112,7 +129,8 @@ InsuranceTarif = R6Class(
       cf = list(
         guaranteed = rep(0, maxlen+1),
         survival = rep(0, maxlen+1),
-        death = rep(0, maxlen+1)
+        death = rep(0, maxlen+1),
+        disease = rep(0, maxlen+1)
       );
       if (self$tariffType == "annuity") {
         # guaranteed payments exist only with annuities (first n years of the payment)
@@ -120,6 +138,8 @@ InsuranceTarif = R6Class(
         cf$survival = c(rep(0, deferral + guaranteed), rep(1, max(0, maxlen - deferral - guaranteed)), 0)
       } else if (self$tariffType == "terme-fix") {
         cf$guaranteed = c(rep(0, policyPeriod), 1);
+      } else if (self$tariffType == "dread-disease") {
+        cf$disease = c(rep(0, deferral), rep(1, maxlen - deferral), 0);
       } else {
         if (self$tariffType == "endowment" || self$tariffType == "pureendowment") {
           cf$survival = c(rep(0, policyPeriod), 1);
@@ -147,6 +167,7 @@ InsuranceTarif = R6Class(
         survival_advance = zeroes,
         survival_arrears = zeroes,
         death_SumInsured = zeroes,
+        disease_SumInsured = zeroes,
         death_GrossPremium = zeroes,
         death_Refund_past = zeroes,
         death_PremiumFree = zeroes,
@@ -172,6 +193,7 @@ InsuranceTarif = R6Class(
 
       # Death Benefits
       cf$death_SumInsured = pad0(basicCashFlows$death, cflen);
+      cf$disease_SumInsured = pad0(basicCashFlows$disease, cflen);
       cf$death_PremiumFree = cf$death_SumInsured;
       # premium refund
       if (self$premiumRefund != 0) {
@@ -212,23 +234,27 @@ InsuranceTarif = R6Class(
     presentValueCashFlows = function(cashflows, age, ..., premiumFrequency = 1, benefitFrequency = 1, maxAge = getOmega(self$mortalityTable)) {
       len = length(cashflows$premiums_advance);
       qq = self$getTransitionProbabilities (age, ...);
-      q = pad0(qq$q, len);
+      qx = pad0(qq$q, len);
+      ix = pad0(qq$i, len);
+      px = pad0(qq$p, len);
       benefitFrequencyCorrection = correctionPaymentFrequency(m = benefitFrequency, i = self$i, order = self$benefitFrequencyOrder);
       premiumFrequencyCorrection = correctionPaymentFrequency(m = premiumFrequency, i = self$i, order = self$premiumFrequencyOrder);
 
-      pvRefund = calculatePVDeath (q, cashflows$death_GrossPremium, v=self$v);
-      pvRefundPast = calculatePVDeath (q, cashflows$death_Refund_past, v=self$v) * (cashflows[,"death_GrossPremium"]-cashflows[,"premiums_advance"]);
+      pvRefund = calculatePVDeath (px, qx, cashflows$death_GrossPremium, v=self$v);
+      pvRefundPast = calculatePVDeath (px, qx, cashflows$death_Refund_past, v=self$v) * (cashflows[,"death_GrossPremium"]-cashflows[,"premiums_advance"]);
 
       pv = cbind(
-        premiums = calculatePVSurvival (q, cashflows$premiums_advance, cashflows$premiums_arrears, m=premiumFrequency, mCorrection=premiumFrequencyCorrection, v=self$v),
-        guaranteed = calculatePVSurvival (q*0, cashflows$guaranteed_advance, cashflows$guaranteed_arrears, m=benefitFrequency, mCorrection=benefitFrequencyCorrection, v=self$v),
-        survival = calculatePVSurvival (q, cashflows$survival_advance, cashflows$survival_arrears, m=benefitFrequency, mCorrection=benefitFrequencyCorrection, v=self$v),
-        death_SumInsured = calculatePVDeath (q, cashflows$death_SumInsured, v=self$v),
+        premiums = calculatePVSurvival (px, qx, cashflows$premiums_advance, cashflows$premiums_arrears, m=premiumFrequency, mCorrection=premiumFrequencyCorrection, v=self$v),
+        guaranteed = calculatePVSurvival (px/px, qx*0, cashflows$guaranteed_advance, cashflows$guaranteed_arrears, m=benefitFrequency, mCorrection=benefitFrequencyCorrection, v=self$v),
+        survival = calculatePVSurvival (px, qx, cashflows$survival_advance, cashflows$survival_arrears, m=benefitFrequency, mCorrection=benefitFrequencyCorrection, v=self$v),
+        death_SumInsured = calculatePVDeath (px, qx, cashflows$death_SumInsured, v=self$v),
+        disease_SumInsured = calculatePVDisease(px, qx, ix, cashflows$disease_SumInsured, v=self$v),
         death_GrossPremium = pvRefund,
         death_Refund_past = pvRefundPast,
         death_Refund_future = pvRefund - pvRefundPast,
-        death_PremiumFree = calculatePVDeath (q, cashflows$death_PremiumFree, v=self$v)
+        death_PremiumFree = calculatePVDeath (px, qx, cashflows$death_PremiumFree, v=self$v)
       );
+
       rownames(pv) <- pad0(rownames(qq), len);
       pv
     },
@@ -236,9 +262,10 @@ InsuranceTarif = R6Class(
     presentValueCashFlowsCosts = function(cashflows, age, ..., maxAge = getOmega(self$mortalityTable)) {
       len = dim(cashflows)[1];
       q = self$getTransitionProbabilities (age, ...);
-      q = pad0(q$q, len);
+      qx = pad0(q$q, len);
+      px = pad0(q$p, len);
 
-      pvc = calculatePVCosts(q, cashflows, v=self$v);
+      pvc = calculatePVCosts(px, qx, cashflows, v=self$v);
       pvc
     },
 
@@ -257,10 +284,11 @@ InsuranceTarif = R6Class(
     getAbsCashFlows = function(cashflows, cashflowsCosts, premiums, sumInsured=1, premiumSum=0, ...) {
       refundAddon = self$premiumRefundLoading;
 
+      # TODO: Set up a nice list with coefficients for each type of cashflow, rather than multiplying each item manually (this also mitigates the risk of forgetting a dimension, because then the dimensions would not match, while here it's easy to overlook a multiplication)
       # Multiply each CF column by the corresponding basis
       cashflows[,c("premiums_advance", "premiums_arrears")] = cashflows[,c("premiums_advance", "premiums_arrears")] * premiums[["gross"]];
-      cashflows[,c("guaranteed_advance", "guaranteed_arrears", "survival_advance", "survival_arrears", "death_SumInsured", "death_PremiumFree")] =
-        cashflows[,c("guaranteed_advance", "guaranteed_arrears", "survival_advance", "survival_arrears", "death_SumInsured", "death_PremiumFree")] * sumInsured;
+      cashflows[,c("guaranteed_advance", "guaranteed_arrears", "survival_advance", "survival_arrears", "death_SumInsured", "death_PremiumFree", "disease_SumInsured")] =
+        cashflows[,c("guaranteed_advance", "guaranteed_arrears", "survival_advance", "survival_arrears", "death_SumInsured", "death_PremiumFree", "disease_SumInsured")] * sumInsured;
       cashflows[,c("death_GrossPremium", "death_Refund_past")] = cashflows[,c("death_GrossPremium","death_Refund_past")] * premiums[["gross"]] * (1+refundAddon);
 
       # Sum all death-related payments to "death"  and remove the death_GrossPremium column
@@ -284,8 +312,8 @@ InsuranceTarif = R6Class(
 
       # Multiply each CF column by the corresponding basis
       pv[,"premiums"] = pv[,"premiums"] * premiums[["gross"]];
-      pv[,c("guaranteed", "survival", "death_SumInsured", "death_PremiumFree")] =
-        pv[,c("guaranteed", "survival", "death_SumInsured", "death_PremiumFree")] * sumInsured;
+      pv[,c("guaranteed", "survival", "death_SumInsured", "disease_SumInsured", "death_PremiumFree")] =
+        pv[,c("guaranteed", "survival", "death_SumInsured", "disease_SumInsured", "death_PremiumFree")] * sumInsured;
       pv[,c("death_GrossPremium", "death_Refund_past", "death_Refund_future")] = pv[,c("death_GrossPremium", "death_Refund_past", "death_Refund_future")] * premiums[["gross"]] * (1+refundAddon);
       pv[,c("benefits", "benefitsAndRefund", "alpha", "Zillmer", "beta", "gamma", "gamma_nopremiums")] =
         pv[,c("benefits", "benefitsAndRefund", "alpha", "Zillmer", "beta", "gamma", "gamma_nopremiums")] * sumInsured;
@@ -301,8 +329,8 @@ InsuranceTarif = R6Class(
     presentValueBenefits = function(presentValues, presentValuesCosts, premiums, sumInsured=1, premiumSum=0, ...) {
       refundAddon = self$premiumRefundLoading;
       # TODO: Here we don't use the securityLoading parameter => Shall it be used or are these values to be understood without additional security loading?
-      benefits    = presentValues[,"survival"] + presentValues[,"death_SumInsured"];
-      allBenefits = presentValues[,"survival"] + presentValues[,"death_SumInsured"] + presentValues[,"death_GrossPremium"] * premiums[["unit.gross"]] * (1+refundAddon);
+      benefits    = presentValues[,"survival"] + presentValues[,"death_SumInsured"] + presentValues[,"disease_SumInsured"];
+      allBenefits = presentValues[,"survival"] + presentValues[,"death_SumInsured"] + presentValues[,"disease_SumInsured"] + presentValues[,"death_GrossPremium"] * premiums[["unit.gross"]] * (1+refundAddon);
 
       benefitsCosts = presentValuesCosts[,,"SumInsured"] +
         presentValuesCosts[,,"SumPremiums"] * premiumSum * premiums[["unit.gross"]] +
@@ -331,6 +359,8 @@ InsuranceTarif = R6Class(
       coeff[["SumInsured"]][["benefits"]][["guaranteed"]]       = 1+securityLoading;
       coeff[["SumInsured"]][["benefits"]][["survival"]]         = 1+securityLoading;
       coeff[["SumInsured"]][["benefits"]][["death_SumInsured"]] = 1+securityLoading;
+      coeff[["SumInsured"]][["benefits"]][["disease_SumInsured"]] = 1+securityLoading;
+
       # Premium refund is handled differently for gross and net premiums, because it is proportional to the gross premium
       if (type == "gross") {
         coeff[["Premium"]][["benefits"]][["death_GrossPremium"]] = -(1+refundAddon) * (1+securityLoading);
@@ -496,7 +526,8 @@ InsuranceTarif = R6Class(
       l = dim(reserves)[[1]];
       premium.savings = getSavingsPremium(reserves[,"Zillmer"], self$v) + getSavingsPremium(reserves[,"gamma"], self$v);
       # TODO: Switch to use the Ziller or net or adequate reserve!
-      premium.risk    = self$v * (cashflows[,"death"] - c(reserves[,"Zillmer"][-1], 0)) * pad0(q$q, l);
+      premium.risk    = self$v * (cashflows[,"death"] - c(reserves[,"Zillmer"][-1], 0)) * pad0(q$q, l) +
+        self$v * (cashflows[,"disease_SumInsured"] - c(reserves[,"Zillmer"][-1], 0)) * pad0(q$i, l);
       # premium.risk    = self$v * (cashflows[,"death"] - c(reserves[,"Zillmer"][-1], 0)) * q$q;