This is an implementation of the formulae of Okada, 1985.

/* formulae (25)-(30) */
static void rectangular_source (const double U[3], double cosd, double sind,
double mulambda, double d,
double psi, double eta, double q,
double u[3])
{
double R = sqrt (psi*psi + eta*eta + q*q);
double X = sqrt (psi*psi + q*q);
double dtilde = eta*sind - q*cosd;
double ytilde = eta*cosd + q*sind;
double atanp = fabs (q) > 1e-6 ? atan (psi*eta/(q*R)) : 0.;

mulambda = mulambda/(1. + mulambda);
double logReta = R + eta > 1e-6 ? log (R + eta) : - log (R - eta);
double Reta = fabs (R + eta) > 1e-6 ? R + eta : 1e30;
double I1, I2, I3, I4, I5;
if (fabs (cosd) > 1e-6) {
/* formula (28) */
I5 = fabs (psi) < 1e-6 ? 0. :
mulambda*2./cosd*atan ((eta*(X + q*cosd) +
X*(R + X)*sind)/(psi*(R + X)*cosd));
I4 = mulambda/cosd*(log (R + dtilde) - sind*logReta);
I3 = mulambda*(1./cosd*ytilde/(R + dtilde) - logReta) + sind/cosd*I4;
I2 = mulambda*(- logReta) - I3;
I1 = mulambda*(-1./cosd*psi/(R + dtilde)) - sind/cosd*I5;
}
else {
/* formula (29) */
double R1 = R + dtilde;
I1 = - mulambda/2.*psi*q/(R1*R1);
I3 = mulambda/2.*(eta/R1 + ytilde*q/(R1*R1) - logReta);
I2 = mulambda*(- logReta) - I3;
I4 = - mulambda*q/R1;
I5 = - mulambda*psi*sind/R1;
}

/* strike-slip, formula (25) */
if (U[0] != 0.) {
double U1pi = U[0]/(2.*M_PI);
u[0] -= U1pi*(psi*q/(R*Reta) + atanp + I1*sind);
u[1] -= U1pi*(ytilde*q/(R*Reta) + q*cosd/Reta + I2*sind);
u[2] -= U1pi*(dtilde*q/(R*Reta) + q*sind/Reta + I4*sind);
}

/* dip-slip, formula (26) */
if (U[1] != 0.) {
double U2pi = U[1]/(2.*M_PI);
u[0] -= U2pi*(q/R - I3*sind*cosd);
u[1] -= U2pi*(ytilde*q/(R*(R + psi)) + cosd*atanp - I1*sind*cosd);
u[2] -= U2pi*(dtilde*q/(R*(R + psi)) + sind*atanp - I5*sind*cosd);
}

/* tensile, formula (27) */
if (U[2] != 0.) {
double U3pi = U[2]/(2.*M_PI);
u[0] += U3pi*(q*q/(R*Reta) - I3*sind*sind);
u[1] += U3pi*(-dtilde*q/(R*(R + psi)) -
sind*(psi*q/(R*Reta) - atanp) - I1*sind*sind);
u[2] += U3pi*(ytilde*q/(R*(R + psi)) +
cosd*(psi*q/(R*Reta) - atanp) - I5*sind*sind);
}
}

/* formula (24) */
static void okada_rectangular_source (const double U[3],
double L, double W, double d,
double delta, double mulambda,
double x, double y,
double u[3])
{
double cosd = cos (delta), sind = sin (delta);
double p = y*cosd + d*sind;
double q = y*sind - d*cosd;

u[0] = u[1] = u[2] = 0.;
rectangular_source (U, cosd, sind, mulambda, d,
x, p, q,
u);
rectangular_source (U, cosd, sind, mulambda, d,
x - L, p - W, q,
u);

double u1[3] = {0., 0., 0.};
rectangular_source (U, cosd, sind, mulambda, d,
x, p - W, q,
u1);
rectangular_source (U, cosd, sind, mulambda, d,
x - L, p, q,
u1);
u[0] -= u1[0];
u[1] -= u1[1];
u[2] -= u1[2];
}

static double dtheta (double theta1, double theta2)
{
double d = theta1 - theta2;
if (d > 180.) d -= 360.;
if (d < -180.) d += 360.;
return d;
}

