diff --git a/GeoTools/flatten_stack.py b/GeoTools/flatten_stack.py
new file mode 100755
index 0000000..2e4e95f
--- /dev/null
+++ b/GeoTools/flatten_stack.py
@@ -0,0 +1,397 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Model coefficient estimation and removal for wrapped interferograms
+
+Usage:
+    flatten_stack.py <phase> <phase_filt> <model> [--outfile=<outfile>] [--nreg=<float>] [--thresh_amp=<float>] [--thresh_cohreg=<float>] [--thresh_model=<float>] [--thresh_std_model=<float>] [--thresh_min_pixel=<int>] [--cyclmax=<float>] [--plot=<yes/no>] [--plot_reg=<yes/no>]
+    flatten_stack.py -h | --help
+    flatten_stack.py estimate <phase_filt> <model> [--outfile=<outfile>] [--nreg=<nreg>] [--thresh_amp=<float>] [--thresh_cohreg=<float>] [--thresh_model=<float>] [--thresh_std_model=<float>] [--thresh_min_pixel=<int>] [--cyclmax=<float>] [--plot_reg=<yes/no>]
+    flatten_stack.py add <unwrapped> <model> --coeff=<coeff> [--outfile=<outfile>]
+    flatten_stack.py remove <phase> <model> --coeff=<coeff> [--outfile=<outfile>]
+
+Options:
+    -h --help           Display command details
+    phase               None filtered interferogram (.int)       
+    phase_filt          Filtered interferogram (.int)
+    model               Model for the propotionality estimation (.r4)
+    outfile             Outfile name
+    nreg                Number of region in range, put 1 for 1 global region (default: 28)
+    thresh_amp          Minimal value of amplitude for pixel selection (default: 0.1)
+    thresh_cohreg       Minimal value of coherence for region selection (default: 0.2)
+    thresh_model        Minimal value of model for pixel selection (default: 0.4)
+    thresh_std_model    Minimal standard deviation within a window for estimation (default: 0.3)
+    thresh_min_pixel    Minimal percent of pixel per windows for estimation (default: 10)
+    cyclmax             Maximum number of phase cycle (default: 3)
+    plot                plot intermediate results (default: no)
+    plot_reg            plot each phase model relation per region (default: no)
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import least_squares
+from osgeo import gdal
+from numba import jit
+import docopt
+import sys
+import os
+from scipy.optimize import minimize
+gdal.UseExceptions()
+
+def open_gdal(file, band=1, supp_ndv=None, complex=False):
+    """
+    Use GDAL to open band as real value
+    """
+    print('-----')
+    print(file)
+    if not os.path.isfile(file):
+        raise FileNotFoundError('File does not exists: {}'.format(file))
+    ds = gdal.Open(file)
+
+    print('dims', ds.RasterXSize, ds.RasterYSize, ds.RasterCount)
+    band = ds.GetRasterBand(band)
+    ndv = band.GetNoDataValue()
+    data = band.ReadAsArray()
+    Xsize = ds.RasterXSize
+    Ysize = ds.RasterYSize
+    if complex:
+        amp = np.absolute(data)
+        phi = np.angle(data)
+        data = [amp, phi]
+        ndv = 0.0
+        if ndv is not None and ndv != np.nan:
+            amp[amp==ndv] = np.nan
+            phi[phi==ndv] = np.nan
+        if supp_ndv is not None and supp_ndv != np.nan:
+            amp[amp==supp_ndv] = np.nan
+            phi[phi==supp_ndv] = np.nan
+    else:
+        if ndv is not None and ndv != np.nan:
+            data[data==ndv] = np.nan
+        if supp_ndv is not None and supp_ndv != np.nan:
+            data[data==supp_ndv] = np.nan
+    return data, Xsize, Ysize
+
+def save_gdal(outfile, data, template=None, band=1, alt_ndv=None, roi_pac=True):
+    """
+    Create a new raster with same dims as the template
+    """
+    if outfile == 'nomodel':
+        pos_rlks = template.rfind('rlks')
+        pos_underscore = template.rfind('_', 0, pos_rlks)
+        outfile = template[:pos_underscore+1] + outfile + '_' + template[pos_underscore+1:]
+    if template is not None:
+        ds_template = gdal.Open(template)
+        if ds_template is None:
+            raise RuntimeError(f"Imposible to open the template {template}")
+        driver_name = "ROI_PAC" if roi_pac else "GTiff"
+        driver = gdal.GetDriverByName(driver_name)
+        ds = driver.CreateCopy(outfile, ds_template)
+    else:
+        ds = gdal.OpenEx(outfile, gdal.GA_Update)
+    band = ds.GetRasterBand(band)
+    ndv = band.GetNoDataValue()
+    if alt_ndv is not None:
+        ndv = alt_ndv
+        band.SetNoDataValue(ndv)
+    if ndv is not None:
+        data[np.isnan(data)] = ndv
+    band.WriteArray(data)
+    ds.FlushCache()
+    print('save : ', outfile)
+
+
+def former_flatten_stack(phase, phase_filt, model, **kwargs):
+    """
+    Remove a pattern from a wrapped signal, as in flatten_stack.f from NSBAS
+    """
+
+    # nregs = [1, 5, 10, 20, 40, 80, 160, 200]
+    # thresh_model = [-0.6, -0.3, 0, 0.3, 0.4, 0.6]
+    # #nregs = [1, 200]
+    # #thresh_model = [0, 0.6]
+    # results_mediane = {t: [] for t in thresh_model}
+    # for nreg in nregs:
+    #     kwargs["nreg"] = nreg
+    #     for thresh in thresh_model:
+    #         kwargs["thresh_model"] = thresh
+    #         median = estimate_coeff(phase_filt, model, **kwargs)
+    #         results_mediane[thresh].append(median)
+    # colors = plt.cm.viridis(np.linspace(0, 1, len(thresh_model)))
+    # for color, thresh in zip(colors, thresh_model):
+    #     plt.plot(nregs, results_mediane[thresh], '-o', color=color,
+    #             label=f'médiane (thresh={thresh})', linewidth=2)
+        
+    # plt.xlabel('nreg', fontsize=12)
+    # plt.ylabel('Coefficient', fontsize=12)
+    # plt.title('Évolution des coefficients selon nreg et thresh_model', fontsize=13)
+    # plt.grid(alpha=0.3)
+    # plt.legend(frameon=False)
+    # plt.tight_layout()
+    # plt.show()
+
+    coeff = estimate_coeff(phase_filt, model, **kwargs)
+    nomodel_phifilt = remove_model(phase_filt[1], model, coeff)
+    nomodel_phi = remove_model(phase[1], model, coeff)
+
+    nomodel_phifilt[np.isnan(nomodel_phifilt)] = 0.0
+    nomodel_phi[np.isnan(nomodel_phi)] = 0.0
+    phase_filt[0][np.isnan(phase_filt[0])] = 0.0
+    phase[0][np.isnan(phase[0])] = 0.0
+
+    complex_nomodel_phifilt = phase_filt[0] * np.exp(1j*nomodel_phifilt)
+    complex_nomodel_phi = phase[0] * np.exp(1j*nomodel_phi)
+
+    if kwargs['plot'] == 'yes':
+        fig, axes = plt.subplots(1, 3, figsize=(12, 5))
+
+        im0 = axes[0].imshow(phase_filt[1], cmap='twilight')
+        axes[0].set_title('phase')
+        plt.colorbar(im0, ax=axes[0])
+
+        im1 = axes[1].imshow(nomodel_phifilt, cmap='twilight')
+        axes[1].set_title('corrected_phase')
+        plt.colorbar(im1, ax=axes[1])
+        
+        im1 = axes[2].imshow(model, cmap='RdBu_r',vmin=np.nanpercentile(model, 2), vmax=np.nanpercentile(model, 98))
+        axes[2].set_title('Model')
+        plt.colorbar(im1, ax=axes[2])
+        plt.show()
+    
+    print('')
+    save_gdal(kwargs['outfile'], complex_nomodel_phifilt, template=arguments["<phase_filt>"])
+    save_gdal(kwargs['outfile'], complex_nomodel_phi, template=arguments["<phase>"])
+
+    dates = next(part for part in kwargs['outfile'].split('_') if '-' in part)
+    with open(f"{dates}.stack", 'w') as f:
+        f.write(f"{coeff:.10f}\n")
+
+def estimate_coeff(phase_filt, model, **kwargs) -> float:
+    """
+    Estimate the proportionnality coefficient between the model and the wrapped phase values
+    minimizing the complex product (maximizing the coherence)
+    """
+    amp_filt = np.copy(phase_filt[0])
+    phi_filt = np.copy(phase_filt[1])
+    model_copy = np.copy(model)
+    
+    index = np.isnan(amp_filt)| np.isnan(phi_filt) | np.isnan(model) | (amp_filt<kwargs['thresh_amp']) | (model<kwargs['thresh_model'])
+
+    amp_filt[index] = np.nan
+    phi_filt[index] = np.nan
+    model_copy[index] = np.nan
+
+    ## Creation of sub-area
+    nreg = kwargs["nreg"]
+    print('---')
+    print('Number of region :', nreg)
+    if nreg != 1:
+        amp_filt = reg_pad(amp_filt, nreg)
+        phi_filt = reg_pad(phi_filt, nreg)
+        model_split = reg_pad(model_copy, nreg)
+    else:
+        amp_filt = np.array([[amp_filt]])
+        phi_filt = np.array([[phi_filt]])
+        model_split = np.array([[model_copy]])
+        kwargs['thresh_cohreg'] = 0.           
+
+    # Initialisation
+    cyclmax = kwargs['cyclmax']
+    list_coeff = np.full(shape=(model_split.shape[0], model_split.shape[1]), fill_value=np.nan)
+    list_coh = np.full(shape=(model_split.shape[0], model_split.shape[1]), fill_value=np.nan)
+
+    print('\nStarting maximisation')
+    for i in range(model_split.shape[0]):     
+        for j in range(model_split.shape[1]):
+            # for each region
+            phi_filt_region = phi_filt[i,j]
+            model_region = model_split[i,j]
+
+            if np.count_nonzero(~np.isnan(phi_filt_region)) < kwargs['thresh_min_pixel']:
+                # print("skip not enough px", i, j)
+                continue
+            elif np.nanstd(model_region) < kwargs['thresh_std_model']:
+                # print("skip not enough std", i, j)
+                continue
+
+            res = minimize(misfit, x0=[0.0], args=(phi_filt_region, model_region), 
+                          bounds=[(-cyclmax, cyclmax) ] )
+            coeff = res.x[0]
+
+            coh = np.abs(np.nanmean(np.exp(1j * phi_filt_region) * np.exp(-1j * (coeff * model_region ))))
+            
+            if coh > kwargs['thresh_cohreg']:
+                list_coeff[i,j] = coeff
+                list_coh[i,j] = coh 
+
+            if kwargs['plot_reg'] == 'yes':
+                plt.figure(figsize=(5, 3))
+                plt.scatter(np.angle(np.exp(1j*model_region)), np.angle(np.exp(1j*phi_filt_region)), s=0.1, c='k', alpha=0.05, label="Data", rasterized=True)
+                
+                x_fit = np.linspace(np.nanmin(np.angle(np.exp(1j*model_region))), np.nanmax(np.angle(np.exp(1j*model_region))), 100)
+                y_fit = np.angle(np.exp(1j*coeff*x_fit))
+
+                x = np.angle(np.exp(1j*model_region))
+                y = np.angle(np.exp(1j*phi_filt_region))
+
+                y_fit2 = np.angle(np.exp(1j*coeff*model_region))
+                cst =  np.nanmean(y - y_fit2)
+
+                plt.plot(x_fit, y_fit + cst, 'r-', lw=2, label=f"Fit: coeff={list_coeff[i,j]:.3f}, coh={list_coh[i,j]:.3f}")
+                plt.xlabel("Model")
+                plt.ylabel("Phase filtered")
+                plt.legend()
+                plt.title(f"Window ({i}, {j})")
+                plt.grid(True, ls='--', alpha=0.4)
+                plt.tight_layout()
+                plt.show()
+    
+    med = np.nanmedian(list_coeff)
+    print(f'--> coeff median: {med:.3f}')
+
+    if kwargs['plot'] == 'yes':
+        fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 7))
+        im1 = ax1.imshow(list_coeff, cmap='RdBu_r',
+                    vmax=cyclmax,
+                    vmin=-cyclmax)
+        ax1.set_title('Coefficients')
+        plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)
+        
+        im2 = ax2.imshow(list_coh, cmap='viridis')
+        ax2.set_title('Coherence')
+        plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04)
+
+        im3 = ax3.imshow(model_copy, cmap='RdBu_r', vmin=np.nanpercentile(model, 2), vmax=np.nanpercentile(model, 98))
+        ax3.set_title('Model')
+        plt.colorbar(im3, ax=ax3, fraction=0.046, pad=0.04)
+        plt.tight_layout()
+        plt.show()
+
+    return med
+
+def misfit(params, phi_filt_region, model_region):
+    c = params  
+    coh = np.nanmean(np.exp(1j * phi_filt_region) *
+            np.exp(-1j * (c * model_region )))
+    return 1./np.abs(coh)
+    
+def reg_pad(array, block_nb_x):
+    """
+    Regionalize an array from a given block number in the x dimension.
+    The squared block size is determined from the x dimension, by padding the array
+    to ensure nx % block_size = 0.
+    The number of blocks in the y dimension is determined by conserving the block size 
+    constant between x and y axis.
+    """
+    ny, nx = array.shape
+    # taille de bloc à nombre de blocs fixé
+    dxy = nx // block_nb_x + 1
+    pad_x = dxy * block_nb_x - nx
+
+    # Nombre de blocs à taille de bloc fixée
+    nby = ny // dxy + 1
+    pad_y = dxy * nby - ny
+    
+    # Pad to ensure that it can perfectly be divided into blocks
+    pad_array = np.pad(array, ((0, pad_y), (0, pad_x)), constant_values=np.nan)
+    # Cut the different blocks into different sub-arrays
+    y_blocks = np.array(np.array_split(pad_array, nby, axis=0))
+    blocks = np.array([np.array_split(b, block_nb_x, axis=1) for b in y_blocks])
+    
+    return blocks
+
+#@jit(nopython=True)
+def remove_model(phase, model, coeff):
+    """
+    Remove a model from wrapped values with a given proportionnality coefficient
+    """
+    corrected_phase = np.angle(np.exp(1j*phase) * np.exp(-1j*(coeff*model)))
+    return corrected_phase
+
+@jit
+def weighted_median(values, weights):
+    mask = ~np.isnan(values)
+    values = values[mask] ; weights = weights[mask]
+
+    sorter = np.argsort(values)
+    values = values[sorter]
+    weights = weights[sorter]
+
+    cumulative_weight = np.cumsum(weights)
+
+    return values[cumulative_weight >= cumulative_weight[-1] / 2.0][0]
+
+def add_model(phase, model, coeff):
+    """
+    Add a model to values given a proportionnality coefficient
+    """
+    return (phase + coeff*model)
+
+def arg2value(value, conversion=None, default=None):
+    """Convert a string argument if exists otherwise use default value"""
+    if value is None:
+        return default
+    elif conversion is not None:
+        return conversion(value)
+    return value
+
+
+if __name__ == '__main__':
+    arguments = docopt.docopt(__doc__)
+    print(arguments)
+    
+    # Additionnal parameters
+    outfile = arg2value(arguments["--outfile"], str, 'nomodel')
+    nreg = arg2value(arguments["--nreg"], int, 28)
+    thresh_amp = arg2value(arguments["--thresh_amp"], float, 0.1)
+    thresh_cohreg = arg2value(arguments["--thresh_cohreg"], float, 0.2)
+    thresh_model = arg2value(arguments["--thresh_model"], float, 0.4)
+    thresh_std_model = arg2value(arguments["--thresh_std_model"], float, 0.3)
+    thresh_min_pixel = arg2value(arguments["--thresh_min_pixel"], float, 10) # nb px (TODO en %)
+    cyclmax = arg2value(arguments["--cyclmax"], float, 3.)
+    plot = arg2value(arguments["--plot"], str, 'no')
+    plot_reg = arg2value(arguments["--plot_reg"], str, 'no')
+
+    if arguments["estimate"]:
+        # Find the optimal coefficient of proportionality between the phase and the model
+        phase, pXsize, pYsize = open_gdal(arguments["<phase>"], band=1, complex=True)
+        phase_filt, pfXsize, pfYsize = open_gdal(arguments["<phase_filt>"], band=1, complex=True)
+        model, mXsize, mYsize = open_gdal(arguments["<model>"])
+        if not (pXsize == pfXsize == mXsize and pYsize == pfYsize == mYsize):
+            sys.exit("Error: input rasters (phase, phase_filt, model) do not have the same dimensions.")
+        coeff = estimate_coeff(phase_filt, model, nreg=nreg, thresh_amp=thresh_amp, thresh_cohreg=thresh_cohreg,thresh_model=thresh_model, 
+                       thresh_std_model=thresh_std_model, thresh_min_pixel=thresh_min_pixel, cyclmax=cyclmax, plot=plot, plot_reg=plot_reg)
+        print(coeff)
+    elif arguments["add"]:
+        # Add the model to an unwrapped file
+        unwrapped, unwXsize, unwYsize = open_gdal(arguments["<unwrapped>"], band=2)
+        model, mXsize, mYsize = open_gdal(arguments["<model>"])
+        if not (unwXsize == mXsize and unwYsize == mYsize):
+            sys.exit("Error: input rasters (unw, model) do not have the same dimensions.")
+        coeff = float(arguments["--coeff"])
+        outfile = arguments["--outfile"]
+        unw_plus_model = add_model(unwrapped, model, coeff)
+        save_gdal(outfile, unw_plus_model, template=arguments["<unwrapped>"], band=2)
+    elif arguments["remove"]:
+        # Remove the model from wrapped signal
+        phase, pXsize, pYsize = open_gdal(arguments["<phase>"], band=1, complex=True)
+        model, mXsize, mYsize = open_gdal(arguments["<model>"])
+        if not (pXsize == mXsize and pYsize == mYsize):
+            sys.exit("Error: input rasters (phase, model) do not have the same dimensions.")
+        
+        coeff = float(arguments["<coeff>"])
+        outfile = arguments["--outfile"]
+        phase_minus_model = remove_model(phase.astype(np.float64), modelastype(np.float64), float(coeff))
+        save_gdal(outfile, phase, phase_minus_model, band=2)
+    else:
+        # Remove a pattern from a wrapped signal, as in flatten_stack.f from NSBAS
+        phase, pXsize, pYsize = open_gdal(arguments["<phase>"], complex=True)
+        phase_filt, pfXsize, pfYsize = open_gdal(arguments["<phase_filt>"], complex=True)
+        model, mXsize, mYsize = open_gdal(arguments["<model>"])
+
+        if not (pXsize == pfXsize == mXsize and pYsize == pfYsize == mYsize):
+            sys.exit("Error: input rasters (phase, phase_filt, model) do not have the same dimensions.")
+        phase_minus_model = former_flatten_stack(phase, phase_filt, model, 
+                                outfile=outfile, nreg=nreg, thresh_amp=thresh_amp, thresh_cohreg=thresh_cohreg, thresh_model=thresh_model, 
+                                thresh_std_model=thresh_std_model, thresh_min_pixel=thresh_min_pixel, cyclmax=cyclmax, plot=plot, plot_reg=plot_reg)
+
diff --git a/NSBAS-playground/scripts/nsb_filtflatunw.py b/NSBAS-playground/scripts/nsb_filtflatunw.py
index 46abfca..1270fc4 100755
--- a/NSBAS-playground/scripts/nsb_filtflatunw.py
+++ b/NSBAS-playground/scripts/nsb_filtflatunw.py
@@ -1186,6 +1186,12 @@ def flat_model(config,kk):
         copyrsc(inrsc,outrsc)
         copyrsc(inrsc,filtoutrsc)
 
