diff --git a/data_analysis/diode.py b/data_analysis/diode.py
index 4c4a55d256f8377dc39176f38def027730516d6b..5d451047e338236c44cdb1356cb5af498da656cb 100644
--- a/data_analysis/diode.py
+++ b/data_analysis/diode.py
@@ -1,17 +1,119 @@
 import numpy as np
 import semiconductor as sc
+from scipy.optimize import fsolve
+def ideal_diode(Vbias,Is,n=1, temp=300):
+  
+    """ Function to calculate the current in an ideal diode 
+        
+        args:
+            - Vbias (scalar or sequence): the bias voltage of the diode
+            - Is (scalar): the saturation current of the diode
+            - n (scalar): the ideality factor of the diode, 1 for radiative recombination, 2 for SRH recombination
+            - temp (scalar): the temperature
+                
+        return:
+            - a scalar or sequence with same dimension as Vbias
+    """ 
+    
+    kB = 1.38e-23 # J/K
+    q = 1.602e-19 # C
+
+    return Is * (np.exp(q * Vbias / (n * kB * temp)) -1 )
+
+
+def two_diodes(Vbias,Is1,Is2,n1=1,n2=2, temp=300):
+    """ Function to calculate the current for a two diodes model 
+        
+        args:
+            - Vbias (scalar or sequence): the bias voltage of the diode
+            - Is1 (scalar): the saturation current of the first diode
+            - Is2 (scalar): the saturation current of the second diode
+            - n1 (scalar): the ideality factor of the first diode, 1 for radiative recombination, 2 for SRH recombination
+            - n1 (scalar): the ideality factor of the second diode, 1 for radiative recombination, 2 for SRH recombination
+            - temp (scalar): the temperature
+                
+        return:
+            - a scalar or sequence with same dimension as Vbias
+    """ 
+    kB = 1.38e-23 # J/K
+    q = 1.602e-19 # C
+
+    return Is1 * (np.exp( q * Vbias / (n1 * kB * temp)) -1 ) + Is2 * (np.exp( q * Vbias / (n2 * kB * temp)) -1 )
+
+def two_diodes_with_resistances(Vbias,Is1,Is2,n1=1,n2=2, temp=300, Rs=0, Rsh=float("inf")):
+    """ Function to calculate the current for a two diodes model by taking into account the series and shunt resistances
+        
+        args:
+            - Vbias (scalar or sequence): the bias voltage of the diode
+            - Is1 (scalar): the saturation current of the first diode
+            - Is2 (scalar): the saturation current of the second diode
+            - n1 (scalar): the ideality factor of the first diode, 1 for radiative recombination, 2 for SRH recombination
+            - n1 (scalar): the ideality factor of the second diode, 1 for radiative recombination, 2 for SRH recombination
+            - temp (scalar): the temperature
+            - Rs (scalar): the serie resistance 
+            - Rsh (scalar): the shunt resistance 
+                
+        return:
+            - a scalar or sequence with same dimension as Vbias
+    """ 
+    
+    kB = 1.38e-23 # J/K
+    q = 1.602e-19 # C
+    
+
+    if isinstance(Vbias, (int,float)):
+        I = fsolve(lambda x:Is1 * (np.exp( q * (Vbias - Rs * x) / (n1 * kB * temp)) -1 ) + Is2 * (np.exp( q * (Vbias - Rs * x) / (n2 * kB * temp)) -1 ) 
+           + (Vbias - Rs * x ) / Rsh - x,x0=two_diodes(Vbias,Is1,Is2,n1,n2, temp))
+    else:
+        I=np.zeros(len(Vbias))
+        i=0
+        for v in Vbias:
+            I[i] = fsolve(lambda x:Is1 * (np.exp( q * (v - Rs * x) / (n1 * kB * temp)) -1 ) + Is2 * (np.exp( q * (v - Rs * x) / (n2 * kB * temp)) -1 ) 
+               + (v - Rs * x ) / Rsh - x,x0=two_diodes(v,Is1,Is2,n1,n2, temp))
+    return I
 
 def depletion_length(doping_in, doping_out, Vbias=0,temp=300):
