Wrap for Geopsy-GPDC

gpdcwrap.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# -------------------------------------------------------------------
# Filename: gpdcwrap.py
#   Author: Damien Pageot
# ------------------------------------------------------------------
"""
Functions to use the Geopsy-gpdc engine.
"""

import numpy as np
from ctypes import CDLL, c_int, c_float, byref, POINTER, c_double
from numpy.ctypeslib import ndpointer, load_library

QGPCOREWAVE_PATH = "/home/geopsy/Codes/geopsy-install/lib/libQGpCoreWave.so"

import matplotlib.pyplot as plt

libCoreWave = load_library('libQGpCoreWave', QGPCOREWAVE_PATH)

def dispersion_curve_init(verbose):
    """
    Initialize the dispersion curve calculation.

    :param verbose: integer, 0 minimal ouput, 1 verbose output
    """
    libCoreWave.dispersion_curve_init_.argtypes = [POINTER(c_int)]
    libCoreWave.dispersion_curve_init_(byref(c_int(verbose)))
    return

def dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group):
    """
    Calculate the Rayleigh theoretical dispersion curve.

    :param nLayers: integer, number of layers
    :param h: double, thickness of layers (m)
    :param vp: double, P-wave velocity in each layer (m/s)
    :param vs: double, S-wave velocity in each layer (m/s)
    :param rho: double, density in each layer (kg/m3)
    :param nSamples: integer, number of frequency samples
    :param omega: double, angular frequencies (rad/s)
    :param nModes: integer, number of modes including fundamental
    :param slowness: double, output of slowness values
    :param group: integer, 0 for phase, 1 for group
    """
    libCoreWave.dispersion_curve_rayleigh_.argtypes = [ POINTER(c_int),
                                                        ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                        ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                        ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                        ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                        POINTER(c_int),
                                                        ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                        POINTER(c_int),
                                                        ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                        POINTER(c_int)]

    libCoreWave.dispersion_curve_rayleigh_(byref(c_int(nLayers)),
                                           h,
                                           vp,
                                           vs,
                                           rho,
                                           byref(c_int(nSamples)),
                                           omega,
                                           byref(c_int(nModes)),
                                           slowness,
                                           byref(c_int(group)))
    return

def dispersion_curve_love(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group):
    """
    Calculate the Love theoretical dispersion curve.

    :param nLayers: integer, number of layers
    :param h: double, thickness of layers (m)
    :param vp: double, P-wave velocity in each layer (m/s)
    :param vs: double, S-wave velocity in each layer (m/s)
    :param rho: double, density in each layer (kg/m3)
    :param nSamples: integer, number of frequency samples
    :param omega: double, angular frequencies (rad/s)
    :param nModes: integer, number of modes including fundamental
    :param slowness: double, output of slowness values
    :param group: integer, 0 for phase, 1 for group
    """
    libCoreWave.dispersion_curve_love_.argtypes = [ POINTER(c_int),
                                                    ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                    ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                    ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                    ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                    POINTER(c_int),
                                                    ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                    POINTER(c_int),
                                                    ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
                                                    POINTER(c_int)]

    libCoreWave.dispersion_curve_love_(byref(c_int(nLayers)),
                                       h,
                                       vp,
                                       vs,
                                       rho,
                                       byref(c_int(nSamples)),
                                       omega,
                                       byref(c_int(nModes)),
                                       slowness,
                                       byref(c_int(group)))
    return
    
   
if __name__ == "__main__":
    
    # Import numpy
    import numpy as np
    import matplotlib.pyplot as plt
    
    # Initialize dispersion curve calculation
    dispersion_curve_init(0)
    
    # Poisson ratio
    nu = 0.25
    
    # Number of modes (fundamental only = 1)
    nModes = 1
    
    # Phase velocity (group velocity =1)
    group = 0
    
    # Number of layers
    nLayers = 4
    
    # Layer thickness array
    h = np.array([2., 4., 8., 0.], dtype=np.float64)
    
    # P and S velocity arrays
    # With Vp = func(Vs, nu)
    vs = np.array([250., 500., 750., 1500.], dtype=np.float64)
    vp = vs * np.sqrt(2.*(1.-nu)/(1.-2.*nu))
    
    # Density array
    rho = np.array([1900., 1900., 1900., 1900.], dtype=np.float64)
    
    # Frequency and angular frequency arrays
    freq = np.arange(0.5, 80.5, 0.5, dtype=np.float64)
    omega = 2.*np.pi*freq 
    nSamples = len(freq)
    
    # Initialize the slowness array (output) 
    slowness = np.zeros(len(freq), dtype=np.float64)

    # Calculate the dispersion curve    
    dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group)
    
    
    # Plot
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Phase velocity (m/s)")
    plt.plot(freq, 1./slowness)
    plt.show()