+        #Parameters
+        nregion = int(1)
+        thresh_amp = 0.1
+        thresh_cohreg = 0.6
+        thresh_model = 0.4
+
         if force:
             rm(outfile); rm(param)
             try:
@@ -1196,8 +1202,8 @@ def flat_model(config,kk):
             do = checkoutfile(config,outfile)
             if do:
                 try:
-                    run("flatten_stack "+str(infile)+" "+str(filtfile)+" "+str(config.model)+" "+str(outfile)+" "+str(filtout)\
-                    +" "+str(config.thresh_amp_atmo)+" > log_flatmodel.txt")
+                    run(f"flatten_stack.py {infile} {filtfile} {config.model} --outfile={outfile} --nreg={nregion} --thresh_amp={thresh_amp} --thresh_cohreg={thresh_cohreg} --thresh_model={thresh_model}")
+                    #run("flatten_stack.py "+str(infile)+" "+str(filtfile)+" "+str(config.model)+" "+"--outfile="+str(outfile)+" "+"--nreg="+str(nregion)+ "")
                     # move param file into a file name independent of prefix and suffix
                     force_link(param,newparam)
                 except Exception as e:
@@ -1688,16 +1694,20 @@ def add_model_back(config,kk):
             outfile = config.stack.getfiltROI(kk) + '.unw'
             outrsc = outfile + '.rsc'
             copyrsc(unwrsc,outrsc)
