# Source code for geophpy.processing.magnetism

    # -*- coding: utf-8 -*-
    '''
        geophpy.processing.magnetism
        ----------------------------

        DataSet Object general magnetism processing routines.

        :copyright: Copyright 2014 Lionel Darras, Philippe Marty and contributors, see AUTHORS.
        :license: GNU GPL v3.

    '''
    import numpy as np
    from scipy.linalg import lu_solve, lu_factor
    from geophpy.operation.general import *

    def logtransform(dataset, multfactor=5):
       '''
       To transform the data in logarihtmics values

       Parameters :

       :dataset: DataSet Object to do the treatment

       :multfactor: multiplying factor

       '''

       val = dataset.data.z_image

       # Calculation of lines and columns numbers
       ny = len(val)
       nx = len(val[0])

       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if (val[iy][ix] != np.nan) :                          # if valid value
                val[iy][ix] = val[iy][ix] * multfactor             # multiplying

       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((val[iy][ix] != np.nan) and (val[iy][ix] < -1)):
                val[iy][ix] = -np.log10(-val[iy][ix])

       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((val[iy][ix] != np.nan) and (val[iy][ix] > 1)):
                val[iy][ix] = np.log10(val[iy][ix])

       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((val[iy][ix] != np.nan) and (val[iy][ix] > -1) and (val[iy][ix] < 1)) :
                val[iy][ix] = np.nan

       dataset.data.z_image = val

    def polereduction(dataset, apod=0, inclineangle = 65, alphaangle = 0):
       '''
       To reduce at the magnetic pôle from a DataSet Object.

       Parameters:

       :dataset: DataSet Object to do the treatment

       :apod: apodisation factor (%), to reduce bord effects

       :inclinangle: incline angle of the magnetic field

       :alphaangle: angle of magnetic field and y axe

       '''

       val = dataset.data.z_image

       # Calculation of lines and columns numbers
       ny = len(val)
       nx = len(val[0])
       # Calculation of total values number
       nt = nx * ny

       nan_indexes = np.isnan(val)                                 # searches the 'nan' indexes in the initial array
       if(np.isnan(val).any()):                                    # if at least 1 'nan' is detected in the array,

          x = np.linspace(dataset.info.x_min, dataset.info.x_max, nx, endpoint=True)
          y = np.linspace(dataset.info.y_min, dataset.info.y_max, ny, endpoint=True)
          _fillnanvalues(x, y, val.T)                              # fills the 'nan' values by interpolated values

       if (apod > 0):
          apodisation2d(val, apod)

       # Meaning calculation
       mean = val.mean()

       # complex conversion, val[][] -> cval[][]
       cval = np.array(val, dtype=complex)

       # fast fourier series calculation
       cvalfft = np.fft.fft2(cval)

       # filter application
       # deg->rad angle conversions
       rinc = (inclineangle*np.pi)/180 # inclination in radians
       ralpha = (alphaangle*np.pi)/180  # alpha angle in radians

       cosu = np.absolute(np.cos(rinc) * np.cos(ralpha))
       cosv = np.absolute(np.cos(rinc) * np.sin(ralpha))
       cosz = np.sin(rinc)

       deltax = dataset.info.x_gridding_delta
       deltay = dataset.info.y_gridding_delta

       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                cz = np.sqrt(np.square(cu) + np.square(cv))
                cred = complex(cz*cosz, cu*cosu + cv*cosv)
                cvalfft[iy][ix] = ((cvalfft[iy][ix] * np.square(cz))/np.square(cred))

       # FFT inversion
       icvalfft = np.fft.ifft2(cvalfft)
       val = icvalfft.real + mean                                  # real part added to mean calculation

       val[nan_indexes] = np.nan                                   # set the initial 'nan' values, as in the initial array

       dataset.data.z_image = val                                  # to get the 'z_image' array in (x,y) format.

    def continuation(dataset, prosptech, apod, downsensoraltitude, upsensoraltitude, continuationflag, continuationvalue):
       '''
       Continuation of the magnetic field

       Parameters :

       :dataset: DataSet Object to do the treatment

       :prosptech: Prospection technical

       :apod: apodisation factor (%), to reduce bord effects

       :downsensoraltitude: Altitude of the down magnetic sensor

       :upsensoraltitude: Altitude of the up magnetic sensor

       :continuationflag: Continuation flag, False if not continuation

       :continuationvalue: Continuation altitude, greater than 0 if upper earth ground, lower than 0 else

       '''

       val = dataset.data.z_image

       # Calculation of lines and columns numbers
       ny = len(val)
       nx = len(val[0])
       # Calculation of total values number
       nt = nx * ny

       nan_indexes = np.isnan(val)                                 # searches the 'nan' indexes in the initial array
       if(np.isnan(val).any()):                                    # if at least 1 'nan' is detected in the array,

          x = np.linspace(dataset.info.x_min, dataset.info.x_max, nx, endpoint=True)
          y = np.linspace(dataset.info.y_min, dataset.info.y_max, ny, endpoint=True)
          _fillnanvalues(x, y, val.T)                                # fills the 'nan' values by interpolated values

       if (apod > 0):
          apodisation2d(val, apod)

       dz = downsensoraltitude - upsensoraltitude
       zp = - continuationvalue + downsensoraltitude

       # complex conversion, val[][] -> cval[][]
       cval = np.array(val, dtype=complex)

       # fast fourier series calculation
       cvalfft = np.fft.fft2(cval)

       deltax = dataset.info.x_gridding_delta
       deltay = dataset.info.y_gridding_delta

       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                cz = np.sqrt(np.square(cu) + np.square(cv)) * (2. * np.pi)
                if (continuationflag == True):
                   ce = np.exp(zp * cz)
                   cvalfft[iy][ix] = (cvalfft[iy][ix] * complex(ce, 0.))
                if (prosptech == prosptechlist[1]):                # if prospection technic is magnetic field gradient
                   cdz = 1. - np.exp(dz * cz)
                   cvalfft[iy][ix] = (cvalfft[iy][ix] / complex(cdz, 0.))

       # FFT inversion
       icvalfft = np.fft.ifft2(cvalfft)
       val = icvalfft.real

       val[nan_indexes] = np.nan                                   # set the initial 'nan' values, as in the initial array

       dataset.data.z_image = val                                  # to get the 'z_image' array in (x,y) format.

    def eulerdeconvolution(dataset, xmin, xmax, ymin, ymax, apod=0, nflag=False, nvalue=0):
       '''
       cf. dataset
       '''

       val = dataset.data.z_image

       # Calculation of lines and columns numbers
       ny = len(val)
       nx = len(val[0])
       # Calculation of total values number
       nt = nx * ny

       x = np.linspace(dataset.info.x_min, dataset.info.x_max, nx, endpoint=True)
       y = np.linspace(dataset.info.y_min, dataset.info.y_max, ny, endpoint=True)

       nan_indexes = np.isnan(val)                                 # searches the 'nan' indexes in the initial array
       if(np.isnan(val).any()):                                    # if at least 1 'nan' is detected in the array,
          _fillnanvalues(x, y, val.T)                                # fills the 'nan' values by interpolated values

       if (apod > 0):
          apodisation2d(val, apod)

       # complex conversion, val[][] -> cval[][]
       cval = np.array(val, dtype=complex)
       val = cval.real

       deltax = dataset.info.x_gridding_delta
       deltay = dataset.info.y_gridding_delta

       # fast fourier series calculation
       cvalfft = np.fft.fft2(cval)
       cvalfftdx = np.array(np.zeros((ny,nx)), dtype=complex)
       cvalfftdy = np.array(np.zeros((ny,nx)), dtype=complex)
       cvalfftdz = np.array(np.zeros((ny,nx)), dtype=complex)

       # derivations calculation
       for iy in range(ny):                                        # for each column
          for ix in range(nx):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                cz = np.sqrt(np.square(cu) + np.square(cv))
                # x derivation
                cvalfftdx[iy][ix] = cvalfft[iy][ix]*complex(0., cu*np.pi*2.)
                # y derivation
                cvalfftdy[iy][ix] = cvalfft[iy][ix]*complex(0., cv*np.pi*2.)
                # z derivation
                cvalfftdz[iy][ix] = cvalfft[iy][ix]*cz*np.pi*2.
             else:
                cvalfftdx[0][0] = complex(0., 0.)
                cvalfftdy[0][0] = complex(0., 0.)
                cvalfftdz[0][0] = complex(0., 0.)

       # FFT inversion
       valdx = (np.fft.ifft2(cvalfftdx)).real
       valdy = (np.fft.ifft2(cvalfftdy)).real
       valdz = (np.fft.ifft2(cvalfftdz)).real

       ixmin, xnearestmin = _searchnearest(x, xmin)
       ixmax, xnearestmax = _searchnearest(x, xmax)

       iymin, ynearestmin = _searchnearest(y, ymin)
       iymax, ynearestmax = _searchnearest(y, ymax)

       if ((nflag == 0) or (nvalue==0)):   # automatic calculation structural index
                                           # 4 unknowns equation to solve
          b = np.zeros(4)
          A = np.zeros((4,4))
          for l in range(iymin, iymax+1):
             for c in range(ixmin, ixmax+1):
                coef = x[c] * valdx[l][c] + y[l] * valdy[l][c]
                b[0] += coef * val[l][c]
                b[1] += coef * valdx[l][c]
                b[2] += coef * valdy[l][c]
                b[3] += coef * valdz[l][c]

                A[0][0] += np.square(val[l][c])
                A[0][1] += val[l][c] * valdx[l][c]
                A[0][2] += val[l][c] * valdy[l][c]
                A[0][3] += val[l][c] * valdz[l][c]
                A[1][1] += np.square(valdx[l][c])
                A[1][2] += valdx[l][c] * valdy[l][c]
                A[1][3] += valdx[l][c] * valdz[l][c]
                A[2][2] += np.square(valdy[l][c])
                A[2][3] += valdy[l][c] * valdz[l][c]
                A[3][3] += np.square(valdz[l][c])

          # symmetry
          A[1][0] = A[0][1]
          A[2][0] = A[0][2]
          A[3][0] = A[0][3]
          A[2][1] = A[1][2]
          A[3][1] = A[1][3]
          A[3][2] = A[2][3]

          x = _gaussresolution(A, b)
          nvalue = -x[0]
          xpt = x[1]
          ypt = x[2]
          depth = x[3]

       else :                              # fixed structural index
                                           # 3 unknowns equation to solve
          b = np.zeros(3)
          A = np.zeros((3,3))
          for l in range(iymin, iymax+1):
             for c in range(ixmin, ixmax+1):
                coef = nvalue * val[l][c] + x[c] * valdx[l][c] + y[l] * valdy[l][c]
                b[0] += coef * valdx[l][c]
                b[1] += coef * valdy[l][c]
                b[2] += coef * valdz[l][c]

                A[0][0] += np.square(valdx[l][c])
                A[0][1] += valdx[l][c] * valdy[l][c]
                A[0][2] += valdx[l][c] * valdz[l][c]
                A[1][1] += np.square(valdy[l][c])
                A[1][2] += valdy[l][c] * valdz[l][c]
                A[2][2] += np.square(valdz[l][c])

          # symmetry
          A[1][0] = A[0][1]
          A[2][0] = A[0][2]
          A[2][1] = A[1][2]

          x = _gaussresolution(A, b)
          xpt = x[0]
          ypt = x[1]
          depth = x[2]

       return nvalue, xpt, ypt, depth, xnearestmin, ynearestmin, xnearestmax, ynearestmax

    def _searchnearest(array, val):
       '''
       Search the nearset value of val, in the array

       Parameters :

       :array:  1D array ordered to find the nearest value of val

       :val: value to find

       Returns :

       :index: index in the array of the nearest value found

       :valnearest: value of the nearest value found

       '''
       n = len(array)
       for i in range(n):
          gap = abs(val-array[i])
          if ((i==0) or (gap<mingap)):
             mingap = gap
             index = i
       return index, array[index]

    def _gaussresolution(A,b):
       '''
       Solving Ax=b equation with gauss pivot and LU factorisation
       '''

       n=len(A)
       # gauss pivot calculation
       for i1 in range(n-1):
          # find the partial pivot
          l=i1
          xam = np.absolute(A[i1][i1])
          for j in range(i1+1,n):
             x = np.absolute(A[j][i1])
             if (x>xam):
                xam = x
                l=j

          # set pivot at its place, swapping pivot and current lines
          if (l>i1):
             aux = b[l]
             b[l] = b[i1]
             b[i1] = aux
             for j in range(n):
                aux = A[l][j]
                A[l][j] = A[i1][j]
                A[i1][j] = aux

       # solves Ax=b using LU factorisation
       x = lu_solve(lu_factor(A), b)

       return x

    def analyticsignal(dataset, apod=0):
       '''
       cf. dataset
       '''

       val = dataset.data.z_image

       # Calculation of lines and columns numbers
       ny = len(val)
       nx = len(val[0])
       # Calculation of total values number
       nt = nx * ny

       nan_indexes = np.isnan(val)                                 # searches the 'nan' indexes in the initial array
       if(np.isnan(val).any()):                                    # if at least 1 'nan' is detected in the array,

          x = np.linspace(dataset.info.x_min, dataset.info.x_max, nx, endpoint=True)
          y = np.linspace(dataset.info.y_min, dataset.info.y_max, ny, endpoint=True)
          _fillnanvalues(x, y, val.T)                              # fills the 'nan' values by interpolated values

       if (apod > 0):
          apodisation2d(val, apod)

       # complex conversion, val[][] -> cval[][]
       cval = np.array(val, dtype=complex)

       deltax = dataset.info.x_gridding_delta
       deltay = dataset.info.y_gridding_delta

       # fast fourier series calculation
       cvalfft = np.fft.fft2(cval)
       cvalfftdx = np.fft.fft2(cval)
       cvalfftdy = np.fft.fft2(cval)
       cvalfftdz = np.fft.fft2(cval)

       # dx derivations calculation
       cvalfftdx[0][0] = complex(0., 0.)
       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                # x derivation
                cvalfftdx[iy][ix] = cvalfft[iy][ix]*complex(0., cu*np.pi*2)

       # FFT inversion
       icvalfft = np.fft.ifft2(cvalfftdx)
       val = np.square(icvalfft.real)

       # dy derivations calculation
       cvalfftdy[0][0] = complex(0., 0.)
       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                # y derivation
                cvalfftdy[iy][ix] = cvalfft[iy][ix]*complex(0., cv*np.pi*2)

       # FFT inversion
       icvalfft = np.fft.ifft2(cvalfftdy)
       val = val + np.square(icvalfft.real)

       # dz derivations calculation
       cvalfftdz[0][0] = complex(0., 0.)
       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                cz = np.sqrt(np.square(cu) + np.square(cv))
                # z derivation
                cvalfftdz[iy][ix] = cvalfft[iy][ix]*cz*np.pi*2

       # FFT inversion
       icvalfft = np.fft.ifft2(cvalfftdz)
       val = np.sqrt(val + np.square(icvalfft.real))

       val[nan_indexes] = np.nan                                   # set the initial 'nan' values, as in the initial array

       dataset.data.z_image = val                                  # to get the 'z_image' array in (x,y) format.

    def gradmagfieldconversion(dataset, prosptechused, prosptechsim, apod=0, downsensoraltused = 0.3, upsensoraltused = 1.0, downsensoraltsim = 0.3, upsensoraltsim = 1.0, inclineangle = 65, alphaangle = 0):
       '''
       cf. dataset
       '''

       case, sense = _definecaseandsense(prosptechused, prosptechsim)

       if ((case != None) and (sense != None)):

          val = dataset.data.z_image

          # Calculation of lines and columns numbers
          ny = len(val)
          nx = len(val[0])
          # Calculation of total values number
          nt = nx * ny

          nan_indexes = np.isnan(val)                                 # searches the 'nan' indexes in the initial array
          if(np.isnan(val).any()):                                    # if at least 1 'nan' is detected in the array,

             x = np.linspace(dataset.info.x_min, dataset.info.x_max, nx, endpoint=True)
             y = np.linspace(dataset.info.y_min, dataset.info.y_max, ny, endpoint=True)
             _fillnanvalues(x, y, val.T)                              # fills the 'nan' values by interpolated values

          if (apod > 0):
             apodisation2d(val, apod)

          # complex conversion, val[][] -> cval[][]
          cval = np.array(val, dtype=complex)

          deltax = dataset.info.x_gridding_delta
          deltay = dataset.info.y_gridding_delta

          # fast fourier series calculation
          cvalfft = np.fft.fft2(cval)

          if (case >= 2):
             # deg->rad angle conversions
             rinc = (inclineangle*np.pi)/180 # inclination in radians
             ralpha = (alphaangle*np.pi)/180  # alpha angle in radians
             cosu = np.absolute(np.cos(rinc) * np.cos(ralpha))
             cosv = np.absolute(np.cos(rinc) * np.sin(ralpha))
             cosz = np.sin(rinc)

          for ix in range(nx):                                        # for each column
             for iy in range(ny):                                     # for each line
                if ((ix != 0) or (iy != 0)):                          # if not first point of array
                   cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                   cz = np.sqrt(np.square(cu) + np.square(cv))

                   if (case == 1):
                      c1 = (downsensoraltsim - downsensoraltused)*2.*np.pi*cz
                      c2 = (upsensoraltsim - downsensoraltsim)*2.*np.pi*cz
                      c = np.exp(-c1) - np.exp(-c2)
                   elif (case == 2):
                      c1 = (downsensoraltsim - downsensoraltused)*2.*np.pi*cz
                      c2 = (upsensoraltsim - downsensoraltsim)*2.*np.pi*cz
                      cred = complex(cz*cosz, cu*cosu + cv*cosv)
                      c = ((np.exp(-c1) - np.exp(-c2))*cz)/cred
                   else :                  # case = 3
                      cred = complex(cz*cosz, cu*cosu + cv*cosv)
                      c = cz/cred

                   if (sense == 0):
                      coef = c
                   else:
                      coef = 1./c

                   cvalfft[iy][ix] = cvalfft[iy][ix] * coef

          # FFT inversion
          icvalfft = np.fft.ifft2(cvalfft)
          val = icvalfft.real

          val[nan_indexes] = np.nan                                   # set the initial 'nan' values, as in the initial array

          dataset.data.z_image = val                                  # to get the 'z_image' array in (x,y) format.

    def _definecaseandsense( prosptechused, prosptechsim):
       '''
       define case and sense of conversion
       case 1 = "Magnetic field" <-> "Magnetic field gradient"
       case 2 = "Magnetic field" <-> "Vertical component gradient"
       case 3 = "Magnetic field gradient" <-> "Vertical component gradient"
       '''

       case = None
       sense = None
       if ((prosptechused == prosptechlist[0]) and (prosptechsim == prosptechlist[1])) :
          case = 1
          sense = 0
       elif ((prosptechused == prosptechlist[0]) and (prosptechsim == prosptechlist[2])) :
          case = 2
          sense = 0
       elif ((prosptechused == prosptechlist[1]) and (prosptechsim == prosptechlist[0])) :
          case = 1
          sense = 1
       elif ((prosptechused == prosptechlist[1]) and (prosptechsim == prosptechlist[2])) :
          case = 3
          sense = 0
       elif ((prosptechused == prosptechlist[2]) and (prosptechsim == prosptechlist[0])) :
          case = 2
          sense = 1
       elif ((prosptechused == prosptechlist[2]) and (prosptechsim == prosptechlist[1])) :
          case = 3
          sense = 1

       return case, sense

    def susceptibility(dataset, prosptech, apod=0, downsensoraltitude = 0.3, upsensoraltitude = 1.0, calculationdepth=.0, stratumthickness=1.0, inclineangle = 65, alphaangle = 0):
       '''
       cf. dataset
       '''
       val = dataset.data.z_image

       # Calculation of lines and columns numbers
       ny = len(val)
       nx = len(val[0])
       # Calculation of total values number
       nt = nx * ny

       nan_indexes = np.isnan(val)                                 # searches the 'nan' indexes in the initial array
       if(np.isnan(val).any()):                                    # if at least 1 'nan' is detected in the array,

          x = np.linspace(dataset.info.x_min, dataset.info.x_max, nx, endpoint=True)
          y = np.linspace(dataset.info.y_min, dataset.info.y_max, ny, endpoint=True)
          _fillnanvalues(x, y, val.T)                                # fills the 'nan' values by interpolated values

       if (apod > 0):
          apodisation2d(val, apod)

       # complex conversion, val[][] -> cval[][]
       coef = 1./(400.*np.pi)
       cval = np.array(val*coef, dtype=complex)

       # fast fourier series calculation
       cvalfft = np.fft.fft2(cval)

       # filter application
       # deg->rad angle conversions
       rinc = (inclineangle*np.pi)/180 # inclination in radians
       ralpha = (alphaangle*np.pi)/180  # alpha angle in radians

       cosu = np.absolute(np.cos(rinc) * np.cos(ralpha))
       cosv = np.absolute(np.cos(rinc) * np.sin(ralpha))
       cosz = np.sin(rinc)

       dz = downsensoraltitude - upsensoraltitude
       zp = calculationdepth + downsensoraltitude

       deltax = dataset.info.x_gridding_delta
       deltay = dataset.info.y_gridding_delta

       if (inclineangle == 0):
          teta = np.pi/2
       elif (inclineangle == 90) :    # pi/2
          teta = np.pi
       else :
          teta = np.pi + np.arctan(-2/np.tan(rinc))

       # Field Module, medium field along I inclination calcultation if bipolar field
       FM = 14.722 * (4*(np.square(np.cos(teta)) + np.square(np.sin(teta))))

       cvalfft[0][0] = 0.
       for ix in range(nx):                                        # for each column
          for iy in range(ny):                                     # for each line
             if ((ix != 0) or (iy != 0)):                          # if not first point of array
                cu, cv = _freq(ix, iy, deltax, deltay, nx, ny)     # complex u and v calculation
                # continuation
                cz = np.sqrt(np.square(cu) + np.square(cv))
                cvalfft[iy][ix] = cvalfft[iy][ix] * complex(np.exp(2*np.pi*cz*zp), 0.)
                # pole reduction with potential calculation(and not with field as standard pole reduction)
                cred = complex(cz*cosz, cu*cosu + cv*cosv)
                cvalfft[iy][ix] = (cvalfft[iy][ix] * complex(cz,0))/(2*np.pi*np.square(cred))

                if (prosptech == prosptechlist[1]):                         # if prospection technic is magnetic field gradient
                   cvalfft[iy][ix] = cvalfft[iy][ix] / (1-np.exp(dz*cz))

       # FFT inversion
       icvalfft = np.fft.ifft2(cvalfft)

       # Equivalent stratum thickness
       val = (icvalfft.real*2*100000)/(FM*stratumthickness)         # 100000 because susceptibiliy in 10^-5 SI

       val[nan_indexes] = np.nan                                   # set the initial 'nan' values, as in the initial array

       dataset.data.z_image = val                                  # to get the 'z_image' array in (x,y) format.

    prosptechlist = ["Magnetic field", "Magnetic field gradient", "Vertical component gradient"]

    def getprosptechlist():
       '''
       Get list of prospection technicals availables

       Returns : list of prospection technicals

       '''

       return prosptechlist

    def _freq(ix, iy, deltax, deltay, nc, nl):
       '''
       Calculation of spatials frequencies u and v
       '''
       nyd = 1 + nl/2
       nxd = 1 + nc/2

       if (iy < nyd):
          cv = (float(iy))/((nl-1)*deltay)
       else:
          cv = float(iy-nl)/((nl-1)*deltay)

       if (ix < nxd):
          cu = (float(ix))/((nc-1)*deltax)
       else:
          cu = float(ix-nc)/((nc-1)*deltax)

       return cu, cv

    def _fillnanvalues(xi, yi, val):
       '''
       Fills the 'nan' values by interpolated values using simple spline interpolation method !
       '''
       for profile in val:                                         # for each profile,
          if (np.isnan(profile).any()):                            # if one 'nan' at least in the profile, does the completion
             nan_indexes = np.isnan(profile)
             data_indexes = np.logical_not(nan_indexes)
             valid_data = profile[data_indexes]
             interpolated = np.interp(nan_indexes.nonzero()[0], data_indexes.nonzero()[0], valid_data)
             profile[nan_indexes] = interpolated