diff --git a/src/spicelib/devices/bsim4/Makefile.am b/src/spicelib/devices/bsim4/Makefile.am
index 14b01d542..ddaf18e75 100644
--- a/src/spicelib/devices/bsim4/Makefile.am
+++ b/src/spicelib/devices/bsim4/Makefile.am
@@ -39,6 +39,8 @@ endif
 AM_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/src/include
 AM_CFLAGS = $(STATIC)
 
+AM_LDFLAGS = -lgsl -lgslcblas
+
 MAINTAINERCLEANFILES = Makefile.in
 
 EXTRA_DIST = B4TERMS_OF_USE
diff --git a/src/spicelib/devices/bsim4/b4reliability.c b/src/spicelib/devices/bsim4/b4reliability.c
index a18996868..f7e7a05a0 100644
--- a/src/spicelib/devices/bsim4/b4reliability.c
+++ b/src/spicelib/devices/bsim4/b4reliability.c
@@ -7,6 +7,8 @@ Author: 2015 Francesco Lannutti - July 2015
 #include "bsim4def.h"
 #include "ngspice/sperror.h"
 
+#include <gsl/gsl_fit.h>
+
 static int
 BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mode)
 {
@@ -19,7 +21,7 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
     // Determine if the transistor is ON or OFF
     vds = ckt->CKTstate0 [here->BSIM4vds] ;
     vgs = ckt->CKTstate0 [here->BSIM4vgs] ;
-    von = model->BSIM4type * here->BSIM4von ;
+    von = here->BSIM4von ;
     if (vds >= 0)
     {
 //        printf ("VDS >= 0\tBSIM4type: %d\tBSIM4instance: %s\tVgs: %-.9g\tVon: %-.9g\t", model->BSIM4type, here->BSIM4name, vgs, von) ;
@@ -118,16 +120,16 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
                 // Until now, the device was OFF - Calculate recovery
                 delta = ckt->CKTtime - here->relStruct->time ;
 
+                // Update time and flag - Stress begins
+                here->relStruct->time = ckt->CKTtime ;
+                here->relStruct->IsON = 1 ;
+
                 // Calculate Aging - Giogio Liatis' Model
                 ret = RELMODELcalculateAging ((GENinstance *)here, here->BSIM4modPtr->BSIM4modType, delta, 0) ;
                 if (ret == 1)
                 {
                     return (E_INTERN) ;
                 }
-
-                // Update time and flag - Stress begins
-                here->relStruct->time = ckt->CKTtime ;
-                here->relStruct->IsON = 1 ;
             } else {
                 fprintf (stderr, "Reliability Analysis Error\n") ;
             }
@@ -137,16 +139,16 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
                 // Until now, the device was ON - Calculate stress
                 delta = ckt->CKTtime - here->relStruct->time ;
 
+                // Update time and flag - Recovery begins
+                here->relStruct->time = ckt->CKTtime ;
+                here->relStruct->IsON = 0 ;
+
                 // Calculate Aging - Giorgio Liatis' Model
                 ret = RELMODELcalculateAging ((GENinstance *)here, here->BSIM4modPtr->BSIM4modType, delta, 1) ;
                 if (ret == 1)
                 {
                     return (E_INTERN) ;
                 }
-
-                // Update time and flag - Recovery begins
-                here->relStruct->time = ckt->CKTtime ;
-                here->relStruct->IsON = 0 ;
             } else if (here->relStruct->IsON == 0) {
                 // Until now, the device was OFF - Do NOTHING
                 delta = -1 ;
@@ -155,32 +157,36 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
             }
         }
     } else if (mode == 1) {
-        // In this mode, it doesn't matter if NOW the device is in stress or in recovery, since it's the last timestep
+        // In this mode, it doesn't matter if NOW the device is in stress or in recovery, since it's the latest timestep
         if (here->relStruct->IsON == 1)
         {
             // Calculate stress
             delta = ckt->CKTtime - here->relStruct->time ;
+
+            // Update time and flag - Maybe Optional
+            here->relStruct->time = ckt->CKTtime ;
+            here->relStruct->IsON = 1 ;
+
+            // Calculate Aging - Giorgio Liatis' Model
             ret = RELMODELcalculateAging ((GENinstance *)here, here->BSIM4modPtr->BSIM4modType, delta, 1) ;
             if (ret == 1)
             {
                 return (E_INTERN) ;
             }
-
-            // Update time and flag - Maybe Optional
-            here->relStruct->time = ckt->CKTtime ;
-            here->relStruct->IsON = 1 ;
         } else if (here->relStruct->IsON == 0) {
             // Calculate recovery
             delta = ckt->CKTtime - here->relStruct->time ;
+
+            // Update time and flag - Maybe Optional
+            here->relStruct->time = ckt->CKTtime ;
+            here->relStruct->IsON = 0 ;
+
+            // Calculate Aging - Giorgio Liatis' Model
             ret = RELMODELcalculateAging ((GENinstance *)here, here->BSIM4modPtr->BSIM4modType, delta, 0) ;
             if (ret == 1)
             {
                 return (E_INTERN) ;
             }
-
-            // Update time and flag - Maybe Optional
-            here->relStruct->time = ckt->CKTtime ;
-            here->relStruct->IsON = 0 ;
         } else {
             fprintf (stderr, "Reliability Analysis Error\n") ;
         }
@@ -188,6 +194,50 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
         printf ("DeltaVth: %-.9gmV\t\t", here->relStruct->deltaVth * 1000) ;
         printf ("Device Name: %s\t\t", here->BSIM4name) ;
         printf ("Device Type: %s\n\n", model->BSIM4modName) ;
+
+
+        /* Calculate fitting */
+        double c0, c1, cov00, cov01, cov11, chisq ;
+
+        /* Count how many deltaVth we have */
+        unsigned int i ;
+        RELMODELrelList *current ;
+
+        i = 0 ;
+        current = here->relStruct->deltaVthList ;
+        while (current != NULL)
+        {
+            i++ ;
+            current = current->next ;
+        }
+
+        /* Assign list members to vectors */
+        double *timeFit, *deltaVthFit ;
+
+        timeFit = TMALLOC (double, i) ;
+        deltaVthFit = TMALLOC (double, i) ;
+
+        i = 0 ;
+        current = here->relStruct->deltaVthList ;
+        while (current != NULL)
+        {
+            timeFit [i] = current->time ;
+            deltaVthFit [i] = exp (current->deltaVth) ;
+            printf ("Time: %-.9g\n", current->time) ;
+            printf ("DeltaVth: %-.9g\n", current->deltaVth) ;
+            printf ("DeltaVth (exp): %-.9g\n", exp (current->deltaVth)) ;
+            i++ ;
+            current = current->next ;
+        }
+
+        /* Execute Fitting */
+        gsl_fit_linear (timeFit + 1, 1, deltaVthFit + 1, 1, i - 1, &c0, &c1, &cov00, &cov01, &cov11, &chisq) ;
+
+        /* Perform Extrapolation */
+        double yFit, yFitErr ;
+        gsl_fit_linear_est (2e-6, c0, c1, cov00, cov01, cov11, &yFit, &yFitErr) ;
+
+        printf ("\n\tFitting -> %-.9gs\t%-.9gmV\n\n", 31536000.0, log (yFit) * 1000) ;
     } else {
         fprintf (stderr, "Reliability Analysis Error\n") ;
     }