-
+            print(param)
             if force:
                 rm(outfile)
             if config.model != None:
                 if path.exists(param) is True:
+                    with open(param, "r") as f:
+                        coeff = float(f.read())
                     do = checkoutfile(config,outfile)
                     if do:
                         try:
+                            print(coeff)
                             run("length.pl "+str(unwfile))
-                            run("unflatten_stack "+str(unwfile)+" "+str(outfile)+" "+str(config.model)+" "+str(param)+" >> log_flatmodel.txt")
+                            run(f"flatten_stack.py add "+str(unwfile)+" "+str(config.model)+" "+"--coeff="+str(coeff)+" "+"--outfile="+str(outfile))
+                            #run("unflatten_stack "+str(unwfile)+" "+str(outfile)+" "+str(config.model)+" "+str(param)+" >> log_flatmodel.txt")
                         except Exception as e:
                             logger.critical(e)
                             logger.critical("Unflatten model failed for int. {0} Failed!".format(unwfile))
diff --git a/TimeSeries/invers_disp_pixel.py b/TimeSeries/invers_disp_pixel.py
index 70c3226..eb200d1 100755
--- a/TimeSeries/invers_disp_pixel.py
+++ b/TimeSeries/invers_disp_pixel.py
@@ -12,6 +12,7 @@
 
 Usage: invers_disp_pixel.py --cols=<values> --lines=<values> [--cube=<path>] [--list_images=<path>] [--windowsize=<value>] [--windowrefsize=<value>]  [--lectfile=<path>] [--aps=<path>] \
 [--linear=<value>] [--steps=<value>] [--postseismic=<value>] [--seasonal=<yes/no>] [--seasonal_increase=<yes/no>] [--vector=<path>] [--info=<path>]\
