Adding C-kernel

.gitignore

@@ -26,3 +26,7 @@ empymod/version.py
 build/
 dist/
 empymod.egg-info/
+
+# C
+*.o
+*.so

Makefile

@@ -14,9 +14,11 @@ help:
 	@echo ""
 
 install:
+	cd empymod && make && cd ..
 	python -m pip install -e .
 
 dev-install:
+	cd empymod && make && cd ..
 	python -m pip install -e .[all]
 
 .ONESHELL:
@@ -46,6 +48,7 @@ linkcheck:
 	cd docs && make linkcheck
 
 clean:
+	cd empymod && make deepclean && cd ..
 	python -m pip uninstall empymod -y
 	rm -rf build/ dist/ .eggs/ empymod.egg-info/ empymod/version.py  # build
 	rm -rf */__pycache__/ */*/__pycache__/      # python cache

empymod/Makefile

@@ -0,0 +1,32 @@
+#CC= /opt/local/bin/gcc-mp-13
+
+all: fields.so reflections.so greenfct.so wavenumber.so
+
+fields.so: fields.o
+	gcc -shared -o $@ fields.o -lm
+
+reflections.so: reflections.o
+	gcc -shared -o $@ reflections.o -lm
+
+greenfct.so: greenfct.o reflections.o fields.o
+	gcc -shared -o $@ $^ -lm
+
+wavenumber.so: wavenumber.o greenfct.o reflections.o fields.o
+	gcc -shared -o $@ $^ -lm
+
+#CFLAGS= -fPIC -Ofast -ftree-loop-vectorize -fopt-info-optimized  -fopt-info-vec-missed  -march=znver2
+#-fopt-info-optimized -fopenmp
+CFLAGS= -fPIC -Ofast -ftree-loop-vectorize -ffast-math -fopt-info-optimized -funroll-all-loops 
+CFLAGS= -fPIC -O3 -ftree-loop-vectorize -funroll-all-loops 
+#CFLAGS= -fPIC -O0 -g
+
+.SUFFIXES : .o .c 
+.c.o    :
+	$(CC) -c $(CFLAGS) $(OPTC) $<
+
+clean:
+	rm -f *.o
+
+deepclean:
+	rm -rf *.o *.so __pycache__
+

empymod/_ckernel.py

@@ -0,0 +1,133 @@
+import os
+import ctypes as ct
+
+import numpy as np
+
+
+C_DOUBLEP = ct.POINTER(ct.c_double)
+path = os.path.dirname(os.path.realpath(__file__))
+cwavenumber = np.ctypeslib.load_library("wavenumber", path)
+cgreenfct = np.ctypeslib.load_library("greenfct", path)
+creflections = np.ctypeslib.load_library("reflections", path)
+cfields = np.ctypeslib.load_library("fields", path)
+
+
+def wavenumber(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd,
+               ab, xdirect, msrc, mrec):
+    """C-Wrapper for wavenumber."""
+    nfreq, nlayer = etaH.shape
+    noff, nlambda = lambd.shape
+    PJ0 = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    PJ1 = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    PJ0b = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    cwavenumber.wavenumber(
+        int(nfreq),
+        int(noff),
+        int(nlayer),
+        int(nlambda),
+        ct.c_double(zsrc),
+        ct.c_double(zrec),
+        int(lsrc),
+        int(lrec),
+        depth.ctypes.data_as(C_DOUBLEP),
+        etaH.ctypes.data_as(C_DOUBLEP),
+        etaV.ctypes.data_as(C_DOUBLEP),
+        zetaH.ctypes.data_as(C_DOUBLEP),
+        zetaV.ctypes.data_as(C_DOUBLEP),
+        lambd.ctypes.data_as(C_DOUBLEP),
+        int(ab),
+        int(xdirect),
+        int(msrc),
+        int(mrec),
+        PJ0.ctypes.data_as(C_DOUBLEP),
+        PJ1.ctypes.data_as(C_DOUBLEP),
+        PJ0b.ctypes.data_as(C_DOUBLEP),
+    )
+    if ab not in [11, 22, 24, 15, 33]:
+        PJ0 = None
+    if ab not in [11, 12, 21, 22, 14, 24, 15, 25]:
+        PJ0b = None
+    if ab in [33, ]:
+        PJ1 = None
+    return PJ0, PJ1, PJ0b
+
+
+def greenfct(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd,
+             ab, xdirect, msrc, mrec):
+    """C-Wrapper for greenfct."""
+    nfreq, nlayer = etaH.shape
+    noff, nlambda = lambd.shape
+    GTM = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    GTE = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    cgreenfct.greenfct(
+        int(nfreq),
+        int(noff),
+        int(nlayer),
+        int(nlambda),
+        ct.c_double(zsrc),
+        ct.c_double(zrec),
+        int(lsrc),
+        int(lrec),
+        depth.ctypes.data_as(C_DOUBLEP),
+        etaH.ctypes.data_as(C_DOUBLEP),
+        etaV.ctypes.data_as(C_DOUBLEP),
+        zetaH.ctypes.data_as(C_DOUBLEP),
+        zetaV.ctypes.data_as(C_DOUBLEP),
+        lambd.ctypes.data_as(C_DOUBLEP),
+        int(ab),
+        int(xdirect),
+        int(msrc),
+        int(mrec),
+        GTM.ctypes.data_as(C_DOUBLEP),
+        GTE.ctypes.data_as(C_DOUBLEP),
+    )
+    return GTM, GTE
+
+
+def reflections(depth, e_zH, Gam, lrec, lsrc):
+    """C-Wrapper for reflections."""
+    nfreq, noff, nlayer, nlambda = Gam.shape
+    maxl = max([lrec, lsrc])
+    minl = min([lrec, lsrc])
+    nl = maxl-minl+1
+    Rp = np.zeros((nfreq, noff, nl, nlambda), dtype=complex)
+    Rm = np.zeros((nfreq, noff, nl, nlambda), dtype=complex)
+    creflections.reflections(
+        int(nfreq),
+        int(noff),
+        int(nlayer),
+        int(nlambda),
+        depth.ctypes.data_as(C_DOUBLEP),
+        e_zH.ctypes.data_as(C_DOUBLEP),
+        Gam.ctypes.data_as(C_DOUBLEP),
+        int(lrec),
+        int(lsrc),
+        Rp.ctypes.data_as(C_DOUBLEP),
+        Rm.ctypes.data_as(C_DOUBLEP),
+    )
+    return Rp, Rm
+
+
+def fields(depth, Rp, Rm, Gam, lrec, lsrc, zsrc, ab, TM):
+    """C-Wrapper for fields."""
+    nfreq, noff, nlayer, nlambda = Rp.shape
+    Pu = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    Pd = np.zeros((nfreq, noff, nlambda), dtype=complex)
+    cfields.fields(
+        int(nfreq),
+        int(noff),
+        int(nlayer),
+        int(nlambda),
+        depth.ctypes.data_as(C_DOUBLEP),
+        Rp.ctypes.data_as(C_DOUBLEP),
+        Rm.ctypes.data_as(C_DOUBLEP),
+        Gam.ctypes.data_as(C_DOUBLEP),
+        int(lrec),
+        int(lsrc),
+        ct.c_double(zsrc),
+        int(ab),
+        int(TM),
+        Pu.ctypes.data_as(C_DOUBLEP),
+        Pd.ctypes.data_as(C_DOUBLEP),
+    )
+    return Pu, Pd