diff --git a/fdelmodc/Makefile b/fdelmodc/Makefile
index 4214ad0836e2eb3fc6b8d426ffa5ec1138c22c57..f214df9f50540585f5d0f3aa2949d3731d1441b8 100644
--- a/fdelmodc/Makefile
+++ b/fdelmodc/Makefile
@@ -34,6 +34,7 @@ SRCC = $(PRG).c \
getWaveletInfo.c \
getModelInfo.c \
applySource.c \
+ applyRupture.c \
getRecTimes.c \
getBeamTimes.c \
writeSnapTimes.c \
diff --git a/fdelmodc/applyRupture.c b/fdelmodc/applyRupture.c
new file mode 100644
index 0000000000000000000000000000000000000000..12212df88beee083d2fdf11280b4a699faf05c8f
--- /dev/null
+++ b/fdelmodc/applyRupture.c
@@ -0,0 +1,96 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc.h"
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+float stresscalc(float Tu, float Ts, float D, float dD, float D0, float V0);
+void slipcalc(float *D, float *dD, float *vx, float dx, int izsrc, int ixsrc1, int ixsrc2, int n1);
+
+void applyRupture(modPar mod, int itime, float *vx, float *txz, int verbose)
+{
+
+ float *D, *dD, Tout;
+ int irup;
+ static float zrup, xrup1, xrup2, dT, T0, Tu, Ts, D0, V0, dx, trup;
+ static int izrup, ixrup1, ixrup2, itrup, nrup, n1, first=1;
+
+ if (first==1) {
+ if(!getparfloat("zrup", &zrup)) zrup = 0.0;
+ if(!getparfloat("xrup1", &xrup1)) xrup1 = 0.0;
+ if(!getparfloat("xrup2", &xrup2)) xrup2 = 0.0;
+ if(!getparfloat("trup", &trup)) trup = 0.0;
+ if(!getparfloat("V0", &V0)) V0 = 0.0;
+ if(!getparfloat("D0", &D0)) D0 = 0.0;
+ if(!getparfloat("Tu", &Tu)) Tu = 0.0;
+ if(!getparfloat("Ts", &Ts)) Ts = 0.0;
+ if(!getparfloat("dT", &dT)) dT = 0.0;
+
+ dx = mod.dx;
+ n1 = mod.naz;
+
+ ixrup1 = NINT(xrup1/dx);
+ ixrup2 = NINT(xrup2/dx);
+ itrup = NINT(trup/mod.dt);
+ nrup = ixrup2 - ixrup1;
+ first = 0;
+ }
+
+ if (itime>itrup) {
+ D = (float *)calloc(nrup,sizeof(float));
+ dD = (float *)calloc(nrup,sizeof(float));
+
+ slipcalc(D, dD, vx, dx, izrup, ixrup1, ixrup2, n1);
+
+ for (irup=ixrup1; irup<ixrup2; irup++){
+ vmess("D[%d]=%.3e",D[irup-ixrup1]);
+ if (fabs(dD[irup-ixrup1])>0.001) {
+ Tout = stresscalc(Tu, Ts, D[irup-ixrup1], dD[irup-ixrup1], D0, V0);
+ txz[irup*n1+izrup-1] += Tout - T0;
+ txz[irup*n1+izrup] += Tout - T0;
+ }
+ else {
+ vx[irup*n1+izrup] = 0.0;
+ }
+ }
+
+ free(D); free(dD);
+ }
+ else if (itime==itrup) {
+ for (irup=ixrup1; irup<ixrup2; irup++){
+ txz[irup*n1+izrup-1] += dT;
+ txz[irup*n1+izrup] += dT;
+ }
+ }
+
+}
+
+float stresscalc(float Tu, float Ts, float D, float dD, float D0, float V0)
+{
+ float T1, T2, Tout;
+
+ T1 = Tu*(1-D/D0);
+ T2 = Ts*(V0/(V0+dD));
+
+ if (fabs(T1) > fabs(T2)) Tout=T1;
+ else Tout = T2;
+
+ return Tout;
+}
+
+void slipcalc(float *D, float *dD, float *vx, float dx, int izsrc, int ixsrc1, int ixsrc2, int n1)
+{
+ int ix;
+
+ for (ix=ixsrc1; ix<ixsrc2; ix++) {
+ dD[ix-ixsrc1] = vx[ix*n1+izsrc+1] - vx[ix*n1+izsrc-1];
+ }
+
+ for (ix=ixsrc1+1; ix<ixsrc2; ix++) {
+ D[ix-ixsrc1] = (dx/2.0)*(dD[ix-ixsrc1-1]+dD[ix-ixsrc1]);
+ }
+ D[0] = (dx/2.0)*dD[0];
+
+}
\ No newline at end of file
diff --git a/fdelmodc/elastic4.c b/fdelmodc/elastic4.c
index bfaee6e6780cff467a1a910ffb5c8524e0845485..532f1cf967119fd86ea5662ec20446a0a49f5e1c 100644
--- a/fdelmodc/elastic4.c
+++ b/fdelmodc/elastic4.c
@@ -17,6 +17,8 @@ int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float
int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
+void applyRupture(modPar mod, int itime, float *vx, float *txz, int verbose);
+
int elastic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx,
float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
*l2m, float *lam, float *mul, int verbose)
@@ -103,7 +105,7 @@ float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
}
/* Add force source */
- if (src.type > 5) {
+ if (src.type > 5 && src.type!=20) {
applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
}
@@ -144,6 +146,11 @@ float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
if (src.type < 6) {
applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
}
+
+ /* Add rupture */
+ if (src.type==20) {
+ applyRupture(mod, itime, vx, txz, verbose);
+ }
/* check if there are sources placed on the boundaries */
storeSourceOnSurface(mod, src, bnd, ixsrc, izsrc, vx, vz, tzz, txx, txz, verbose);
diff --git a/marchenko3D/zfpmar.h b/marchenko3D/zfpmar.h
new file mode 100644
index 0000000000000000000000000000000000000000..6442357821c1e54a233886c8f02eeb0d377f77f9
--- /dev/null
+++ b/marchenko3D/zfpmar.h
@@ -0,0 +1,65 @@
+#ifndef ZFPMAR_H
+#define ZFPMAR_H
+#define MARBYTES 48
+#define TOPBYTES 80
+
+/* Size of the different types
+long - 8 bytes
+int - 4 bytes
+char - 1 byte
+float - 4 bytes
+*/
+
+/* TYPEDEFS */
+typedef struct { /* zfpmar - headers for the compression and decompression of data for marchenko */
+
+ long nx; /* Number of samples in x-direction */
+
+ long ny; /* Number of samples in y-direction */
+
+ int sx; /* Source coordinate in x-direction */
+
+ int sy; /* Source coordinate in y-direction */
+
+ int sz; /* Source coordinate in z-direction */
+
+ int gx; /* Receiver coordinate in x-direction */
+
+ long gy; /* Receiver coordinate in y-direction */
+
+ long compsize; /* Size of the compressed data */
+
+} zfpmar;
+
+typedef struct { /* zfptop - headers for the compression and decompression of data for marchenko */
+
+ long ndim; /* Number of dimension is between 1 and 4 */
+
+ long nz; /* Number of samples in z-direction */
+
+ long ns; /* Number of shots */
+
+ long nt; /* Number of time samples */
+
+ float dx; /* Sampling distance in x-direction */
+
+ float dy; /* Sampling distance in x-direction */
+
+ float dz; /* Sampling distance in x-direction */
+
+ float scale; /* Scaling of the coordinates and the sampling */
+
+ double tolerance; /* Set the tolerance of the zfp (de)compression */
+
+ float fmin; /* Minimum frequency of the signal */
+
+ float fmax; /* Maximum frequency of the singal */
+
+ float fx; /* First location of receiver in x-direction */
+
+ float fy; /* First location of receiver in y-direction */
+
+ float fz; /* First location of receiver in y-direction */
+
+} zfptop;
+#endif
\ No newline at end of file