+[--linear=<value>] [--acceleration=<value>] [--steps=<value>] [--postseismic=<value>] [--seasonal=<yes/no>] [--seasonal_increase=<yes/no>] [--vector=<path>] [--info=<path>]\
 [--semianual=<yes/no>] [--bianual=<yes/no>] [--degreeday=<values>] [--bperp=<yes/no>] [--imref=<value>] [--cond=<value>] [--slowslip=<value>] [--ineq=<value>] \
 [--name=<value>] [--scale=<value>] [--plot=<yes/no>] [<iref>] [<jref>] [--bounds=<value>] [--dateslim=<values>] [--plot_dateslim=<values>] [--color=<value>] [--fillstyle=<value>] [--ndatasets=<nb_data_sets>] 
 
@@ -24,11 +25,12 @@
 --cube PATH             Path to displacement file [default: depl_cumul_flat]
 --list_images PATH      Path to list images file made of 5 columns containing for each images 1) number 2) Doppler freq (not read) 3) date in YYYYMMDD format 4) numerical date 5) perpendicular baseline [default: images_retenues]
 --windowsize VALUE      Number of pixels around the pixel defining the window [default: 0]
---windowrefsize VALUE      Number of pixels around the referenced pixel defining the window [default: windowsize]
+--windowrefsize VALUE   Number of pixels around the referenced pixel defining the window [default: windowsize]
 --lectfile PATH         Path to the lect.in file (output of invers_pixel) [default: lect.in]
 --aps PATH              Path to the APS file giving the error associated to each dates [default: No weigthing]
