/* Implementation of the formulae of Okada, 1985, "Surface deformation
   due to shear and tensile faults in a half-space", Bulletin of the
   Seismological Society of America, 75:4, 1135-1154, */

/* 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;
}

struct Okada {
  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;
};

void okada (struct Okada p)
{
  // 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) */

  double dtr = pi/180.;
  if (p.rake != nodata) {
    p.vU[0] = p.U*cos (p.rake*dtr);
    p.vU[1] = p.U*sin (p.rake*dtr);
  }
  double sina = sin ((90. - p.strike)*dtr);
  double cosa = cos ((90. - p.strike)*dtr);
  double sind = sin (p.dip*dtr);
  /* depth of the bottom edge */
  double depth = sind > 0. ? p.depth + p.width*sind : p.depth;
  /* origin to the centroid */
  double x0 = p.length/2., y0 = p.width/2.*cos (p.dip*dtr);
  
  foreach() {
    if ( p.flat ) {
      x -= p.x;
      y -= p.y;
    }
    else {
      x = p.R*cos(y*dtr)*dtheta(x, p.x)*dtr;
      y = p.R*dtheta(y, p.y)*dtr;
    }
    double x1 =   cosa*x + sina*y;
    double y1 = - sina*x + cosa*y;
    double oka[3];
    okada_rectangular_source (p.vU, p.length, p.width, depth, 
			      p.dip*dtr,
			      p.mu/p.lambda,
			      x0 + x1, y0 + y1,
			      oka);
    val(p.d,0,0) = oka[2];
  }
}

void fault (struct Okada p)
{
  scalar hold[];
  // save the initial water depth
  scalar_clone (hold, h);
  foreach()
    hold[] = h[];
  boundary ({hold});

  p.d = h;
  do {
    okada (p);
    // 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());
}