+    """ Function to calculate the depletion length in a pn junction
+        
+        args:
+            - doping_in (scalar): the doping in the region for which the depletion length has to be calculated
+            - doping_out (scalar): the doping in the adjacent region for which the depletion length has to be calculated
+            - Vbias (scalar): the bias voltage of the pn junction
+            - temp (scalar): the temperature
+                
+        return:
+            - a scalar with the depletion length calculated in one region
+    """ 
+    
+    kB = 1.38e-23 # J/K
+    q = 1.602e-19 # C
+    epsilon_0 = 8.8542e-12 # F/m
+    epsilon_si = 11.7
+    
     phi_0 = kB * temp / q * np.log(doping_in * doping_out / sc.intrinsic_concentration(temp)**2 )
     return np.sqrt(2 * epsilon_si * epsilon_0 / q * doping_out / doping_in / (doping_in + doping_out) * (phi_0 - Vbias))
 
 
-def j_srh(mu_n,mu_p,ND,NA,tau=1e-7,temp=300,Vbias=0):
-      
+def j_srh(Vbias,ND,NA,tau=1e-7,temp=300):
+    """ Function to calculate the Shockley-Read-Hall contribution to the current density in a pn junction
+        
+        args:
+            - Vbias (scalar): the bias voltage of the pn junction
+            - ND (scalar): the donor doping concentration in the n region 
+            - NA (scalar): the acceptor doping concentraion in the p region 
+            - tau (scalar): the global lifetime associated to the SRH mechanism
+            - temp (scalar): the temperature
+
+        return:
+            - a scalar with the SRH current density calculated 
+    """  
+    kB = 1.38e-23 # J/K
+    q = 1.602e-19 # C
+
     ld_n = depletion_length(ND, NA, Vbias)
     ld_p = depletion_length(NA, ND, Vbias)
     ni = sc.intrinsic_concentration(temp)
-    x=(ld_p + ld_n) # approximation by considering only the depletion region without diffusion mechanism
+    x=(ld_p + ld_n) # approximation by considering only the depletion region without diffusion mechanism, gives an upper limit as the effective length is always below
     
     coeff_SRH = q * ni * x / (2 * tau)
     
@@ -19,8 +121,29 @@ def j_srh(mu_n,mu_p,ND,NA,tau=1e-7,temp=300,Vbias=0):
 
 
 
-def j_radial(mu_n,mu_p,tau_n,tau_p,ND,NA,ln,lp,temp=300,Vbias=0):
-    
+def j_radiative(Vbias,mu_n,mu_p,tau_n,tau_p,ND,NA,ln,lp,temp=300):
+    """ Function to calculate the radial contribution to the current density in a silicon pn junction
+        
+        args:
+            - Vbias (scalar): the bias voltage of the pn junction
+            - mu_n (scalar): the mobility of the electrons
+            - mu_p (scalar): the mobility of the holes
+            - tau_n (scalar): the lifetime of the electrons
+            - tau_p (scalar): the lifetime of the holes
+            - ND (scalar): the donor doping concentration in the n region 
+            - NA (scalar): the acceptor doping concentraion in the p region 
+            - ln (scalar): the length of the cathode (n-doped region)
+            - lp (scalar): the length of the anode (p-doped region)
+            - temp (scalar): the temperature
+
+        return:
+            - a scalar with the radiative current density calculated 
+    """  
+    b_rad = 4.76e-15 # cm3/s - low-impurity value entre 1 et 10
+
+    kB = 1.38e-23 # J/K
+    q = 1.602e-19 # C
+
     Dn = kB * temp / q * mu_n
     Dp = kB * temp / q * mu_p
     
@@ -38,15 +161,3 @@ def j_radial(mu_n,mu_p,tau_n,tau_p,ND,NA,ln,lp,temp=300,Vbias=0):
     
     
     return q * ( coeff_radial * (np.exp(q * Vbias/ ( kB * temp)) - 1 ) )