---linear PATH     Add a linear function to the inversion
---steps PATH        Add heaviside functions to the inversion .Indicate steps time (e.g 2004.,2006.)
+--linear PATH           Add a linear function to the inversion
+--acceleration=PATH     Add an acceleration fonction in the inversion [default: no]
+--steps PATH            Add heaviside functions to the inversion .Indicate steps time (e.g 2004.,2006.)
 --postseismic PATH      Add logarithmic transients to each steps step. Indicate characteristic time of the log function, must be a serie of values of the same lenght than steps (e.g 1.,1.). To not associate postseismic function to a given steps step, put None (e.g None,1.) 
 --slowslip   VALUE      Add slow-slip function in the inversion (as defined by Larson et al., 2004). Indicate median and characteristic time of the events (e.g. 2004.,1,2006,0.5) [default: None] 
 --seasonal PATH         If yes, add seasonal terms in the inversion
@@ -38,15 +40,19 @@
 --bianual PATH          If yes, add bianual  terms in the inversion
 --vector PATH           Path to the vector text files containing a value for each dates [default: None]
 --info PATH             Path to extra file in r4 or tif format to plot is value on the selected pixel, e.g. aspect [default: None].
+<<<<<<< HEAD
 --bperp PATH              If yes, add term proportional to the perpendicular baseline in the inversion