typedef struct {
double depth, x, y;
double strike, dip, rake;
double length, width, U;
double vU[3];
} Fault;

scalar d;
double x, y, depth;
double strike, dip, rake;
double mu, lambda;
double length, width, vU[3], U;
double R;
int (* iterate) (void);
bool flat, centroid;
Fault * faults;
};

{
// default settings
if (p.mu == 0.)     p.mu = 1.;
if (p.lambda == 0.) p.lambda = 1.;
if (p.R == 0.)      p.R = 6371220.; /* Earth radius (metres) */

Fault faults[2] = {0};
if (p.faults == NULL) {
faults[0] = (Fault) {p.depth, p.x, p.y,
p.strike, p.dip, p.rake,
p.length, p.width, p.U};
p.faults = faults;
}
foreach() {
val(p.d) = 0.;
for (Fault * f = p.faults; f && f->depth > 0.; f++) {
double depth = f->depth, dtr = pi/180.;
if (p.centroid)
depth -= f->width*fabs (sin (f->dip*dtr))/2.;
if (f->rake != nodata) {
f->vU[0] = f->U*cos (f->rake*dtr);
f->vU[1] = f->U*sin (f->rake*dtr);
}
double sina = sin ((90. - f->strike)*dtr);
double cosa = cos ((90. - f->strike)*dtr);
double sind = sin (f->dip*dtr);
/* depth of Okada origin */
depth = sind > 0. ? depth + f->width*sind : depth;
/* origin to the centroid */
double x0 = f->length/2., y0 = f->width/2.*cos (f->dip*dtr);

double xc, yc;
if (p.flat) {
xc = x - f->x;
yc = y - f->y;
}
else {
xc = p.R*cos(y*dtr)*dtheta(x, f->x)*dtr;
yc = p.R*dtheta(y, f->y)*dtr;
}
double x1 =   cosa*xc + sina*yc;
double y1 = - sina*xc + cosa*yc;
double oka[3];
f->dip*dtr,
p.mu/p.lambda,
x0 + x1, y0 + y1,
oka);
val(p.d) += oka[2];
}
}
}

## User interface

Use function fault() to alter water depth h (where h > dry) according to the fault parameters:

• x, y: coordinates of the fault centroid (see boolean flat for coordinate type).
• depth: depth of the top edge of the fault (see also centroid).
• strike, dip, rake: fault parameters in degrees. (0 <= strike < 360, -90 <= dip <= 90, -90 <= rake <= 90 where rake = 90 degs and dip > 0 is reverse faulting. NB: Okada defines normal faulting by rake = 90 deg and dip < 0 whereas the seismological convention now is generally 0 <= dip <= 90 and rake = -90 for normal faulting).
• mu, lambda: only the ratio is important and default is mu/lambda = 1.
• length, width, U: length and width of the fault plane and slip on the fault plane (generally in meters).
• vU[3]: displacement vector. Note that vU[0] and vU[1] can be calculated from rake and slip (U), vU[2] is tensile opening of a fault so default is 0. (not needed as long as U and rake are given).
• R: is the radius of the earth (for when x, y are in longitude and latitude i.e. flat = false).
• iterate: is the function to use to iterate.
• flat: true assumes x, y cartesian, false assumes longitude and latitude (default).
• centroid: assumes that depth is measured to the centroid of the fault, not the top edge.
• faults: if non-NULL, this defines an array of several fault parameters which will be used instead of the parameters for a single fault above. This is useful to efficiently define a deformation composed of many Okada subfaults. Note that the array must be terminated by a “dummy fault” of depth smaller than or equal to zero.
void fault (struct Okada p)
{
scalar hold[];
// save the initial water depth
scalar_clone (hold, h);
foreach()
hold[] = h[];

p.d = h;
int nitermax = 20;
do {
// h[] now contains the Okada vertical displacement
foreach() {
// deformation is added to hold[] (water depth) only in wet areas
h[] = hold[] > dry ? max (0., hold[] + h[]) : hold[];
eta[] = zb[] + h[];
}
} while (p.iterate && p.iterate() && nitermax--);
}

## References

 [okada1985] Yoshimitsu Okada. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the seismological society of America, 75(4):1135–1154, 1985.