diff --git a/src/spicelib/devices/bsim4/b4set.c b/src/spicelib/devices/bsim4/b4set.c
index 687aad5b3..311d50f56 100644
--- a/src/spicelib/devices/bsim4/b4set.c
+++ b/src/spicelib/devices/bsim4/b4set.c
@@ -2663,6 +2663,7 @@ do { if((here->ptr = SMPmakeElt(matrix,here->first,here->second))==(double *)NUL
                 here->relStruct->deltaVth = 0 ;
                 here->relStruct->t_star = 0 ;
                 here->relStruct->IsON = -1 ;
+                here->relStruct->deltaVthList = NULL ;
             }
 #endif
 
diff --git a/src/spicelib/devices/relmodel/relmodelcalcaging.c b/src/spicelib/devices/relmodel/relmodelcalcaging.c
index 8d30c515c..97b0b1adf 100644
--- a/src/spicelib/devices/relmodel/relmodelcalcaging.c
+++ b/src/spicelib/devices/relmodel/relmodelcalcaging.c
@@ -13,6 +13,39 @@ Author: Francesco Lannutti - July 2015
 //#define CONSTepsZero (8.854214871e-12)   /* epsilon zero F/m */
 //#define CONSTepsSiO2 (3.4531479969e-11)  /* epsilon SiO2 F/m */
 
+static int
+listInsert (RELMODELrelList **list, double time, double deltaVth)
+{
+    RELMODELrelList *current ;
+    RELMODELrelList *previous ;
+
+    /* Loop until the end of the list */
+    previous = NULL ;
+    current = *list ;
+    while (current != NULL)
+    {
+        previous = current ;
+        current = current->next ;
+    }
+
+    /* Insert the new element into the list */
+    if (previous == NULL)
+    {
+        *list = TMALLOC (RELMODELrelList, 1) ;
+        current = *list ;
+    } else {
+        current = TMALLOC (RELMODELrelList, 1) ;
+        previous->next = current ;
+    }
+
+    /* Populate the element */
+    current->time = time ;
+    current->deltaVth = deltaVth ;
+    current->next = NULL ;
+
+    return 0 ;
+}
+
 int
 RELMODELcalculateAging (GENinstance *inInstance, int modType, double t_aging, unsigned int stress_or_recovery)
 {
@@ -50,6 +83,16 @@ RELMODELcalculateAging (GENinstance *inInstance, int modType, double t_aging, un
         here->relStruct->deltaVth = here->relStruct->deltaVth * log (1 + (1.718 / (1 + pow ((t_aging / relmodel->RELMODELtau_e), relmodel->RELMODELbeta1)))) ;
     }
 
+    /* Insert 'here->relStruct->deltaVth' into the list for the later fitting */
+    listInsert (&(here->relStruct->deltaVthList), here->relStruct->time, here->relStruct->deltaVth) ;
+//    RELMODELrelList *temp ;
+//    temp = TMALLOC (RELMODELrelList, 1) ;
+//    temp->time = here->relStruct->time ;
+//    temp->deltaVth = here->relStruct->deltaVth ;
+//    temp->next = here->relStruct->deltaVthList ;
+//    here->relStruct->deltaVthList = temp ;
+//    printf ("QUI\n\n") ;
+
     if (!stress_or_recovery)
     {
         here->relStruct->t_star = pow ((exp (sqrt (here->relStruct->deltaVth / (Nt * A))) - 1), (1 / relmodel->RELMODELbeta)) * relmodel->RELMODELtau_0 ;
diff --git a/src/spicelib/devices/relmodel/relmodeldefs.h b/src/spicelib/devices/relmodel/relmodeldefs.h
index a32ca9f59..38d636452 100644
--- a/src/spicelib/devices/relmodel/relmodeldefs.h
+++ b/src/spicelib/devices/relmodel/relmodeldefs.h
@@ -11,11 +11,18 @@ Author: Francesco Lannutti - July 2015
 #include "ngspice/complex.h"
 #include "ngspice/noisedef.h"
 
+typedef struct sRELMODELrelList {
+    double time ;
+    double deltaVth ;
+    struct sRELMODELrelList *next ;
+} RELMODELrelList ;
+
 typedef struct sRELMODELrelStruct {
     double time ;
     double deltaVth ;
     double t_star ;
     int IsON ;
+    RELMODELrelList *deltaVthList ;
 } RELMODELrelStruct ;
 
 typedef struct sRELMODELmodel