+=======
+--bperp PATH            If yes, add term proportional to the perpendicular baseline in the inversion
+>>>>>>> 06ce944 (Add an acceleration fonction in the inversion [default: no])
 --imref VALUE           Reference image number [default: 1]
 --cond VALUE            Condition value for optimization: Singular value smaller than cond are considered zero [default: 1.e-6]
 --ineq VALUE            If yes, add inequality constrained in the inversion: use least square result to iterate the inversion. Force postseismic to be the  same sign than steps [default: no].       
 --name Value            Name output figures [default: None] 
---scale                Scaling value between input data and desired output [default: 1]
+--scale                 Scaling value between input data and desired output [default: 1]
 --plot                  Display results [default: yes]            
-iref                  colum numbers of the reference pixel [default: None] 
-jref                  lign number of the reference pixel [default: None]
+iref                    colum numbers of the reference pixel [default: None] 
+jref                    lign number of the reference pixel [default: None]
 --bounds                yMin,yMax time series plots 
 --dateslim              Datemin,Datemax time series  
 --plot_dateslim         Datemin,Datemax time series for plot only 
@@ -98,7 +104,6 @@
 # docopt (command line parser)
 import docopt
 
-
 ########################################################################
 # Define basis functions
 ########################################################################
@@ -152,6 +157,15 @@ def __init__(self,name,reduction,date):
     def g(self,t):
         func=(t-self.to)
         return func
