lc = 25e3;
lc_fault = 5000;

Fault_length = 40e3;
Fault_width = 20e3;
Fault_dip = 60*Pi/180.;

//Nucleation in X,Z local coordinates
X_nucl = 0e3;
Width_nucl = 0.5*Fault_width;
R_nucl = 2e3;
lc_nucl = 100;

xc = 614520.00;
yc = 5173600.0;
dx = 250e3;
Xmin = xc-dx;
Xmax = xc+dx;
Ymin = yc-dx;
Ymax = yc+dx;
Zmin = -200e3;


//Create the Volume
Point(1) = {Xmin, Ymin, 10e3, lc};
Point(2) = {Xmin, Ymax, 10e3, lc};
Point(3) = {Xmax, Ymax, 10e3, lc};
Point(4) = {Xmax, Ymin, 10e3, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {1,2,3,4};
Plane Surface(1) = {5};
Extrude {0,0, Zmin} { Surface{1}; }

// Epicentre
Point(newp) = {583805.8942, 5170385.0381, -10840.0000, lc_fault};

Macro Fault
  x = X[t];
  y = Y[t];
  depth = -D[t]*1000;
  s = strike[t]*Pi/180.;
  d = dip[t]*Pi/180.;
  r = 0.0;
  l = L[t]*1000.;
  w = -W[t]*1000.;
  sx = Sin(s)*Cos(r) - Cos(s)*Cos(d)*Sin(r);
  sy = Cos(s)*Cos(r) + Sin(s)*Cos(d)*Sin(r);
  sz = Sin(d)*Sin(r);
  dx = -Sin(s)*Sin(r) - Cos(s)*Cos(d)*Cos(r);
  dy = -Cos(s)*Sin(r) + Sin(s)*Cos(d)*Cos(r);
  dz = Sin(d)*Cos(r);

  p1 = newp; Point(p1) = {x, y, depth, lc_fault};
  p2 = newp; Point(p2) = {x+l*sx, y+l*sy, depth+l*sz, lc_fault};
  p3 = newp; Point(p3) = {x+l*sx+w*dx, y+l*sy+w*dy, depth+l*sz+w*dz, lc_fault};
  p4 = newp; Point(p4) = {x+w*dx, y+w*dy, depth+w*dz, lc_fault};

  l1 = newl; Line(l1) = {p1, p2};
  /*If (depth == 0)*/
    /*Line{l1} In Surface{1};*/
  /*EndIf*/
  l2 = newl; Line(l2) = {p2, p3};
  l3 = newl; Line(l3) = {p3, p4};
  l4 = newl; Line(l4) = {p4, p1};

  ll = newll;
  Line Loop(ll) = {l1,l2,l3,l4};
  ps = news;
  Plane Surface(ps) = {ll};

  If (t == 4)
    p5 = newp; Point(p5) = {x + l*sx*7/8+w*dx*3/4, y + l*sy*7/8+w*dy*3/4, depth + l*sz*7/8+w*dz*3/4, lc_fault};
  EndIf
Return

X[] = {614520.00, 604700.00, 586000.00, 603220.00, 588800.00, 576891.69, 630915.38};
Y[] = {5173600.0, 5173000.0, 5172040.0, 5174400.0, 5172650.0, 5177979.0, 5175496.579};
D[] = {0, 0, 0, 1, 1, 0.5, 3.5};
L[] = {12, 20, 12, 8, 8, 8, 19};
W[] = {12, 12, 12, 8, 8, 8, 13};
strike[] = {269, 265, 302, 24, 31, 214, 60};
dip[] = {87, 87, 87, 65, 75, 65, 60};
rake[] = {180, 180, 180, 90, 90, 90, 145};

For t In {0:6}
  Call Fault;
EndFor