-
-
-#################################################################
-#
-#                          CONSTANTS
-#
-#################################################################
-b_rad = 4.76e-15 # cm3/s - low-impurity value entre 1 et 10
-kB = 1.38e-23 # J/K
-q = 1.602e-19 # C
-epsilon_0 = 8.8542e-12 # F/m
-epsilon_si = 11.7
\ No newline at end of file
diff --git a/data_analysis/semiconductor.py b/data_analysis/semiconductor.py
index 93121b6d74657d46e3d9f9224109c4def5ca5cfe..53570f554775aa8ebe0509c32a313dabcb866f5c 100644
--- a/data_analysis/semiconductor.py
+++ b/data_analysis/semiconductor.py
@@ -1,5 +1,6 @@
 import numpy as np
 import mechanics as mec
+import scipy.interpolate as interp
 
 def mobility_impurity(carrier="n",Ni=1e15,temp=300, dopant="phosphorus"):
     """ Function to calculate the silicon mobility according to Masetti relation (1983)
@@ -166,3 +167,30 @@ def piezoresitivity_strain(strain_tensor, pi11, pi12, pi44):
     stress_tensor= mec.stress_from_strain(strain_tensor)
     
     return piezo_tensor @ stress_tensor
+
+def piezo_ratio_temperature(temp,carrier="n"):
+    """ Function to calculate the correcting ratio due to temperature variations for the piezoresistive coefficients. The data are taken from Kanda (1982) between 200 K and 400K.
+        The correction is refered to the 300K coefficients. The new piezoresistive coefficients can be calculated by piezo_coefficient(temp) = ratio * piezo_coefficient(300K)
+        
+        args:
+            - temp (scalar or sequence): the temperatures for which the correting ratio has to be calculated
+            - carrier (string): "n" for electrons, "p" for holes
+                
+        return:
+            - a scalar or sequence with the same dimension as temp with the correcting ratio
+    """  
+
+    temp_kanda1982=np.linspace(125, -75,9)+273
+    dpi_kanda1982_n=np.array([0.7547318611987381,0.8067823343848579,0.8611987381703468,0.9298107255520504,1.0007886435331228,1.0977917981072554,1.2113564668769716,1.3438485804416402,1.5165615141955835])
+    dpi_kanda1982_p=np.array([0.7523342983957861,0.8037360206464299,0.8598127702627387,0.922894558591104,0.995324258126102,1.0934548187864221,1.205601887444947,1.3411219494071012,1.5093457676819353])
+    
+    fn=interp.interp1d(temp_kanda1982,dpi_kanda1982_n,kind="linear",fill_value="extrapolate" )
+    fp=interp.interp1d(temp_kanda1982,dpi_kanda1982_p,kind="linear",fill_value="extrapolate" )
+
+    if carrier=="n":
+        return fn(temp)
+    elif carrier =="p":
+        return fp(temp)
+    else:
+        return 0
+    
diff --git a/data_analysis/transistor.py b/data_analysis/transistor.py
index aa6d59a2c87d934b8fca9b4ec0be426717ad6eed..3976e5b4b7ba8af482052f5092547bb584e50797 100644
--- a/data_analysis/transistor.py
+++ b/data_analysis/transistor.py
@@ -2,10 +2,25 @@ import data_processing as proc
 import numpy as np
 import scipy.interpolate as interp
 
-def threshold_voltage_extraction(i,v,smoothing=True,accuracy=6): # second-derivative method
-    v_step=np.mean(v[1:]-v[:-1])
+def threshold_voltage_extraction(ids,vgs,accuracy=4): 
+    """ Function to find the threshold voltage of a transistor defined as the voltage where the transconductance gm is maximum, i.e., when the first derivative of the current is maximum or when the second derivative is zero.
+        The method is based on the second-derivative method from Wong HS, White MH, Krutsick TJ, Booth RV. "Modeling of transconductance degradation and extraction of threshold voltage in thin oxide MOSFET’s." Solid-State Electron (1987).
+        The transistor should be put in linear region (VDS <= VGS-Vth), i.e., with low VDS bias.
+        The implementation is based on central finite difference to limit the impact of noise for the determination of the threshold voltage.
+        The numerical derivative is done on accuracy + 1 points, but accuracy should be even
+        
+        args:
+            - ids (scalar): the source to drain current of the transistor (nmos or pmos)
+            - vgs (scalar): the voltage between the gate and the source of the transistor (nmos or pmos) 
+            - accuracy (integer): related to the number of points taken for the finite difference. For central finite difference, the number of points for a first order derivation is accuracy + 1 but accuracy should be even to guarantee a centered interval.
+            
+        return:
+            - a tuple of two scalars with the threshold voltage found and the transconductance value found at this voltage
+    """  
+
+    v_step=np.mean(vgs[1:]-vgs[:-1])
     
-    gm,v_gm=proc.finite_difference(i, v_step,order=2,accuracy=accuracy)
+    gm,v_gm=proc.finite_difference(ids, v_step,order=1,accuracy=accuracy)
     
     f = interp.InterpolatedUnivariateSpline(v_gm, gm, k=4)
     cr_pts = f.derivative().roots()
@@ -17,7 +32,27 @@ def threshold_voltage_extraction(i,v,smoothing=True,accuracy=6): # second-deriva
     
     return vth, gm_max
     
-def mos_transistor_current(vg,vd,vs=0,vth=0.8,k=1e-4,vea=50,mos_type="nmos",early=True):
+def mos_transistor_current(vg,vd,vs=0,vth=0.8,k=1e-4,vea=50,theta=0,mos_type="nmos",early=True):
+    """ Function to calculate the current in a MOS transistor:
+            cutoff (vgs < vth) : ids = 0
+            triode (vgs - vth >= vds) : ids = k * ( (vgs-vth) * vds - vds**2/2) * (1+(vds)/vea)
+            saturation (vgs - vth < vds)  : ids = k/2 * ( (vgs-vth)**2) * (1+(vds)/vea)
+        with k = "mobility" * "oxide capacitance" * "width" / "length" where the oxide capacitance is given by the ratio between the permittivity and the oxide thickness, i.e., Cox = epsilon_r * epsilon_0 / t_ox
+        The mobility degradation is taken into account by replacing k by k / (1 + theta * (vgs - vth)) in the current equation.
+        
+        args:
+            - vg (scalar or sequence): gate voltage of the transistor. Maximum one parameter among vg, vd and vs can be a sequence.
+            - vd (scala or sequence): drain voltage of the transistor. Maximum one parameter among vg, vd and vs can be a sequence.
+            - vs (scalar or sequence): source voltage of the transistor. Maximum one parameter among vg, vd and vs can be a sequence.
+            - vth (scalar): threshold voltage of the transistor
+            - k (scalar): the gain factor of the transistor
+            - vea (scalar) : early voltage of the transistor
+            - theta (scalar) : mobility degradation factor
+            
+        return:
+            - a scalar or sequence with the value of current calculated
+    """  
+    
     if mos_type=="nmos":
         vth=abs(vth)
         vds=vd-vs
@@ -29,18 +64,20 @@ def mos_transistor_current(vg,vd,vs=0,vth=0.8,k=1e-4,vea=50,mos_type="nmos",earl
         
     if vgs<vth:
         mode="cutoff"
-    elif vgs>=vds:
+    elif vgs-vth>=vds:
         mode="saturation"
-    elif vg-vs<vds:
+    elif vgs-vth<vds:
         mode="triode"
 
+    k_theta=k/(1+theta*(vgs-vth))
+
     if mode=="cutoff":
         jd=0
     elif mode=="saturation":
-        jd=k/2*(vgs-vth)**2
+        jd=k_theta/2*(vgs-vth)**2
     
     elif mode=="triode":
-        jd=k*((vgs-vth)*(vds)-(vds)**2/2)
+        jd=k_theta*((vgs-vth)*(vds)-(vds)**2/2)
 
     if early:
         return jd*(1+(vds)/vea)