+    
+class acceleration(pattern):
+    def __init__(self,name,reduction,date):
+        pattern.__init__(self,name,reduction,date)
+        self.to=date
+
+    def g(self,t):
+        func= (t-self.to)**2
+        return func
 
 class sint(pattern):
     def __init__(self,name,reduction,date):
@@ -356,6 +370,10 @@ def date2dec(dates):
     inter = 'no'
 else:
     inter = arguments["--linear"]
+if arguments["--acceleration"] == None:
+    acc = 'no'
+else:
+    acc = arguments["--acceleration"]
 if arguments["--seasonal"] ==  None:
     seasonal = 'no'
 else:
@@ -654,6 +672,11 @@ def date2dec(dates):
       indexinter=index
       index = index + 1
 
+if acc =='yes':
+    basis.append(acceleration(name='acceleration', reduction='acc', date=datemin))
+    indexacc=index
+    index += 1 
+
 if degreeday=='yes':
    indexdd = index
    basis.append(stefan(name='degree-day',reduction='ddt',date=datemin,tcar1=ddt,tcar2=ddf))
@@ -961,8 +984,10 @@ def eq_cond(x):
     disp_seas = np.zeros((N))
     disp_seast = np.zeros((N))
 
-    if inter=='yes':
+    if inter=='yes' and acc=='no':
         lin[k] = np.dot(G[:,indexinter],m[indexinter])
+    if inter=='yes' and acc=='yes':
+        lin[k] = np.dot(G[:,indexinter],m[indexinter]) + np.dot(G[:, indexacc], m[indexacc])
 
     # plot data and model minus bperp error and seasonal terms
     if seasonal=='yes':
@@ -991,7 +1016,7 @@ def eq_cond(x):
             # convert between 0 and 2pi
             if phit<0: 
                 phit = phit + 2*np.pi
-
+            
             print(phi,phit)
             print(m[indexseas+1]/m[indexseas],m[indexseast+1]/m[indexseast] )
  
@@ -1038,7 +1063,14 @@ def eq_cond(x):
         )
   
     # plot data and model minus bperp error and linear term
+<<<<<<< HEAD
     if inter=='yes':
+=======
+    if inter=='yes' and acc=='no':
+        ax3.plot(x,disp-bperperr-lin,markers[jj],color=color,fillstyle=fillstyle,markersize=4,label='detrended data')
+        ax3.errorbar(x,disp-bperperr-lin,yerr = sigmad, ecolor=color,fmt='none', alpha=0.5)
+    elif inter=='yes' and acc=='yes':
+>>>>>>> 06ce944 (Add an acceleration fonction in the inversion [default: no])
         ax3.plot(x,disp-bperperr-lin,markers[jj],color=color,fillstyle=fillstyle,markersize=4,label='detrended data')
         ax3.errorbar(x,disp-bperperr-lin,yerr = sigmad, ecolor=color,fmt='none', alpha=0.5)
             
@@ -1061,8 +1093,13 @@ def eq_cond(x):
         G[:,Mbasis+l]=np.interp(tdec,tabx,kernels[l].g(k))
     model = np.dot(G,m)
     
-    if inter=='yes':
+    if inter=='yes' and acc=='no':
         model_lin = np.dot(G[:,indexinter],m[indexinter])
+<<<<<<< HEAD
+=======
+    elif inter=='yes' and acc=='yes':
+        model_lin = np.dot(G[:,indexinter],m[indexinter]) + np.dot(G[:, indexacc], m[indexacc])
+>>>>>>> 06ce944 (Add an acceleration fonction in the inversion [default: no])
    
     # we need to substrat bperp model to the whole model 
     if bperp=='yes':
@@ -1096,6 +1133,7 @@ def eq_cond(x):
 
     # plot model
     if inter=='yes':
+<<<<<<< HEAD
         if seasonal=='yes' and seasonalt=='no':
             ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexinter],amp,phi))
             ax3.plot(t,model-model_lin-model_bperp,'-r')
@@ -1105,6 +1143,26 @@ def eq_cond(x):
         elif seasonal=='no' and seasonalt=='yes':
             ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}'.format(m[indexinter],ampt,phit))
             ax3.plot(t,model-model_lin-model_bperp,'-r')
+=======
+        if seasonal=='yes' and seasonalt=='no' and acc=='no':
+            ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexinter],amp,phi))
+            ax3.plot(t,model-model_lin-model_bperp,'-r')
+        elif seasonal=='yes' and seasonalt=='no' and acc=='yes':
+            ax.plot(t,model-model_bperp,'-r',label='Acc: {:.2f}, Rate: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexacc], m[indexinter],amp,phi))
+            ax3.plot(t,model-model_lin-model_bperp,'-r')
+        elif seasonal=='yes' and seasonalt=='yes' and acc=='no':
+            ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexinter],ampt,phit,amp,phi))
+            ax3.plot(t,model-model_lin-model_bperp,'-r')
+        elif seasonal=='yes' and seasonalt=='yes' and acc=='yes':
+            ax.plot(t,model-model_bperp,'-r',label='Acc: {:.2f}, Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexacc], m[indexinter],ampt,phit,amp,phi))
+            ax3.plot(t,model-model_lin-model_bperp,'-r')        
+        elif seasonal=='no' and seasonalt=='yes' and acc=='no':
+            ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}'.format(m[indexinter],ampt,phit))
+            ax3.plot(t,model-model_lin-model_bperp,'-r')
+        elif seasonal=='no' and seasonalt=='yes' and acc=='yes':
+            ax.plot(t,model-model_bperp,'-r',label='Acc: {:.2f}, Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}'.format(m[indexacc], m[indexinter],ampt,phit))
+            ax3.plot(t,model-model_lin-model_bperp,'-r')
+>>>>>>> 06ce944 (Add an acceleration fonction in the inversion [default: no])
         else:
             ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}'.format(m[indexinter]))
             ax3.plot(t,model-model_lin-model_bperp,'-r')
@@ -1152,3 +1210,4 @@ def eq_cond(x):
 if plot == 'yes':
     plt.show()
 sys.exit()
+