[cig-commits] r6774 - in cs/cigma/trunk/sandbox: . okada
okada/disloc3d okada/mcginnis
luis at geodynamics.org
luis at geodynamics.org
Thu May 3 16:33:54 PDT 2007
Author: luis
Date: 2007-05-03 16:33:53 -0700 (Thu, 03 May 2007)
New Revision: 6774
Added:
cs/cigma/trunk/sandbox/okada/
cs/cigma/trunk/sandbox/okada/disloc3d/
cs/cigma/trunk/sandbox/okada/disloc3d/Makefile
cs/cigma/trunk/sandbox/okada/disloc3d/ReverseSlip_ng.m
cs/cigma/trunk/sandbox/okada/disloc3d/StrikeSlip_ng.m
cs/cigma/trunk/sandbox/okada/disloc3d/bm5_analytic_README_2.txt
cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.c
cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.h
cs/cigma/trunk/sandbox/okada/disloc3d/dc3d.h
cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.c
cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.h
cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d_mod2.f
cs/cigma/trunk/sandbox/okada/disloc3d/h5array.c
cs/cigma/trunk/sandbox/okada/disloc3d/h5array.h
cs/cigma/trunk/sandbox/okada/disloc3d/h5points.c
cs/cigma/trunk/sandbox/okada/disloc3d/mexfunction.f
cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.c
cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.m
cs/cigma/trunk/sandbox/okada/disloc3d/points0.h5.gz
cs/cigma/trunk/sandbox/okada/disloc3d/points0.in.gz
cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.c
cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.h
cs/cigma/trunk/sandbox/okada/mcginnis/
cs/cigma/trunk/sandbox/okada/mcginnis/furls
cs/cigma/trunk/sandbox/okada/mcginnis/okada.c
cs/cigma/trunk/sandbox/okada/mcginnis/okada.cleanup.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.disp.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.lame.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.notation.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.response.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.self.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.x.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.y.h
cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.z.h
Log:
Okada elastic dislocation solution, original versions taken from
* http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/utilities/elastic-okada
* http://richter.colorado.edu/okada/
Modified version of disloc3d contains a translation of the matlab wrapper into C.
Added: cs/cigma/trunk/sandbox/okada/disloc3d/Makefile
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/Makefile 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/Makefile 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,34 @@
+CC = gcc
+CFLAGS = -Wall -g -O3
+LIBS = -lm -lgfortran -lhdf5
+
+OBJS = \
+ dc3d.o \
+ disloc3d.o \
+ cervelli.o \
+ txtarray.o \
+ h5array.o \
+
+TARGETS = \
+ okadasoln_ng \
+ h5points \
+
+
+all: $(TARGETS) $(OBJS)
+
+dc3d.o: disloc3d_mod2.f
+ gfortran -c $< -o $@
+
+.c.o:
+ $(CC) $(CFLAGS) -c $<
+
+okadasoln_ng: okadasoln_ng.c $(OBJS)
+ $(CC) $(CFLAGS) $^ $(LIBS) -o $@
+
+h5points: h5points.c $(OBJS)
+ $(CC) $(CFLAGS) $^ $(LIBS) -o $@
+
+clean:
+ rm -f $(TARGETS) $(OBJS)
+
+.PHONY: clean
Added: cs/cigma/trunk/sandbox/okada/disloc3d/ReverseSlip_ng.m
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/ReverseSlip_ng.m 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/ReverseSlip_ng.m 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,36 @@
+function ReverseSlip_ng(varargin)
+%
+% Matlab script to generate the analytic boundary conditions for
+% the Reverse-Slip (no gravity) bench mark
+%
+% See okadasoln_ng for more info.
+
+%% Define number of faults to taper with
+taperd = 10.0;
+n = 12000.0:taperd:(16000.0-taperd);
+
+%% Get a number of faults with tapered slip
+fi = [];
+for i = 1 : 1 : length(n)
+ fx1 = 4000.0;
+ fy1 = -n(i);
+
+ fx2 = 4000.0;
+ fy2 = n(i);
+
+ dip = 45.0;
+ D = n(i);
+ Pr = 0.25;
+ bd = 0;
+ ss = 0;
+ ds = 1/length(n);
+ ts = 0;
+
+ fi = [fi; fx1 fy1 fx2 fy2 dip D Pr bd ss ds ts];
+end
+
+if length(varargin)==0
+ okadasoln_ng(fi)
+else
+ okadasoln_ng(fi, varargin{1})
+end
Added: cs/cigma/trunk/sandbox/okada/disloc3d/StrikeSlip_ng.m
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/StrikeSlip_ng.m 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/StrikeSlip_ng.m 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,36 @@
+function StrikeSlip_ng(varargin)
+%
+% Matlab script to generate the analytic boundary conditions for
+% the Strike-Slip (no gravity) bench mark
+%
+% See okadasoln_ng for more info.
+
+%% Define number of faults to taper with
+taperd = 10.0;
+n = 12000.0:taperd:(16000.0-taperd);
+
+%% Get a number of faults with tapered slip
+fi = [];
+for i = 1 : 1 : length(n)
+ fx1 = 12000.0;
+ fy1 = -n(i);
+
+ fx2 = 12000.0;
+ fy2 = n(i);
+
+ dip = 90.0;
+ D = n(i);
+ Pr = 0.25;
+ bd = 0;
+ ss = -1/length(n);
+ ds = 0;
+ ts = 0;
+
+ fi = [fi; fx1 fy1 fx2 fy2 dip D Pr bd ss ds ts];
+end
+
+if length(varargin)==1
+ okadasoln_ng(fi)
+else
+ okadasoln_ng(fi, varargin{1})
+end
Added: cs/cigma/trunk/sandbox/okada/disloc3d/bm5_analytic_README_2.txt
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/bm5_analytic_README_2.txt 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/bm5_analytic_README_2.txt 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,60 @@
+FEM working group Bench Mark # 5 Analytic Okada solutions info:
+
+To FEMers: I don't mind calculating your disp, strain energy etc
+if you send me files with xyz coord triples ( no headers). It will take me
+only a few seconds if they are in the correct format.
+If you want to do these calcs yourself, read on...
+
+
+Here you should find three things:
+ 1. bm5.m - matlab function, call on matlab command line, with no inputs:
+ >> bm5;
+
+ The only output is some files, see the function description at
+ the beginning of bm5.m for file and other info.
+
+ 2. disloc3d_mod2.f - fortran (mex) source code to make mex function.
+ At matlab prompt:
+ >> mex disloc3d_mod2.f
+ This should produce a mex file which you then can call from
+ the command line (bm5.m calls it for you, all you have to do
+ is create the proper mex function for your platform, e.g., for
+ a Sun Solaris, disloc3d_mod2.f --> disloc3d_mod2.mexsol
+
+ Note: mex functions are system and matlab version dependent.
+ That is, a mex function (*.mex*) created on a Sun
+ will not work on an SGI etc, and a mex function created
+ with Matlab6 will not work with matlab5, etc.
+
+ See online matlab mex documentation:
+ http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_external/matlab_external.shtml
+
+
+ 3. dc3d.f - fortran source code (written by Okada) for his 1992 paper.
+ disloced_mod2.f serves as a wrapper around dc3d.f looping over
+ many fault patches and many stations and an interface for matlab
+ mex stuff.
+
+ I got all the fortran and mex stuff downloaded from Peter Cervelli's
+ website (then at Stanford, now at HVO i think). I loooked very quickly
+ and this Stanford site is no longer. I tell you this because i modified
+ Cervelli's original disloc3d.f (into disloc3d_mod2.f) because there was
+ a mistake: it looped before it added instead of adding before looping
+ (an okada geometry detail i can explain in more detail if you want).
+ Moral of the story: discloc3d_mod2.f works correctly.
+
+ matlab function fault_params_to_okada_form.m is from Brendan Meade, MIT
+
+
+
+
+FYI: the fault tapers (in y and z) every 125 meters,
+which makes for 33 fault patches, each having 1/33 m slip.
+This is hardcoded in bm5.m and can be changed by changing variable "taperd".
+
+Contact: Noah Fay
+ nfay at newberry.uoregon.edu
+ (541) 346-4653
+
+
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,247 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "disloc3d.h"
+#include "cervelli.h"
+
+/* function fault_params_to_okada_form
+ *
+ * Based on fault_params_to_okada_form.m MATLAB script
+ * (comments copied verbatim)
+ *
+ * This function takes fault trace, dip and locking depth information
+ * and calculates the anchor coordinates, length, width and strike
+ * of the fault plane as per Okada (1985).
+ *
+ * It may be beneficial in the future to compute the midpoint base
+ * anchor as well for using Okada (1992) Stanford binary.
+ *
+ * We should also allow for buried faults in the future. This can
+ * be achieved by passing a locking depth and a burial depth. The
+ * only output changed would be the width of the fault plane.
+ *
+ * We may need to be more formal about say endpoint1 < endpoint2
+ * ... maybe
+ *
+ * Arguments:
+ * sx1: x coord of fault trace endpoint1
+ * sy1: y coord of fault trace endpoint1
+ * sx2: x coord of fault trace endpoint2
+ * sy2: y coord of fault trace endpoint2
+ * dip: dip of fault plane [radians]
+ * D : fault locking depth
+ * bd : burial depth (top "locking depth")
+ *
+ * Returned variables:
+ * strike: strike of fault plane
+ * L : fault length
+ * W : fault width
+ * ofx : x coord of fault anchor
+ * ofy : y coord of fault anchor
+ * ofxe : x coord of other buried corner
+ * ofye : y coord of other buried corner
+ * tfx : x coord of fault anchor (top relative)
+ * tfy : y coord of fault anchor (top relative)
+ * tfxe : x coord of other buried corner (top relative)
+ * tfye : y coord of other buried corner (top relative)
+ */
+static void
+fault_params_to_okada_form(double sx1, double sy1,
+ double sx2, double sy2,
+ double dip, double D, double bd,
+ double *strike,
+ double *L, double *W,
+ double *ofx, double *ofy,
+ double *ofxe, double *ofye,
+ double *tfx, double *tfy,
+ double *tfxe, double *tfye)
+{
+ //
+ // Do some basic checks
+ //
+ if (D < bd)
+ {
+ fprintf(stderr, "Burial depth is greater than locking depth!\n");
+ fprintf(stderr, "Fault width will be negative!\n");
+ }
+ if (D == bd)
+ {
+ fprintf(stderr, "Burial depth is equal to locking depth!\n");
+ fprintf(stderr, "Fault width will be zero!\n");
+ }
+
+ //
+ // Calculate needed parameters for Okada input
+ //
+ *strike = atan2(sy1 - sy2, sx1 - sx2) + M_PI;
+ *L = sqrt((sx2 - sx1)*(sx2 - sx1) + (sy2 - sy1)*(sy2 - sy1));
+ *W = (D - bd) / sin(dip);
+
+ //
+ // For older versions without a burial depth option
+ //
+ // *W = D/sin(dip)
+ //
+
+ //
+ // Calculate fault segment anchor and other buried point
+ //
+ *ofx = sx1 + D / tan(dip) + sin(*strike);
+ *ofy = sy1 - D / tan(dip) + cos(*strike);
+ *ofxe = sx2 + D / tan(dip) + sin(*strike);
+ *ofye = sy2 - D / tan(dip) + cos(*strike);
+
+ //
+ // Calculate fault segment anchor and other buried point (top relative)
+ //
+ *tfx = sx1 + bd / tan(dip) + sin(*strike);
+ *tfy = sy1 - bd / tan(dip) + cos(*strike);
+ *tfxe = sx2 + bd / tan(dip) + sin(*strike);
+ *tfye = sy2 - bd / tan(dip) + cos(*strike);
+}
+
+
+
+/* function getM
+ * Convert fault parameters into form needed for Cervelli calculations.
+ *
+ * Inputs: fi fault info
+ * nmod number of fault patches
+ * mu Poisson's ratio
+ *
+ * Outputs: M matrix with fault stuff
+ * (see below for what is on the rows)
+ * (one row per fault patch)
+ */
+static void
+getM(double *fi, int nmod, double mu, double *M)
+{
+ int j;
+ double *fault;
+ double *model;
+
+ double fx1,fy1;
+ double fx2,fy2;
+ double dip;
+ double D;
+ double Pr;
+ double bd; // burial depth
+ double ss;
+ double ds;
+ double ts;
+
+ double strike;
+ double L, W;
+ double ofx, ofy;
+ double ofxe, ofye;
+ double tfx, tfy;
+ double tfxe, tfye;
+
+ double length;
+ double width;
+ double depth;
+ double str;
+ double east;
+ double north;
+
+ for (j = 1; j < nmod; j++)
+ {
+ fault = &fi[11*j];
+ model = &M[10*j];
+
+ fx1 = fault[0];
+ fy1 = fault[1];
+ fx2 = fault[2];
+ fy2 = fault[3];
+ dip = fault[4]; // in degs here
+ D = fault[5];
+ Pr = fault[6];
+ bd = fault[7];
+ ss = fault[8];
+ ds = fault[9];
+ ts = fault[10];
+
+ fault_params_to_okada_form(
+ fx1, fy1, fx2, fy2, dip*M_PI/180.0, D, bd,
+ &strike, &L, &W, &ofx, &ofy, &ofxe, &ofye, &tfx, &tfy, &tfxe, &tfye);
+
+ length = L;
+ width = W;
+
+ if (dip < 0)
+ {
+ east = (L/2) * cos(strike) + tfx; // east component of top center
+ north = (L/2) * sin(strike) + tfy; // north " " " "
+ depth = bd;
+ width = abs(width);
+ ss = -1 * ss;
+ ds = -1 * ds;
+ ts = -1 * ts;
+ }
+ else
+ {
+ east = (L/2) * cos(strike) + ofx; // east component of bottom center
+ north = (L/2) * sin(strike) + ofy;
+ depth = D;
+ }
+
+ // Cervelli takes strike in degs (strike comes out of
+ // fault_params_to_okada_form in rads)
+ str = (180.0/M_PI) * ((M_PI/2) - strike);
+
+ model[0] = length;
+ model[1] = width;
+ model[2] = depth;
+ model[3] = dip;
+ model[4] = str;
+ model[5] = east;
+ model[6] = north;
+ model[7] = ss;
+ model[8] = ds;
+ model[9] = ts;
+ }
+
+}
+
+
+
+/* function calc_disp_cervelli
+ *
+ * Calculate displacement's via Cervelli's Okada (1992) code
+ *
+ * Inputs: fault_info matrix with fault stuff
+ * sx station x coord
+ * sy station y coord
+ * sz station z coord
+ *
+ * Outputs: dispx disp in x direction
+ * dispy disp in y direction
+ * dispz disp in z direction
+ * D matrix with derivatives
+ * S Stresses (not used)
+ *
+ * TODO: update this comment
+ */
+void calc_disp_cervelli(double mu, double nu,
+ double *models, double *fault_info, int nfaults,
+ double *stations, int nstations,
+ double *disp, double *deriv, double *stress,
+ double *flag, double *flag2)
+{
+ //
+ // get the "model matrix" with
+ // fault patch info in the form for cervelli calcs
+ //
+ getM(fault_info, nfaults, mu, models);
+
+ //
+ // call the disloc3d wrapper
+ //
+ disloc3d(models, nfaults,
+ stations, nstations,
+ mu, nu,
+ disp, deriv, stress,
+ flag, flag2);
+
+}
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/cervelli.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,11 @@
+#ifndef __CERVELLI_H__
+#define __CERVELLI_H__
+
+
+void calc_disp_cervelli(double mu, double nu,
+ double *models, double *fault_info, int nfaults,
+ double *stations, int nstations,
+ double *disp, double *deriv, double *stress,
+ double *flag, double *flag2);
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/disloc3d/dc3d.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/dc3d.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/dc3d.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,16 @@
+#ifndef __DC3D_H__
+#define __DC3D_H__
+
+void dc3d_(double *alpha,
+ double *x, double *y, double *z,
+ double *depth, double *dip,
+ double *al1, double *al2,
+ double *aw1, double *aw2,
+ double *disl1, double *disl2, double *disl3,
+ double *ux, double *uy, double *uz,
+ double *uxx, double *uyx, double *uzx,
+ double *uxy, double *uyy, double *uzy,
+ double *uxz, double *uyz, double *uzz,
+ double *iret);
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,191 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "dc3d.h"
+#include "disloc3d.h"
+
+//
+// TODO: change flagout,flagout2 to int arrays
+// TODO: flagout2 currently disabled
+//
+
+void disloc3d(double *models, int nmod,
+ double *stations, int nstat,
+ double mu, double nu,
+ double *uout, double *dout, double *sout,
+ double *flagout, double *flagout2)
+{
+ // input arguments
+ // models [nmod x 10]
+ // stations [nstat x 3]
+ // mu [1 x 1]
+ // nu [1 x 1]
+
+ // matrices for return arguments
+ // uout [nstat x 3]
+ // dout [nstat x 9]
+ // sout [nstat x 6]
+ // flagout [nstat x 1]
+ // flagout2 [nstat x nmod]
+
+ double lambda;
+ double theta;
+
+ double *model;
+ double *stat;
+ double *u, *d, *s;
+ double *flag;
+ //double *flag2;
+ double iret;
+
+ double ux, uy, uz;
+ double uxx, uxy, uxz;
+ double uyx, uyy, uyz;
+ double uzx, uzy, uzz;
+
+ double uxt, uyt, uzt;
+ double uxxt, uxyt, uxzt;
+ double uyxt, uyyt, uyzt;
+ double uzxt, uzyt, uzzt;
+
+ double alpha;
+ double x, y, z;
+
+ double strike;
+ double cs,ss;
+
+ double dip;
+ double cd,sd;
+
+ double depth;
+ double disl1, disl2, disl3;
+ double al1, al2;
+ double aw1, aw2;
+
+ int i, j;
+
+
+ lambda = 2*mu*nu/(1 - 2*nu);
+ alpha = (lambda + mu)/(lambda + 2*mu);
+
+ for (i = 0; i < nstat; i++)
+ {
+ stat = &stations[3*i];
+ flag = &flagout[i];
+ //flag2 = &flagout2[nmod*i];
+
+ *flag = 0;
+
+ if (stat[2] > 0)
+ {
+ *flag += 1;
+ fprintf(stderr, "Warning: Positive depth given. (station %d)\n", i);
+ }
+ else
+ {
+ uxt = uyt = uzt = 0;
+ uxxt = uxyt = uxzt = 0;
+ uyxt = uyyt = uyzt = 0;
+ uzxt = uzyt = uzzt = 0;
+
+ for (j = 0; j < nmod; j++)
+ {
+ model = &models[10*j];
+
+ strike = (model[4] - 90)*(M_PI/180.0);
+ cs = cos(strike);
+ ss = sin(strike);
+ dip = model[3];
+ cd = cos(dip * M_PI/180.0);
+ sd = sin(dip * M_PI/180.0);
+ disl1 = model[7];
+ disl2 = model[8];
+ disl3 = model[9];
+ depth = model[2] - 0.5*model[1]*sd;
+ al1 = model[0]/2;
+ al2 = al1;
+ aw1 = model[1]/2;
+ aw2 = aw1;
+ x = cs * (-model[5] + stat[0]) - ss * (-model[6] + stat[1]);
+ y = ss * (-model[5] + stat[0]) + cs * (-model[6] + stat[1]) - 0.5 * cd * model[1];
+ z = stat[2];
+
+ // generate warnings for unphysical models
+ if ((model[2] - sd*model[1] < 0) ||
+ (model[0] <= 0) ||
+ (model[1] <= 0) ||
+ (model[2] < 0))
+ {
+ *flag += 1;
+ fprintf(stderr, "Warning: Unphysical model (station %d, subfault %d)\n",i,j);
+ }
+ else
+ {
+ dc3d_(&alpha, &x, &y, &z, &depth, &dip,
+ &al1, &al2, &aw1, &aw2, &disl1, &disl2, &disl3,
+ &ux, &uy, &uz,
+ &uxx, &uyx, &uzx,
+ &uxy, &uyy, &uzy,
+ &uxz, &uyz, &uzz,
+ &iret);
+
+ // fill flags - the rows contain flags for each station, one row per model/fault
+ *flag += iret;
+ //flag2[j] = iret;
+
+ //
+ // rotate then add
+ //
+ uxt += cs*ux + ss*uy;
+ uyt += -ss*ux + cs*uy;
+ uzt += uz;
+
+ uxxt += cs*cs*uxx + cs*ss*(uxy + uyx) + ss*ss*uyy; // aaa
+ uxyt += cs*cs*uxy - ss*ss*uyx + cs*ss*(-uxx + uyy); // bbb
+ uxzt += cs*uxz + ss*uyz; // ccc
+
+ uyxt += -(ss*(cs*uxx + ss*uxy)) + cs*(cs*uyx + ss*uyy); // ddd
+ uyyt += ss*ss*uxx - cs*ss*(uxy + uyx) + cs*cs*uyy; // eee
+ uyzt += -(ss*uxz) + cs*uyz; // fff
+
+ uzxt += cs*uzx + ss*uzy; // ggg
+ uzyt += -(ss*uzx) + cs*uzy; // hhh
+ uzzt += uzz; // iii
+ }
+ }
+
+ // jump to appropriate offset
+ u = &uout[3*i];
+ d = &dout[9*i];
+ s = &sout[6*i];
+
+ // assign outputs
+ u[0] = uxt;
+ u[1] = uyt;
+ u[2] = uzt;
+
+ d[0] = uxxt;
+ d[1] = uxyt;
+ d[2] = uxzt;
+
+ d[3] = uyxt;
+ d[4] = uyyt;
+ d[5] = uyzt;
+
+ d[6] = uzxt;
+ d[7] = uzyt;
+ d[8] = uzzt;
+
+ // calculate stresses
+ theta = d[0] + d[4] + d[8];
+ s[0] = lambda*theta + 2*mu*d[0];
+ s[1] = mu*(d[1] + d[3]);
+ s[2] = mu*(d[2] + d[6]);
+ s[3] = lambda*theta + 2*mu*d[4];
+ s[4] = mu*(d[5] + d[7]);
+ s[5] = lambda*theta + 2*mu*d[8];
+
+ } // loop over models (subfaults)
+
+ } // loop over stations
+}
Added: cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,10 @@
+#ifndef __DISLOC3D_H__
+#define __DISLOC3D_H__
+
+void disloc3d(double *models, int nmod,
+ double *stations, int nstat,
+ double mu, double nu,
+ double *uout, double *dout, double *sout,
+ double *flag, double *flag2);
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d_mod2.f
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d_mod2.f 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/disloc3d_mod2.f 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,612 @@
+ SUBROUTINE DC3D(ALPHA,X,Y,Z,DEPTH,DIP, 04610000
+ * AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, 04620002
+ * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET) 04630002
+ IMPLICIT REAL*8 (A-H,O-Z) 04640000
+ REAL*8 ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, 04650000
+ * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET 04660000
+C 04670000
+C******************************************************************** 04680000
+C***** ***** 04690000
+C***** DISPLACEMENT AND STRAIN AT DEPTH ***** 04700000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 04710000
+C***** CODED BY Y.OKADA ... SEP 1991 ***** 04720002
+C***** REVISED Y.OKADA ... NOV 1991 ***** 04730002
+C***** ***** 04740000
+C******************************************************************** 04750000
+C 04760000
+C***** INPUT 04770000
+C***** ALPHA : MEDIUM CONSTANT (LAMBDA+MYU)/(LAMBDA+2*MYU) 04780000
+C***** X,Y,Z : COORDINATE OF OBSERVING POINT 04790000
+C***** DEPTH : SOURCE DEPTH 04800000
+C***** DIP : DIP-ANGLE (DEGREE) 04810000
+C***** AL1,AL2 : FAULT LENGTH (-STRIKE,+STRIKE) 04820000
+C***** AW1,AW2 : FAULT WIDTH ( DOWNDIP, UPDIP) 04830000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 04840000
+C 04850000
+C***** OUTPUT 04860000
+C***** UX, UY, UZ : DISPLACEMENT ( UNIT=(UNIT OF DISL) 04870000
+C***** UXX,UYX,UZX : X-DERIVATIVE ( UNIT=(UNIT OF DISL) / 04880000
+C***** UXY,UYY,UZY : Y-DERIVATIVE (UNIT OF X,Y,Z,DEPTH,AL,AW) )04890000
+C***** UXZ,UYZ,UZZ : Z-DERIVATIVE 04900000
+C***** IRET : RETURN CODE ( =0....NORMAL, =1....SINGULAR ) 04910002
+C 04920000
+ COMMON /C0/DUMMY(5),SD,CD ,Dummy2(5)
+ COMMON /C2/XI2,ET2,Q2,R ,Dummy3(20)
+ DIMENSION U(12),DU(12),DUA(12),DUB(12),DUC(12) 04950000
+ DATA F0/0.D0/ 04960000
+C----- 04970000
+ IF(Z.GT.0.) WRITE(6,'('' ** POSITIVE Z WAS GIVEN IN SUB-DC3D'')') 04980000
+ DO 111 I=1,12 04990000
+ U (I)=F0 05000000
+ DUA(I)=F0 05010000
+ DUB(I)=F0 05020000
+ DUC(I)=F0 05030000
+ 111 CONTINUE 05040000
+ AALPHA=ALPHA 05050000
+ DDIP=DIP 05060000
+ CALL DCCON0(AALPHA,DDIP) 05070000
+C====================================== 05080000
+C===== REAL-SOURCE CONTRIBUTION ===== 05090000
+C====================================== 05100000
+ D=DEPTH+Z 05110000
+ P=Y*CD+D*SD 05120000
+ Q=Y*SD-D*CD 05130000
+ JXI=0 05140000
+ JET=0 05150000
+ IF((X+AL1)*(X-AL2).LE.0.) JXI=1 05160000
+ IF((P+AW1)*(P-AW2).LE.0.) JET=1 05170000
+ DD1=DISL1 05180000
+ DD2=DISL2 05190000
+ DD3=DISL3 05200000
+C----- 05210000
+ DO 223 K=1,2 05220000
+ IF(K.EQ.1) ET=P+AW1 05230000
+ IF(K.EQ.2) ET=P-AW2 05240000
+ DO 222 J=1,2 05250000
+ IF(J.EQ.1) XI=X+AL1 05260000
+ IF(J.EQ.2) XI=X-AL2 05270000
+ CALL DCCON2(XI,ET,Q,SD,CD) 05280002
+ IF(JXI.EQ.1 .AND. Q.EQ.F0 .AND. ET.EQ.F0) GO TO 99 05290002
+ IF(JET.EQ.1 .AND. Q.EQ.F0 .AND. XI.EQ.F0) GO TO 99 05300002
+ CALL UA(XI,ET,Q,DD1,DD2,DD3,DUA) 05310000
+C----- 05320000
+ DO 220 I=1,10,3 05330000
+ DU(I) =-DUA(I) 05340000
+ DU(I+1)=-DUA(I+1)*CD+DUA(I+2)*SD 05350000
+ DU(I+2)=-DUA(I+1)*SD-DUA(I+2)*CD 05360000
+ IF(I.LT.10) GO TO 220 05370000
+ DU(I) =-DU(I) 05380000
+ DU(I+1)=-DU(I+1) 05390000
+ DU(I+2)=-DU(I+2) 05400000
+ 220 CONTINUE 05410000
+ DO 221 I=1,12 05420000
+ IF(J+K.NE.3) U(I)=U(I)+DU(I) 05430000
+ IF(J+K.EQ.3) U(I)=U(I)-DU(I) 05440000
+ 221 CONTINUE 05450000
+C----- 05460000
+ 222 CONTINUE 05470000
+ 223 CONTINUE 05480000
+C======================================= 05490000
+C===== IMAGE-SOURCE CONTRIBUTION ===== 05500000
+C======================================= 05510000
+ ZZ=Z 05520000
+ D=DEPTH-Z 05530000
+ P=Y*CD+D*SD 05540000
+ Q=Y*SD-D*CD 05550000
+ JET=0 05560000
+ IF((P+AW1)*(P-AW2).LE.0.) JET=1 05570000
+C----- 05580000
+ DO 334 K=1,2 05590000
+ IF(K.EQ.1) ET=P+AW1 05600000
+ IF(K.EQ.2) ET=P-AW2 05610000
+ DO 333 J=1,2 05620000
+ IF(J.EQ.1) XI=X+AL1 05630000
+ IF(J.EQ.2) XI=X-AL2 05640000
+ CALL DCCON2(XI,ET,Q,SD,CD) 05650002
+ CALL UA(XI,ET,Q,DD1,DD2,DD3,DUA) 05660000
+ CALL UB(XI,ET,Q,DD1,DD2,DD3,DUB) 05670000
+ CALL UC(XI,ET,Q,ZZ,DD1,DD2,DD3,DUC) 05680000
+C----- 05690000
+ DO 330 I=1,10,3 05700000
+ DU(I)=DUA(I)+DUB(I)+Z*DUC(I) 05710000
+ DU(I+1)=(DUA(I+1)+DUB(I+1)+Z*DUC(I+1))*CD 05720000
+ * -(DUA(I+2)+DUB(I+2)+Z*DUC(I+2))*SD 05730000
+ DU(I+2)=(DUA(I+1)+DUB(I+1)-Z*DUC(I+1))*SD 05740000
+ * +(DUA(I+2)+DUB(I+2)-Z*DUC(I+2))*CD 05750000
+ IF(I.LT.10) GO TO 330 05760000
+ DU(10)=DU(10)+DUC(1) 05770000
+ DU(11)=DU(11)+DUC(2)*CD-DUC(3)*SD 05780000
+ DU(12)=DU(12)-DUC(2)*SD-DUC(3)*CD 05790000
+ 330 CONTINUE 05800000
+ DO 331 I=1,12 05810000
+ IF(J+K.NE.3) U(I)=U(I)+DU(I) 05820000
+ IF(J+K.EQ.3) U(I)=U(I)-DU(I) 05830000
+ 331 CONTINUE 05840000
+C----- 05850000
+ 333 CONTINUE 05860000
+ 334 CONTINUE 05870000
+C===== 05880000
+ UX=U(1) 05890000
+ UY=U(2) 05900000
+ UZ=U(3) 05910000
+ UXX=U(4) 05920000
+ UYX=U(5) 05930000
+ UZX=U(6) 05940000
+ UXY=U(7) 05950000
+ UYY=U(8) 05960000
+ UZY=U(9) 05970000
+ UXZ=U(10) 05980000
+ UYZ=U(11) 05990000
+ UZZ=U(12) 06000000
+ IRET=0 06010002
+ RETURN 06020000
+C======================================= 06030000
+C===== IN CASE OF SINGULAR (R=0) ===== 06040000
+C======================================= 06050000
+ 99 UX=F0 06060000
+ UY=F0 06070000
+ UZ=F0 06080000
+ UXX=F0 06090000
+ UYX=F0 06100000
+ UZX=F0 06110000
+ UXY=F0 06120000
+ UYY=F0 06130000
+ UZY=F0 06140000
+ UXZ=F0 06150000
+ UYZ=F0 06160000
+ UZZ=F0 06170000
+ IRET=1 06180002
+ RETURN 06190000
+ END 06200000
+
+
+
+ SUBROUTINE DCCON0(ALPHA,DIP) 09300000
+ IMPLICIT REAL*8 (A-H,O-Z) 09310000
+C 09320000
+C******************************************************************* 09330000
+C***** CALCULATE MEDIUM CONSTANTS AND FAULT-DIP CONSTANTS ***** 09340000
+C******************************************************************* 09350000
+C 09360000
+C***** INPUT 09370000
+C***** ALPHA : MEDIUM CONSTANT (LAMBDA+MYU)/(LAMBDA+2*MYU) 09380000
+C***** DIP : DIP-ANGLE (DEGREE) 09390000
+C### CAUTION ### IF COS(DIP) IS SUFFICIENTLY SMALL, IT IS SET TO ZERO 09400000
+C 09410000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 09420000
+c everything in CO is used in this routine
+
+ DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/ 09430000
+ DATA EPS/1.D-6/ 09440000
+C----- 09450000
+ ALP1=(F1-ALPHA)/F2 09460000
+ ALP2= ALPHA/F2 09470000
+ ALP3=(F1-ALPHA)/ALPHA 09480000
+ ALP4= F1-ALPHA 09490000
+ ALP5= ALPHA 09500000
+C----- 09510000
+ P18=PI2/360.D0 09520000
+ SD=DSIN(DIP*P18) 09530000
+ CD=DCOS(DIP*P18) 09540000
+ IF(DABS(CD).LT.EPS) THEN 09550000
+ CD=F0 09560000
+ IF(SD.GT.F0) SD= F1 09570000
+ IF(SD.LT.F0) SD=-F1 09580000
+ ENDIF 09590000
+ SDSD=SD*SD 09600000
+ CDCD=CD*CD 09610000
+ SDCD=SD*CD 09620000
+ S2D=F2*SDCD 09630000
+ C2D=CDCD-SDSD 09640000
+ RETURN 09650000
+ END 09660000
+
+
+ SUBROUTINE UA(XI,ET,Q,DISL1,DISL2,DISL3,U) 06210000
+ IMPLICIT REAL*8 (A-H,O-Z) 06220000
+ DIMENSION U(12),DU(12) 06230000
+C 06240000
+C******************************************************************** 06250000
+C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-A) ***** 06260000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 06270000
+C******************************************************************** 06280000
+C 06290000
+C***** INPUT 06300000
+C***** XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM 06310000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 06320000
+C***** OUTPUT 06330000
+C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES 06340000
+C 06350000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 06360000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 06370000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 06380000
+ DATA F0,F2,PI2/0.D0,2.D0,6.283185307179586D0/ 06390000
+C----- 06400000
+ DO 111 I=1,12 06410000
+ 111 U(I)=F0 06420000
+ XY=XI*Y11 06430000
+ QX=Q *X11 06440000
+ QY=Q *Y11 06450000
+C====================================== 06460000
+C===== STRIKE-SLIP CONTRIBUTION ===== 06470000
+C====================================== 06480000
+ IF(DISL1.NE.F0) THEN 06490000
+ DU( 1)= TT/F2 +ALP2*XI*QY 06500000
+ DU( 2)= ALP2*Q/R 06510000
+ DU( 3)= ALP1*ALE -ALP2*Q*QY 06520000
+ DU( 4)=-ALP1*QY -ALP2*XI2*Q*Y32 06530000
+ DU( 5)= -ALP2*XI*Q/R3 06540000
+ DU( 6)= ALP1*XY +ALP2*XI*Q2*Y32 06550000
+ DU( 7)= ALP1*XY*SD +ALP2*XI*FY+D/F2*X11 06560000
+ DU( 8)= ALP2*EY 06570000
+ DU( 9)= ALP1*(CD/R+QY*SD) -ALP2*Q*FY 06580000
+ DU(10)= ALP1*XY*CD +ALP2*XI*FZ+Y/F2*X11 06590000
+ DU(11)= ALP2*EZ 06600000
+ DU(12)=-ALP1*(SD/R-QY*CD) -ALP2*Q*FZ 06610000
+ DO 222 I=1,12 06620000
+ 222 U(I)=U(I)+DISL1/PI2*DU(I) 06630000
+ ENDIF 06640000
+C====================================== 06650000
+C===== DIP-SLIP CONTRIBUTION ===== 06660000
+C====================================== 06670000
+ IF(DISL2.NE.F0) THEN 06680000
+ DU( 1)= ALP2*Q/R 06690000
+ DU( 2)= TT/F2 +ALP2*ET*QX 06700000
+ DU( 3)= ALP1*ALX -ALP2*Q*QX 06710000
+ DU( 4)= -ALP2*XI*Q/R3 06720000
+ DU( 5)= -QY/F2 -ALP2*ET*Q/R3 06730000
+ DU( 6)= ALP1/R +ALP2*Q2/R3 06740000
+ DU( 7)= ALP2*EY 06750000
+ DU( 8)= ALP1*D*X11+XY/F2*SD +ALP2*ET*GY 06760000
+ DU( 9)= ALP1*Y*X11 -ALP2*Q*GY 06770000
+ DU(10)= ALP2*EZ 06780000
+ DU(11)= ALP1*Y*X11+XY/F2*CD +ALP2*ET*GZ 06790000
+ DU(12)=-ALP1*D*X11 -ALP2*Q*GZ 06800000
+ DO 333 I=1,12 06810000
+ 333 U(I)=U(I)+DISL2/PI2*DU(I) 06820000
+ ENDIF 06830000
+C======================================== 06840000
+C===== TENSILE-FAULT CONTRIBUTION ===== 06850000
+C======================================== 06860000
+ IF(DISL3.NE.F0) THEN 06870000
+ DU( 1)=-ALP1*ALE -ALP2*Q*QY 06880000
+ DU( 2)=-ALP1*ALX -ALP2*Q*QX 06890000
+ DU( 3)= TT/F2 -ALP2*(ET*QX+XI*QY) 06900000
+ DU( 4)=-ALP1*XY +ALP2*XI*Q2*Y32 06910000
+ DU( 5)=-ALP1/R +ALP2*Q2/R3 06920000
+ DU( 6)=-ALP1*QY -ALP2*Q*Q2*Y32 06930000
+ DU( 7)=-ALP1*(CD/R+QY*SD) -ALP2*Q*FY 06940000
+ DU( 8)=-ALP1*Y*X11 -ALP2*Q*GY 06950000
+ DU( 9)= ALP1*(D*X11+XY*SD) +ALP2*Q*HY 06960000
+ DU(10)= ALP1*(SD/R-QY*CD) -ALP2*Q*FZ 06970000
+ DU(11)= ALP1*D*X11 -ALP2*Q*GZ 06980000
+ DU(12)= ALP1*(Y*X11+XY*CD) +ALP2*Q*HZ 06990000
+ DO 444 I=1,12 07000000
+ 444 U(I)=U(I)+DISL3/PI2*DU(I) 07010000
+ ENDIF 07020000
+ RETURN 07030000
+ END 07040000
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ SUBROUTINE UB(XI,ET,Q,DISL1,DISL2,DISL3,U) 07050000
+ IMPLICIT REAL*8 (A-H,O-Z) 07060000
+ DIMENSION U(12),DU(12) 07070000
+C 07080000
+C******************************************************************** 07090000
+C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-B) ***** 07100000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 07110000
+C******************************************************************** 07120000
+C 07130000
+C***** INPUT 07140000
+C***** XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM 07150000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 07160000
+C***** OUTPUT 07170000
+C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES 07180000
+C 07190000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 07200000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 07210000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 07220000
+ DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/ 07230000
+C----- 07240000
+ RD=R+D 07250000
+ D11=F1/(R*RD) 07260000
+ AJ2=XI*Y/RD*D11 07270000
+ AJ5=-(D+Y*Y/RD)*D11 07280000
+ IF(CD.NE.F0) THEN 07290000
+ IF(XI.EQ.F0) THEN 07300000
+ AI4=F0 07310000
+ ELSE 07320000
+ X=DSQRT(XI2+Q2) 07330000
+ AI4=F1/CDCD*( XI/RD*SDCD 07340000
+ * +F2*DATAN((ET*(X+Q*CD)+X*(R+X)*SD)/(XI*(R+X)*CD)) ) 07350000
+ ENDIF 07360000
+ AI3=(Y*CD/RD-ALE+SD*DLOG(RD))/CDCD 07370000
+ AK1=XI*(D11-Y11*SD)/CD 07380000
+ AK3=(Q*Y11-Y*D11)/CD 07390000
+ AJ3=(AK1-AJ2*SD)/CD 07400000
+ AJ6=(AK3-AJ5*SD)/CD 07410000
+ ELSE 07420000
+ RD2=RD*RD 07430000
+ AI3=(ET/RD+Y*Q/RD2-ALE)/F2 07440000
+ AI4=XI*Y/RD2/F2 07450000
+ AK1=XI*Q/RD*D11 07460000
+ AK3=SD/RD*(XI2*D11-F1) 07470000
+ AJ3=-XI/RD2*(Q2*D11-F1/F2) 07480000
+ AJ6=-Y/RD2*(XI2*D11-F1/F2) 07490000
+ ENDIF 07500000
+C----- 07510000
+ XY=XI*Y11 07520000
+ AI1=-XI/RD*CD-AI4*SD 07530000
+ AI2= DLOG(RD)+AI3*SD 07540000
+ AK2= F1/R+AK3*SD 07550000
+ AK4= XY*CD-AK1*SD 07560000
+ AJ1= AJ5*CD-AJ6*SD 07570000
+ AJ4=-XY-AJ2*CD+AJ3*SD 07580000
+C===== 07590000
+ DO 111 I=1,12 07600000
+ 111 U(I)=F0 07610000
+ QX=Q*X11 07620000
+ QY=Q*Y11 07630000
+C====================================== 07640000
+C===== STRIKE-SLIP CONTRIBUTION ===== 07650000
+C====================================== 07660000
+ IF(DISL1.NE.F0) THEN 07670000
+ DU( 1)=-XI*QY-TT -ALP3*AI1*SD 07680000
+ DU( 2)=-Q/R +ALP3*Y/RD*SD 07690000
+ DU( 3)= Q*QY -ALP3*AI2*SD 07700000
+ DU( 4)= XI2*Q*Y32 -ALP3*AJ1*SD 07710000
+ DU( 5)= XI*Q/R3 -ALP3*AJ2*SD 07720000
+ DU( 6)=-XI*Q2*Y32 -ALP3*AJ3*SD 07730000
+ DU( 7)=-XI*FY-D*X11 +ALP3*(XY+AJ4)*SD 07740000
+ DU( 8)=-EY +ALP3*(F1/R+AJ5)*SD 07750000
+ DU( 9)= Q*FY -ALP3*(QY-AJ6)*SD 07760000
+ DU(10)=-XI*FZ-Y*X11 +ALP3*AK1*SD 07770000
+ DU(11)=-EZ +ALP3*Y*D11*SD 07780000
+ DU(12)= Q*FZ +ALP3*AK2*SD 07790000
+ DO 222 I=1,12 07800000
+ 222 U(I)=U(I)+DISL1/PI2*DU(I) 07810000
+ ENDIF 07820000
+C====================================== 07830000
+C===== DIP-SLIP CONTRIBUTION ===== 07840000
+C====================================== 07850000
+ IF(DISL2.NE.F0) THEN 07860000
+ DU( 1)=-Q/R +ALP3*AI3*SDCD 07870000
+ DU( 2)=-ET*QX-TT -ALP3*XI/RD*SDCD 07880000
+ DU( 3)= Q*QX +ALP3*AI4*SDCD 07890000
+ DU( 4)= XI*Q/R3 +ALP3*AJ4*SDCD 07900000
+ DU( 5)= ET*Q/R3+QY +ALP3*AJ5*SDCD 07910000
+ DU( 6)=-Q2/R3 +ALP3*AJ6*SDCD 07920000
+ DU( 7)=-EY +ALP3*AJ1*SDCD 07930000
+ DU( 8)=-ET*GY-XY*SD +ALP3*AJ2*SDCD 07940000
+ DU( 9)= Q*GY +ALP3*AJ3*SDCD 07950000
+ DU(10)=-EZ -ALP3*AK3*SDCD 07960000
+ DU(11)=-ET*GZ-XY*CD -ALP3*XI*D11*SDCD 07970000
+ DU(12)= Q*GZ -ALP3*AK4*SDCD 07980000
+ DO 333 I=1,12 07990000
+ 333 U(I)=U(I)+DISL2/PI2*DU(I) 08000000
+ ENDIF 08010000
+C======================================== 08020000
+C===== TENSILE-FAULT CONTRIBUTION ===== 08030000
+C======================================== 08040000
+ IF(DISL3.NE.F0) THEN 08050000
+ DU( 1)= Q*QY -ALP3*AI3*SDSD 08060000
+ DU( 2)= Q*QX +ALP3*XI/RD*SDSD 08070000
+ DU( 3)= ET*QX+XI*QY-TT -ALP3*AI4*SDSD 08080000
+ DU( 4)=-XI*Q2*Y32 -ALP3*AJ4*SDSD 08090000
+ DU( 5)=-Q2/R3 -ALP3*AJ5*SDSD 08100000
+ DU( 6)= Q*Q2*Y32 -ALP3*AJ6*SDSD 08110000
+ DU( 7)= Q*FY -ALP3*AJ1*SDSD 08120000
+ DU( 8)= Q*GY -ALP3*AJ2*SDSD 08130000
+ DU( 9)=-Q*HY -ALP3*AJ3*SDSD 08140000
+ DU(10)= Q*FZ +ALP3*AK3*SDSD 08150000
+ DU(11)= Q*GZ +ALP3*XI*D11*SDSD 08160000
+ DU(12)=-Q*HZ +ALP3*AK4*SDSD 08170000
+ DO 444 I=1,12 08180000
+ 444 U(I)=U(I)+DISL3/PI2*DU(I) 08190000
+ ENDIF 08200000
+ RETURN 08210000
+ END 08220000
+
+
+
+
+
+ SUBROUTINE UC(XI,ET,Q,Z,DISL1,DISL2,DISL3,U) 08230000
+ IMPLICIT REAL*8 (A-H,O-Z) 08240000
+ DIMENSION U(12),DU(12) 08250000
+C 08260000
+C******************************************************************** 08270000
+C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-C) ***** 08280000
+C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** 08290000
+C******************************************************************** 08300000
+C 08310000
+C***** INPUT 08320000
+C***** XI,ET,Q,Z : STATION COORDINATES IN FAULT SYSTEM 08330000
+C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS 08340000
+C***** OUTPUT 08350000
+C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES 08360000
+C 08370000
+ COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D 08380000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 08390000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 08400000
+ DATA F0,F1,F2,F3,PI2/0.D0,1.D0,2.D0,3.D0,6.283185307179586D0/ 08410000
+C----- 08420000
+ C=D+Z 08430000
+ X53=(8.D0*R2+9.D0*R*XI+F3*XI2)*X11*X11*X11/R2 08440000
+ Y53=(8.D0*R2+9.D0*R*ET+F3*ET2)*Y11*Y11*Y11/R2 08450000
+ H=Q*CD-Z 08460000
+ Z32=SD/R3-H*Y32 08470000
+ Z53=F3*SD/R5-H*Y53 08480000
+ Y0=Y11-XI2*Y32 08490000
+ Z0=Z32-XI2*Z53 08500000
+ PPY=CD/R3+Q*Y32*SD 08510000
+ PPZ=SD/R3-Q*Y32*CD 08520000
+ QQ=Z*Y32+Z32+Z0 08530000
+ QQY=F3*C*D/R5-QQ*SD 08540000
+ QQZ=F3*C*Y/R5-QQ*CD+Q*Y32 08550000
+ XY=XI*Y11 08560000
+ QX=Q*X11 08570000
+ QY=Q*Y11 08580000
+ QR=F3*Q/R5 08590000
+ CQX=C*Q*X53 08600000
+ CDR=(C+D)/R3 08610000
+ YY0=Y/R3-Y0*CD 08620000
+C===== 08630000
+ DO 111 I=1,12 08640000
+ 111 U(I)=F0 08650000
+C====================================== 08660000
+C===== STRIKE-SLIP CONTRIBUTION ===== 08670000
+C====================================== 08680000
+ IF(DISL1.NE.F0) THEN 08690000
+ DU( 1)= ALP4*XY*CD -ALP5*XI*Q*Z32 08700000
+ DU( 2)= ALP4*(CD/R+F2*QY*SD) -ALP5*C*Q/R3 08710000
+ DU( 3)= ALP4*QY*CD -ALP5*(C*ET/R3-Z*Y11+XI2*Z32) 08720000
+ DU( 4)= ALP4*Y0*CD -ALP5*Q*Z0 08730000
+ DU( 5)=-ALP4*XI*(CD/R3+F2*Q*Y32*SD) +ALP5*C*XI*QR 08740000
+ DU( 6)=-ALP4*XI*Q*Y32*CD +ALP5*XI*(F3*C*ET/R5-QQ) 08750000
+ DU( 7)=-ALP4*XI*PPY*CD -ALP5*XI*QQY 08760000
+ DU( 8)= ALP4*F2*(D/R3-Y0*SD)*SD-Y/R3*CD 08770000
+ * -ALP5*(CDR*SD-ET/R3-C*Y*QR) 08780000
+ DU( 9)=-ALP4*Q/R3+YY0*SD +ALP5*(CDR*CD+C*D*QR-(Y0*CD+Q*Z0)*SD) 08790000
+ DU(10)= ALP4*XI*PPZ*CD -ALP5*XI*QQZ 08800000
+ DU(11)= ALP4*F2*(Y/R3-Y0*CD)*SD+D/R3*CD -ALP5*(CDR*CD+C*D*QR) 08810000
+ DU(12)= YY0*CD -ALP5*(CDR*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD) 08820000
+ DO 222 I=1,12 08830000
+ 222 U(I)=U(I)+DISL1/PI2*DU(I) 08840000
+ ENDIF 08850000
+C====================================== 08860000
+C===== DIP-SLIP CONTRIBUTION ===== 08870000
+C====================================== 08880000
+ IF(DISL2.NE.F0) THEN 08890000
+ DU( 1)= ALP4*CD/R -QY*SD -ALP5*C*Q/R3 08900000
+ DU( 2)= ALP4*Y*X11 -ALP5*C*ET*Q*X32 08910000
+ DU( 3)= -D*X11-XY*SD -ALP5*C*(X11-Q2*X32) 08920000
+ DU( 4)=-ALP4*XI/R3*CD +ALP5*C*XI*QR +XI*Q*Y32*SD 08930000
+ DU( 5)=-ALP4*Y/R3 +ALP5*C*ET*QR 08940000
+ DU( 6)= D/R3-Y0*SD +ALP5*C/R3*(F1-F3*Q2/R2) 08950000
+ DU( 7)=-ALP4*ET/R3+Y0*SDSD -ALP5*(CDR*SD-C*Y*QR) 08960000
+ DU( 8)= ALP4*(X11-Y*Y*X32) -ALP5*C*((D+F2*Q*CD)*X32-Y*ET*Q*X53) 08970000
+ DU( 9)= XI*PPY*SD+Y*D*X32 +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53) 08980000
+ DU(10)= -Q/R3+Y0*SDCD -ALP5*(CDR*CD+C*D*QR) 08990000
+ DU(11)= ALP4*Y*D*X32 -ALP5*C*((Y-F2*Q*SD)*X32+D*ET*Q*X53) 09000000
+ DU(12)=-XI*PPZ*SD+X11-D*D*X32-ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53) 09010000
+ DO 333 I=1,12 09020000
+ 333 U(I)=U(I)+DISL2/PI2*DU(I) 09030000
+ ENDIF 09040000
+C======================================== 09050000
+C===== TENSILE-FAULT CONTRIBUTION ===== 09060000
+C======================================== 09070000
+ IF(DISL3.NE.F0) THEN 09080000
+ DU( 1)=-ALP4*(SD/R+QY*CD) -ALP5*(Z*Y11-Q2*Z32) 09090000
+ DU( 2)= ALP4*F2*XY*SD+D*X11 -ALP5*C*(X11-Q2*X32) 09100000
+ DU( 3)= ALP4*(Y*X11+XY*CD) +ALP5*Q*(C*ET*X32+XI*Z32) 09110000
+ DU( 4)= ALP4*XI/R3*SD+XI*Q*Y32*CD+ALP5*XI*(F3*C*ET/R5-F2*Z32-Z0)09120000
+ DU( 5)= ALP4*F2*Y0*SD-D/R3 +ALP5*C/R3*(F1-F3*Q2/R2) 09130000
+ DU( 6)=-ALP4*YY0 -ALP5*(C*ET*QR-Q*Z0) 09140000
+ DU( 7)= ALP4*(Q/R3+Y0*SDCD) +ALP5*(Z/R3*CD+C*D*QR-Q*Z0*SD) 09150000
+ DU( 8)=-ALP4*F2*XI*PPY*SD-Y*D*X32 09160000
+ * +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53) 09170000
+ DU( 9)=-ALP4*(XI*PPY*CD-X11+Y*Y*X32) 09180000
+ * +ALP5*(C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)+XI*QQY) 09190000
+ DU(10)= -ET/R3+Y0*CDCD -ALP5*(Z/R3*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD) 09200000
+ DU(11)= ALP4*F2*XI*PPZ*SD-X11+D*D*X32 09210000
+ * -ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53) 09220000
+ DU(12)= ALP4*(XI*PPZ*CD+Y*D*X32) 09230000
+ * +ALP5*(C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)+XI*QQZ) 09240000
+ DO 444 I=1,12 09250000
+ 444 U(I)=U(I)+DISL3/PI2*DU(I) 09260000
+ ENDIF 09270000
+ RETURN 09280000
+ END 09290000
+
+
+
+
+
+
+
+ SUBROUTINE DCCON2(XI,ET,Q,SD,CD) 10170002
+ IMPLICIT REAL*8 (A-H,O-Z) 10180000
+C 10190000
+C********************************************************************** 10200000
+C***** CALCULATE STATION GEOMETRY CONSTANTS FOR FINITE SOURCE ***** 10210000
+C********************************************************************** 10220000
+C 10230000
+C***** INPUT 10240000
+C***** XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM 10250000
+C***** SD,CD : SIN, COS OF DIP-ANGLE 10260000
+C 10270000
+C### CAUTION ### IF XI,ET,Q ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZER010280000
+C 10290000
+ COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32, 10300000
+ * EY,EZ,FY,FZ,GY,GZ,HY,HZ 10310000
+ DATA F0,F1,F2,EPS/0.D0,1.D0,2.D0,1.D-6/ 10320000
+C----- 10330000
+ IF(DABS(XI).LT.EPS) XI=F0 10340000
+ IF(DABS(ET).LT.EPS) ET=F0 10350000
+ IF(DABS( Q).LT.EPS) Q=F0 10360000
+ XI2=XI*XI 10370000
+ ET2=ET*ET 10380000
+ Q2=Q*Q 10390000
+ R2=XI2+ET2+Q2 10400000
+ R =DSQRT(R2) 10410000
+ IF(R.EQ.F0) RETURN 10420000
+ R3=R *R2 10430000
+ R5=R3*R2 10440000
+ Y =ET*CD+Q*SD 10450000
+ D =ET*SD-Q*CD 10460000
+C----- 10470000
+ IF(Q.EQ.F0) THEN 10480000
+ TT=F0 10490000
+ ELSE 10500000
+ TT=DATAN(XI*ET/(Q*R)) 10510000
+ ENDIF 10520000
+C----- 10530000
+ IF(XI.LT.F0 .AND. Q.EQ.F0 .AND. ET.EQ.F0) THEN 10540002
+ ALX=-DLOG(R-XI) 10550000
+ X11=F0 10560000
+ X32=F0 10570000
+ ELSE 10580000
+ RXI=R+XI 10590002
+ ALX=DLOG(RXI) 10600000
+ X11=F1/(R*RXI) 10610000
+ X32=(R+RXI)*X11*X11/R 10620002
+ ENDIF 10630000
+C----- 10640000
+ IF(ET.LT.F0 .AND. Q.EQ.F0 .AND. XI.EQ.F0) THEN 10650002
+ ALE=-DLOG(R-ET) 10660000
+ Y11=F0 10670000
+ Y32=F0 10680000
+ ELSE 10690000
+ RET=R+ET 10700002
+ ALE=DLOG(RET) 10710000
+ Y11=F1/(R*RET) 10720000
+ Y32=(R+RET)*Y11*Y11/R 10730002
+ ENDIF 10740000
+C----- 10750000
+ EY=SD/R-Y*Q/R3 10760000
+ EZ=CD/R+D*Q/R3 10770000
+ FY=D/R3+XI2*Y32*SD 10780000
+ FZ=Y/R3+XI2*Y32*CD 10790000
+ GY=F2*X11*SD-Y*Q*X32 10800000
+ GZ=F2*X11*CD+D*Q*X32 10810000
+ HY=D*Q*X32+XI*Q*Y32*SD 10820000
+ HZ=Y*Q*X32+XI*Q*Y32*CD 10830000
+ RETURN 10840000
+ END 10850000
Added: cs/cigma/trunk/sandbox/okada/disloc3d/h5array.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/h5array.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/h5array.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,95 @@
+#include <stdlib.h>
+#include "hdf5.h"
+#include "h5array.h"
+
+int h5array_read(hid_t loc_id, const char *name,
+ double **data, int *m, int *n)
+{
+ hid_t dataset_id;
+ hid_t dataspace_id;
+ hsize_t dims[2];
+ herr_t status;
+
+ int rank;
+
+
+ dataset_id = H5Dopen(loc_id, name);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+
+ dataspace_id = H5Dget_space(dataset_id);
+ rank = H5Sget_simple_extent_ndims(dataspace_id);
+ if (rank != 2)
+ {
+ H5Dclose(dataset_id);
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+ status = H5Sclose(dataspace_id);
+
+ *m = dims[0];
+ *n = dims[1];
+
+ *data = (double *)malloc((dims[0] * dims[1]) * sizeof(double));
+ if (*data == NULL)
+ {
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, *data);
+ if (status < 0)
+ {
+ H5Dclose(dataset_id);
+ return -4;
+ }
+
+ status = H5Dclose(dataset_id);
+
+ return 0;
+}
+
+int h5array_write(hid_t loc_id, const char *name,
+ double *data, int m, int n)
+{
+ hsize_t dims[2];
+ hid_t dataspace_id;
+ hid_t dataset_id;
+ herr_t status;
+
+ dims[0] = m;
+ dims[1] = n;
+ dataspace_id = H5Screate_simple(2, dims, NULL);
+ if (dataspace_id < 0)
+ {
+ return -1;
+ }
+
+ dataset_id = H5Dcreate(loc_id, name, H5T_NATIVE_DOUBLE, dataspace_id,
+ H5P_DEFAULT);
+ if (dataset_id < 0)
+ {
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+ status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, data);
+ if (status < 0)
+ {
+ H5Sclose(dataspace_id);
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ H5Sclose(dataspace_id);
+ H5Dclose(dataset_id);
+
+ return 0;
+}
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/h5array.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/h5array.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/h5array.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,11 @@
+#ifndef __H5ARRAY_H__
+#define __H5ARRAY_H__
+#include "hdf5.h"
+
+int h5array_read(hid_t loc_id, const char *name,
+ double **data, int *m, int *n);
+
+int h5array_write(hid_t loc_id, const char *name,
+ double *data, int m, int n);
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/disloc3d/h5points.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/h5points.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/h5points.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,77 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <hdf5.h>
+#include "txtarray.h"
+#include "h5array.h"
+
+int main(int argc, char *argv[])
+{
+ char infile[100];
+ char outfile[100];
+
+ hid_t file_id;
+ herr_t status;
+
+ int nstations;
+ double *stations;
+
+
+
+ if (argc == 1)
+ {
+ fprintf(stderr, "Usage: %s points.in points.out\n", argv[0]);
+ return -1;
+ }
+
+ if (argc > 1)
+ {
+ snprintf(infile, 100, argv[1]);
+ }
+ else
+ {
+ snprintf(infile, 100, "points.in");
+ }
+
+ if (argc > 2)
+ {
+ snprintf(outfile, 100, argv[2]);
+ }
+ else
+ {
+ snprintf(outfile, 100, "points.out");
+ }
+
+
+ txtarray_read_stations(infile, &stations, &nstations); // allocates stations
+
+
+ file_id = H5Fcreate(outfile, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write to file %s\n", outfile);
+ return -2;
+ }
+
+
+ status = h5array_write(file_id, "/stations", stations, nstations, 3);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset /stations\n");
+ return -3;
+ }
+
+
+ status = H5Fclose(file_id);
+
+ free(stations);
+
+
+ printf("Processed %d stations into dataset /stations in %s\n", nstations, outfile);
+
+
+ return EXIT_SUCCESS;
+}
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/mexfunction.f
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/mexfunction.f 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/mexfunction.f 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,260 @@
+ SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
+
+C USE MSFLIB
+
+C These are arrays of pointers. Size needs to match size of pointers.
+C Use INTEGER*8 on 64-bit machines, INTEGER*4 on 32-bit machines.
+ INTEGER*4 PLHS(*), PRHS(*)
+
+ INTEGER NLHS, NRHS
+
+C These are functions that return pointers. Size needs to match size
+C of pointers. Use INTEGER*8 on 64-bit machines, INTEGER*4 on
+C 32-bit machines.
+ INTEGER*4 MXCREATEFULL, MXGETPR
+
+ INTEGER MXGETM, MXGETN
+C
+C KEEP THE ABOVE SUBROUTINE, ARGUMENT, AND FUNCTION DECLARATIONS FOR USE
+C IN ALL YOUR FORTRAN MEX FILES.
+C---------------------------------------------------------------------
+C
+
+C These are actually pointers. Size needs to match size of pointers.
+C Use INTEGER*8 on 64-bit machines, INTEGER*4 on 32-bit machines.
+ INTEGER*4 pModel, pStat, pMu, pNU
+ INTEGER*4 UOUT, DOUT, SOUT, FLAGOUT, FLAGOUT2
+
+ INTEGER M, N, NSTAT, NMOD, I, J
+ REAL*8 MODEL(10), STAT(3), MU, LAMBDA, NU
+ REAL*8 ALPHA, STRIKE, DIP, DEG2RAD, SD, CD, SS, CS, DEPTH
+ REAL*8 X,Y,Z,DISL1,DISL2,DISL3,AL1,AL2,AW1,AW2
+ REAL*8 U(3),D(9),S(6),UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ
+ REAL*8 UXT,UYT,UZT,UXXT,UYXT,UZXT,UXYT,UYYT,UZYT,UXZT,UYZT,UZZT
+ REAL*8 FLAG
+ REAL*8 A, B, C, AAA, BBB, CCC, DDD, EEE, FFF, GGG, HHH, III
+ REAL*8 FLAG2(5000,5000)
+
+C For Windows only!
+C This resets the floating point exception to allow divide by zero,
+C overflow and invalid numbers.
+C
+C INTEGER(2) CONTROL
+C CALL GETCONTROLFPQQ(CONTROL)
+C CONTROL = CONTROL .OR. FPCW$ZERODIVIDE
+C CONTROL = CONTROL .OR. FPCW$INVALID
+C CONTROL = CONTROL .OR. FPCW$OVERFLOW
+C CALL SETCONTROLFPQQ(CONTROL)
+
+C
+C CHECK FOR PROPER NUMBER OF ARGUMENTS
+C
+ IF (NRHS .NE. 4) THEN
+ CALL MEXERRMSGTXT('Usage: [U,D,S,flag]=disloc3d(m,x,mu,nu)')
+ ENDIF
+C
+C CHECK THE ARGUMENTS
+C
+ M = MXGETM(PRHS(1))
+ NMOD = MXGETN(PRHS(1))
+ IF (M .NE. 10) THEN
+ CALL MEXERRMSGTXT('m must be 10x1 model vector')
+ ENDIF
+
+ M = MXGETM(PRHS(2))
+ NSTAT = MXGETN(PRHS(2))
+ IF (M .NE. 3) THEN
+ CALL MEXERRMSGTXT('x must be 3xn.')
+ ENDIF
+
+ M = MXGETM(PRHS(3))
+ N = MXGETN(PRHS(3))
+ IF ((M .NE. 1) .OR. (N .NE. 1)) THEN
+ CALL MEXERRMSGTXT('mu must be a scalar.')
+ ENDIF
+
+ M = MXGETM(PRHS(4))
+ N = MXGETN(PRHS(4))
+ IF ((M .NE. 1) .OR. (N .NE. 1)) THEN
+ CALL MEXERRMSGTXT('nu must be a scalar.')
+ ENDIF
+
+C
+C CREATE A MATRIX FOR RETURN ARGUMENT
+C
+ PLHS(1) = MXCREATEFULL(3,NSTAT,0)
+ PLHS(2) = MXCREATEFULL(9,NSTAT,0)
+ PLHS(3) = MXCREATEFULL(6,NSTAT,0)
+ PLHS(4) = MXCREATEFULL(NSTAT,1,0)
+ PLHS(5) = MXCREATEFULL(NMOD,NSTAT,0)
+C
+C ASSIGN POINTERS TO THE VARIOUS PARAMETERS
+C
+ UOUT = MXGETPR(PLHS(1))
+ DOUT = MXGETPR(PLHS(2))
+ SOUT = MXGETPR(PLHS(3))
+ FLAGOUT = MXGETPR(PLHS(4))
+ FLAGOUT2= MXGETPR(PLHS(5))
+
+C
+ pMODEL = MXGETPR(PRHS(1))
+ pSTAT = MXGETPR(PRHS(2))
+ pMU = MXGETPR(PRHS(3))
+ pNU = MXGETPR(PRHS(4))
+C
+C COPY RIGHT HAND ARGUMENTS TO LOCAL ARRAYS OR VARIABLES
+C
+ CALL MXCOPYPTRTOREAL8(pMU, MU, 1)
+ CALL MXCOPYPTRTOREAL8(pNU, NU, 1)
+ LAMBDA = 2*MU*NU/(1-2*NU)
+
+C LOOP OVER STATIONS
+
+ DO 111 I=1,NSTAT
+
+ CALL MXCOPYPTRTOREAL8(pSTAT+(I-1)*24, STAT, 3)
+
+C GENERATE WARNINGS FOR POSITIVE DEPTHS
+
+ IF (STAT(3) .GT. 0) THEN
+ CALL MEXEVALSTRING("warning('Positive depth given.')")
+ ELSE
+
+C LOOP OVER MODELS
+
+ DO 222 J=1,NMOD
+
+ CALL MXCOPYPTRTOREAL8(pMODEL+(J-1)*80, MODEL, 10)
+
+ DEG2RAD=3.14159265358979/180
+ ALPHA = (LAMBDA+MU)/(LAMBDA+2*MU)
+ STRIKE = (MODEL(5)-90)*DEG2RAD
+ CS = DCOS(STRIKE)
+ SS = DSIN(STRIKE)
+ DIP = MODEL(4)
+ CD = DCOS(DIP*DEG2RAD)
+ SD = DSIN(DIP*DEG2RAD)
+ DISL1=MODEL(8)
+ DISL2=MODEL(9)
+ DISL3=MODEL(10)
+ DEPTH=MODEL(3)-0.5*MODEL(2)*SD
+ AL1=MODEL(1)/2
+ AL2=AL1
+ AW1=MODEL(2)/2
+ AW2=AW1
+ X=CS*(-MODEL(6) + STAT(1)) - SS*(-MODEL(7) + STAT(2))
+ Y=-0.5*CD*MODEL(2) + SS*(-MODEL(6) + STAT(1)) + CS*(-MODEL(7) +
+ - STAT(2))
+ Z=STAT(3)
+
+C GENERATE WARNINGS FOR UNPHYSICAL MODELS
+
+ IF ((MODEL(3)-SD*MODEL(2) .LT. 0) .OR.
+ - (MODEL(1) .LE. 0) .OR.
+ - (MODEL(2) .LE. 0) .OR.
+ - (MODEL(3) .LT. 0)) THEN
+ CALL MEXEVALSTRING("warning('Unphysical model.')")
+ ELSE
+
+ CALL DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
+ * AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
+ * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,FLAG)
+
+C fill flag2 - the rows contain flags for each station, one row per model/fault
+
+
+ FLAG2(J,I) = FLAG + 5 + 5*FLAG
+
+
+
+C ROTATE THEN ADD
+ A = CS*UX + SS*UY
+ B = -SS*UX + CS*UY
+ C = UZ
+ AAA = CS**2*UXX + CS*SS*(UXY + UYX) + SS**2*UYY
+ BBB = CS**2*UXY - SS**2*UYX + CS*SS*(-UXX + UYY)
+ CCC = CS*UXZ + SS*UYZ
+ DDD = -(SS*(CS*UXX + SS*UXY)) + CS*(CS*UYX + SS*UYY)
+ EEE = SS**2*UXX - CS*SS*(UXY + UYX) + CS**2*UYY
+ FFF = -(SS*UXZ) + CS*UYZ
+ GGG = CS*UZX + SS*UZY
+ HHH = -(SS*UZX) + CS*UZY
+ III = UZZ
+
+
+
+
+
+ UXT=UXT + A
+ UYT=UYT + B
+ UZT=UZT + C
+
+ UXXT=UXXT + AAA
+ UXYT=UXYT + BBB
+ UXZT=UXZT + CCC
+ UYXT=UYXT + DDD
+ UYYT=UYYT + EEE
+ UYZT=UYZT + FFF
+ UZXT=UZXT + GGG
+ UZYT=UZYT + HHH
+ UZZT=UZZT + III
+ ENDIF
+
+ 222 CONTINUE
+
+C ASSIGN OUTPUTS
+
+ U(1)=UXT
+ U(2)=UYT
+ U(3)=UZT
+
+ D(1)=UXXT
+ D(2)=UXYT
+ D(3)=UXZT
+ D(4)=UYXT
+ D(5)=UYYT
+ D(6)=UYZT
+ D(7)=UZXT
+ D(8)=UZYT
+ D(9)=UZZT
+
+
+
+C CALCULATE STRESSES
+
+ THETA=D(1)+D(5)+D(9)
+
+ S(1)=LAMBDA*THETA+2*MU*D(1)
+ S(2)=MU*(D(2)+D(4))
+ S(3)=MU*(D(3)+D(7))
+ S(4)=LAMBDA*THETA+2*MU*D(5)
+ S(5)=MU*(D(6)+D(8))
+ S(6)=LAMBDA*THETA+2*MU*D(9)
+
+C COPY TO MATLAB
+
+ CALL MXCOPYREAL8TOPTR(U, UOUT+(I-1)*24, 3)
+ CALL MXCOPYREAL8TOPTR(D, DOUT+(I-1)*72, 9)
+ CALL MXCOPYREAL8TOPTR(S, SOUT+(I-1)*48, 6)
+ CALL MXCOPYREAL8TOPTR(FLAG, FLAGOUT+(I-1)*8, 1)
+
+ UXT=0
+ UYT=0
+ UZT=0
+
+ UXXT=0
+ UXYT=0
+ UXZT=0
+ UYXT=0
+ UYYT=0
+ UYZT=0
+ UZXT=0
+ UZYT=0
+ UZZT=0
+ ENDIF
+ 111 CONTINUE
+
+ CALL MXCOPYREAL8TOPTR(FLAG2, FLAGOUT2, NSTAT*NMOD)
+C
+ RETURN
+ END
Added: cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,501 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "cervelli.h"
+#include "txtarray.h"
+#include "h5array.h"
+
+
+static void linspace(double **x, double a, double b, int n)
+{
+ int i;
+ double dx = (b - a)/n;
+ double *y;
+
+ y = (double *)malloc(n * sizeof(double));
+
+ for (i = 0; i < n; i++)
+ {
+ y[i] = a + i*dx;
+ }
+
+ *x = y;
+}
+
+
+void strike_slip_ng(double *fault_info, int nfaults)
+{
+ int i;
+ double *n;
+ double taperd;
+
+ double *fi;
+ double fx1, fy1;
+ double fx2, fy2;
+ double dip, D, Pr, bd;
+ double ss, ds, ts;
+
+ //
+ // define number of faults to taper with
+ //
+ taperd = 10.0; // (16000.0 - 12000.0)/400
+ linspace(&n, 12000.0, 16000.0, nfaults);
+
+ for (i = 0; i < nfaults; i++)
+ {
+ fi = &fault_info[11*i];
+
+ fx1 = 12000.0;
+ fy1 = -n[i];
+
+ fx2 = 12000.0;
+ fy2 = n[i];
+
+ dip = 90.0;
+ D = n[i];
+ Pr = 0.25;
+ bd = 0;
+ ss = -1.0/nfaults;
+ ds = 0;
+ ts = 0;
+
+ fi[0] = fx1;
+ fi[1] = fy1;
+ fi[2] = fx2;
+ fi[3] = fy2;
+ fi[4] = dip;
+ fi[5] = D;
+ fi[6] = Pr;
+ fi[7] = bd;
+ fi[8] = ss;
+ fi[9] = ds;
+ fi[10] = ts;
+ }
+
+ free(n);
+}
+
+void reverse_slip_ng(double *fault_info, int nfaults)
+{
+ int i;
+ double *n;
+ double taperd;
+
+ double *fi;
+ double fx1, fy1;
+ double fx2, fy2;
+ double dip, D, Pr, bd;
+ double ss, ds, ts;
+
+ //
+ // define number of faults to taper with
+ //
+ taperd = 10.0; // (16000.0 - 12000.0)/400
+ linspace(&n, 12000.0, 16000.0, nfaults);
+
+ for (i = 0; i < nfaults; i++)
+ {
+ fi = &fault_info[11*i];
+
+ fx1 = 4000.0;
+ fy1 = -n[i];
+
+ fx2 = 4000.0;
+ fy2 = n[i];
+
+ dip = 90.0;
+ D = n[i];
+ Pr = 0.25;
+ bd = 0;
+ ss = 0;
+ ds = 1.0/nfaults;
+ ts = 0;
+
+ fi[0] = fx1;
+ fi[1] = fy1;
+ fi[2] = fx2;
+ fi[3] = fy2;
+ fi[4] = dip;
+ fi[5] = D;
+ fi[6] = Pr;
+ fi[7] = bd;
+ fi[8] = ss;
+ fi[9] = ds;
+ fi[10] = ts;
+ }
+
+ free(n);
+}
+
+
+/*
+ * Note: Be careful with flag2. For a model with
+ * 400 subfaults and 918980 stations, this function would try
+ * to allocate
+ *
+ * (918980*(3+3+9+6+1+400) + 400*(10+11)) * 8 bytes = 2958.8 MiB
+ *
+ * disabling flag2 gives us
+ *
+ * (918980*(3+3+9+6+1) + 400*(10+11)) * 8 bytes = 154.3 MiB
+ *
+ */
+void okada_alloc(int nmod, int nstat,
+ double **fault_info,
+ double **models,
+ double **stations,
+ double **displacement,
+ double **derivatives,
+ double **stress,
+ double **flag,
+ double **flag2)
+{
+ if (fault_info != NULL)
+ *fault_info = (double *)malloc((nmod * 11) * sizeof(double));
+
+ if (models != NULL)
+ *models = (double *)malloc((nmod * 10) * sizeof(double));
+
+ if (stations != NULL)
+ *stations = (double *)malloc((nstat * 3) * sizeof(double));
+
+ if (displacement != NULL)
+ *displacement = (double *)malloc((nstat * 3) * sizeof(double));
+
+ if (derivatives != NULL)
+ *derivatives = (double *)malloc((nstat * 9) * sizeof(double));
+
+ if (stress != NULL)
+ *stress = (double *)malloc((nstat * 6) * sizeof(double));
+
+ if (flag != NULL)
+ *flag = (double *)malloc((nstat * 1) * sizeof(double));
+
+ if (flag2 != NULL)
+ *flag2 = (double *)malloc((nstat * nmod) * sizeof(double));
+}
+
+void okada_free(double *fault_info, double *models, double *stations,
+ double *displacement, double *derivatives, double *stress,
+ double *flag, double *flag2)
+{
+ if (fault_info != NULL) free(fault_info);
+ if (models != NULL) free(models);
+ if (stations != NULL) free(stations);
+ if (displacement != NULL) free(displacement);
+ if (derivatives != NULL) free(derivatives);
+ if (stress != NULL) free(stress);
+ if (flag != NULL) free(flag);
+ if (flag2 != NULL) free(flag2);
+}
+
+
+
+int txt_main(char *infile, char *outfile)
+{
+ int nsubfaults;
+ double *subfaults;
+ double *models;
+
+ int nstations;
+ double *stations;
+ double *displacement;
+ double *derivatives;
+ double *stress;
+
+ double *flag;
+ //double *flag2;
+
+
+ // elastic parameters;
+ double mu = 0.25;
+ double nu = 0.25;
+
+
+
+ /* determine number of subfaults */
+ nsubfaults = 400;
+
+
+ /* read station coordinates */
+ txtarray_read_stations(infile, &stations, &nstations);
+
+ printf("Found %d stations\n", nstations);
+
+ /*
+ * Allocate memory. Skip stations array, since it has
+ * already been allocated above. Also, do not allocate
+ * flag2 array, since it is typically too large.
+ */
+ okada_alloc(nsubfaults, nstations,
+ &subfaults, &models,
+ NULL, &displacement, &derivatives, &stress,
+ &flag, NULL);
+
+
+ /*
+ * Pick one of the following models
+ */
+ strike_slip_ng(subfaults, nsubfaults);
+ //reverse_slip_ng(subfaults, nsubfaults);
+
+
+ /*
+ * Perform the actual calculation
+ */
+ calc_disp_cervelli(mu, nu,
+ models, subfaults, nsubfaults,
+ stations, nstations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+
+ /*
+ * Write all data to a text file
+ */
+ txtarray_write_all(outfile,
+ "sx sy sz "
+ "ux uy uz "
+ "uxx uxy uxz uyx uyy uyz uzx uzy uzz "
+ "sxx syy szz sxy sxz syz",
+ nstations,
+ stations,
+ displacement,
+ derivatives,
+ stress);
+
+
+ /* free the allocated memory */
+ okada_free(subfaults, models, stations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+ return 0;
+}
+
+
+
+int split_filepath(const char *filepath, char *filename, char *path, int maxlen)
+{
+ char *x, *y;
+ int i;
+
+ y = strchr(filepath, ':');
+
+ if (y == NULL)
+ {
+ strncpy(filename, filepath, maxlen);
+ path[0] = '\0';
+ return 1;
+ }
+
+ for (i = 0, x = (char *)filepath; (x != y) && (i < maxlen); i++, x++)
+ {
+ filename[i] = *x;
+ }
+ filename[i] = '\0';
+
+ for(i = 0, x = y+1; (*x != '\0') && (i < maxlen); i++, x++)
+ {
+ path[i] = *x;
+ }
+ path[i] = '\0';
+
+ return 0;
+}
+
+int h5_main(char *filepath1, char *filepath2)
+{
+ hid_t file_id;
+ //hid_t group_id;
+ herr_t status;
+
+ char infile[150];
+ char outfile[150];
+ char stations_path[150];
+ char output_group[150];
+
+ int nsubfaults;
+ double *subfaults;
+ double *models;
+
+ int nstations;
+ double *stations;
+ double *displacement;
+ double *derivatives;
+ double *stress;
+
+ double *flag;
+ //double flag2;
+
+ int components;
+
+ // elastic parameters
+ double mu = 0.25;
+ double nu = 0.25;
+
+
+ /* determine number of subfaults */
+ nsubfaults = 400;
+
+
+ /*
+ * Process filepaths
+ */
+
+ split_filepath(filepath1, infile, stations_path, 150);
+ if ('\0' == stations_path[0])
+ {
+ fprintf(stderr, "Error: Could not split file and path from \"%s\"\n", filepath1);
+ return -1;
+ }
+
+ split_filepath(filepath2, outfile, output_group, 150);
+ if (output_group[0] != '\0')
+ {
+ fprintf(stderr, "Error: Writing to an output group currently not supported\n");
+ return -2;
+ }
+
+
+ /*
+ * Read station coordinates
+ */
+
+
+ file_id = H5Fopen(infile, H5F_ACC_RDONLY, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ fprintf(stderr, "Error: Could not open file %s\n", infile);
+ return -3;
+ }
+
+ status = h5array_read(file_id, stations_path, &stations, &nstations, &components);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ return -4;
+ }
+ if (components != 3)
+ {
+ return -5;
+ }
+
+ status = H5Fclose(file_id);
+
+
+ /*
+ * Allocate memory. Skip stations array, since it has
+ * already been allocated above. Also, do not allocate
+ * flag2 array, since it is typically too large.
+ */
+ okada_alloc(nsubfaults, nstations,
+ &subfaults, &models,
+ NULL, &displacement, &derivatives, &stress,
+ &flag, NULL);
+
+ /*
+ * Pick one of the following models
+ */
+ strike_slip_ng(subfaults, nsubfaults);
+ //reverse_slip_ng(subfaults, nsubfaults);
+
+
+ /*
+ * Perform the actual calculation
+ */
+ calc_disp_cervelli(mu, nu,
+ models, subfaults, nsubfaults,
+ stations, nstations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+
+ /*
+ * Write all data to an HDF5 file (to root group for now)
+ *
+ * TODO: Specify the parent group of these datasets.
+ * Essentially, I need a function that creates
+ * all the necessary groups in a path.
+ */
+
+ file_id = H5Fcreate(outfile, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write to file %s\n", outfile);
+ return -6;
+ }
+
+ status = h5array_write(file_id, "/stations", stations, nstations, 3);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/stations");
+ }
+
+ status = h5array_write(file_id, "/displacement", displacement, nstations, 3);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/displacement");
+ }
+
+ status = h5array_write(file_id, "/derivatives", derivatives, nstations, 9);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/derivatives");
+ }
+
+ status = h5array_write(file_id, "/stress", stress, nstations, 6);
+ if (status < 0)
+ {
+ H5Fclose(file_id);
+ fprintf(stderr, "Error: Could not write dataset %s\n", "/stress");
+ }
+
+ status = H5Fclose(file_id);
+
+
+ /* free the allocated memory */
+ okada_free(subfaults, models, stations,
+ displacement, derivatives, stress,
+ flag, NULL);
+
+ return 0;
+}
+
+
+int main(int argc, char *argv[])
+{
+ char arg1[300];
+ char arg2[300];
+
+ if (argc < 4)
+ {
+ fprintf(stderr, "Usage: %s [hdf5 | txt] points.in points.out\n", argv[0]);
+ return -1;
+ }
+
+ snprintf(arg1, 300, argv[2]); // points.in or points.h5:/okada/stations
+ snprintf(arg2, 300, argv[3]); // points.out or data.h5
+
+ if (strcmp(argv[1], "hdf5") == 0)
+ {
+ return h5_main(arg1, arg2);
+ }
+ else if (strcmp(argv[1], "txt") == 0)
+ {
+ return txt_main(arg1, arg2);
+ }
+ else
+ {
+ fprintf(stderr, "Error: Could not understand option \"%s\"\n", argv[3]);
+ }
+
+
+ return EXIT_FAILURE;
+}
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.m
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.m 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/okadasoln_ng.m 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,412 @@
+function okadasoln_ng(varargin)
+%
+% Matlab script to generate the analytic boundary conditions for
+% the Reverse-Slip and Strike-Slip (no gravity) benchmarks.
+%
+% This script is essentially bm5.m by Noah Fay (the version for the
+% old benchmark 5), all comments below this are his (slip taper was
+% modified so that zero slip at the edge of the fault plane in the
+% FEMs is also zero slip in this code)...
+%
+% ARGUMENTS:
+% #1 - Subfaults information [REQUIRED]
+% #2 - Name of input file [OPTIONAL]
+
+%% Inputs:
+%% points.in file which must exist in the local directory
+%% or somewhere in path. this file contains
+%% all the xyz triples at which displacement
+%% are calculated. No headerlines.
+%%
+%% Outputs:
+%% disp.out file with 6 columns:
+%% xcoord ycoord zcoord ux uy uz
+%%
+%% Notes: --See bm5_analytic_README.txt for further discussion
+%% --All units (.in and .out files) are SI (use meters)
+%% --Some things are hardcoded:
+%% -Taper depth from 12 to 16km
+%% -Taper step, taperd = 10m
+%% -Dip angle, dip = 45
+%% -xcoord of fault, fx1 = fx2 = 4km
+%% -Poisson's ratio for elastic solid, Pr = 0.25
+%% -Elastic constants lamda = mu = 30e9 Pa
+%% Credits: --The mex function (disloc3d_mod2) is a modified (and
+%% corrected) version of P.Cervelli's disloc3d.f
+%% --function fault_params_to_okada_form.m is courtesy of
+%% Brendan Meade, MIT
+%%
+%%
+%% Noah Fay, 2003
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Get the stations (anywhwere) to calculate displacement,
+%% strain, strain energy, and total strain energy (sum for those pts)
+
+if length(varargin) < 1
+ error('usage: okadasoln(subfaults [, filenameIn])')
+ return
+else
+ subfaults = varargin{1};
+ if length(varargin)==2
+ filein = varargin{2};
+ fileout = sprintf('%s.disp.out',filein);
+ else
+ filein = 'points.in';
+ fileout = 'disp.out';
+ end
+end
+
+
+[sx,sy,sz] = textread(filein,'%f %f %f');
+
+
+%% Read in node coords for calculating displacement at
+%% the sides of the box for the analytic bc
+%% [node, sx, sy, sz] = textread('node_coords.txt','%f %f %f %f');
+%% Read in the face of interest
+%% face = 'back';
+%% fid = fopen([face,'.nodes']);
+%% i = fscanf(fid, '%f');
+%% Get the stations that belong to that face
+%% sx = sx(i);
+%% sy = sy(i);
+%% sz = sz(i);
+
+
+%% Calculate displacements, derivatives
+[dispx,dispy,dispz,D,S] = calc_disp_cervelli(sx,sy,sz,subfaults);
+
+% $$$ %% Calculate strains strij:
+% $$$ strxx = D(1,:);
+% $$$ strxy = 0.5 * (D(2,:) + D(4,:));
+% $$$ strxz = 0.5 * (D(3,:) + D(7,:));
+% $$$
+% $$$ stryx = strxy;
+% $$$ stryy = D(5,:);
+% $$$ stryz = 0.5 * (D(6,:) + D(8,:));
+% $$$
+% $$$ strzx = strxz;
+% $$$ strzy = stryz;
+% $$$ strzz = D(9,:);
+% $$$
+% $$$ %% Calculate stress tauij
+% $$$ lamda = 30e9;
+% $$$ mu = lamda;
+% $$$ strkk = sum([strxx; stryy; strzz]);
+% $$$
+% $$$ tauxx = lamda*strkk + 2*mu*strxx;
+% $$$ tauxy = 2*mu*strxy;
+% $$$ tauxz = 2*mu*strxz;
+% $$$
+% $$$ tauyx = 2*mu*stryx;
+% $$$ tauyy = lamda*strkk + 2*mu*stryy;
+% $$$ tauyz = 2*mu*stryz;
+% $$$
+% $$$ tauzx = 2*mu*strzx;
+% $$$ tauzy = 2*mu*strzy;
+% $$$ tauzz = lamda*strkk + 2*mu*strzz;
+% $$$
+% $$$
+% $$$ %% Calculate strain energy
+% $$$ Exx = strxx.*tauxx;
+% $$$ Exy = strxy.*tauxy;
+% $$$ Exz = strxz.*tauxz;
+% $$$
+% $$$ Eyx = stryx.*tauyx;
+% $$$ Eyy = stryy.*tauyy;
+% $$$ Eyz = stryz.*tauyz;
+% $$$
+% $$$ Ezx = strzx.*tauzx;
+% $$$ Ezy = strzy.*tauzy;
+% $$$ Ezz = strzz.*tauzz;
+% $$$
+% $$$ Etot = 0.5 * sum([Exx; Exy; Exz; Eyx; Eyy; Eyz; Ezx; Ezy; Ezz]);
+
+
+%% Print a file with displacements
+fid = fopen(fileout,'w');
+for i = 1:1:length(sx)
+ fprintf(fid,'%0.12e %0.12e %0.12e %0.12e %0.12e %0.12e\n',...
+ sx(i),sy(i),sz(i),dispx(i),dispy(i),dispz(i));
+% $$$ fprintf(fid, '%e\t', sx(i));
+% $$$ fprintf(fid, '%e\t', sy(i));
+% $$$ fprintf(fid, '%e\t', sz(i));
+% $$$
+% $$$ fprintf(fid, '%0.12e\t', dispx(i));
+% $$$ fprintf(fid, '%0.12e\t', dispy(i));
+% $$$ fprintf(fid, '%0.12e\n', dispz(i));
+end
+fclose(fid);
+
+% $$$ %% Print a file with strains
+% $$$ fid = fopen('strain.out','w');
+% $$$ for i = 1:1:length(sx)
+% $$$ fprintf(fid, '%e\t', sx(i));
+% $$$ fprintf(fid, '%e\t', sy(i));
+% $$$ fprintf(fid, '%e\t', sz(i));
+% $$$
+% $$$ fprintf(fid, '%e\t', strxx(i));
+% $$$ fprintf(fid, '%e\t', strxy(i));
+% $$$ fprintf(fid, '%e\t', strxz(i));
+% $$$ fprintf(fid, '%e\t', stryx(i));
+% $$$ fprintf(fid, '%e\t', stryy(i));
+% $$$ fprintf(fid, '%e\t', stryz(i));
+% $$$ fprintf(fid, '%e\t', strzx(i));
+% $$$ fprintf(fid, '%e\t', strzy(i));
+% $$$ fprintf(fid, '%e\n', strzz(i));
+% $$$ end
+% $$$ fclose(fid);
+% $$$
+% $$$ %% Print a file with strain energy (9) and total
+% $$$ fid = fopen('E.out','w');
+% $$$ for i = 1:1:length(sx)
+% $$$ fprintf(fid, '%e\t', sx(i));
+% $$$ fprintf(fid, '%e\t', sy(i));
+% $$$ fprintf(fid, '%e\t', sz(i));
+% $$$
+% $$$ %%fprintf(fid, '%e\t', Exx(i));
+% $$$ %%fprintf(fid, '%e\t', Exy(i));
+% $$$ %%fprintf(fid, '%e\t', Exz(i));
+% $$$ %%fprintf(fid, '%e\t', Eyx(i));
+% $$$ %%fprintf(fid, '%e\t', Eyy(i));
+% $$$ %%fprintf(fid, '%e\t', Eyz(i));
+% $$$ %%fprintf(fid, '%e\t', Ezx(i));
+% $$$ %%fprintf(fid, '%e\t', Ezy(i));
+% $$$ %%fprintf(fid, '%e\t', Ezz(i));
+% $$$ fprintf(fid, '%e\n', Etot(i));
+% $$$ end
+% $$$ fclose(fid);
+% $$$
+% $$$ %% Print a file with sum of all Etot (total strain energy for all pts)
+% $$$ fid = fopen('Esum.out','w');
+% $$$ fprintf(fid, '%e', sum(Etot));
+% $$$ fclose(fid);
+
+%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-
+function [dispx, dispy, dispz,D,S] = ...
+ calc_disp_cervelli(sx, sy, sz, fault_info)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% function calc_disp_cervelli
+%% Calculate dipslacements via Cervelli's Okada (1992) code
+%%
+%% Inputs: sx station x coord
+%% sy station y coord
+%% fault_info matrix with fault stuff
+%%
+%% Outputs: dispx disp in x direction
+%% dispy disp in y direction
+%% dispz disp in z direciton
+%% D matrix with derivatives
+%% S Stresses (do not use)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Introduce a dummy station at the beginning
+%% of the list so the cervelli fortran code
+%% doesn't freak out...remove it later
+ sx = [-rand; sx];
+ sy = [-rand; sy];
+ sz = [-rand; sz];
+
+%% Set station depth and elastic parameters
+ mu = 0.25; %% poissons ratio
+ nu = 0.25;
+
+%% get the "model matrix" with fault patch
+%% info in the form for cervelli calcs
+ M = getM(fault_info,mu);
+
+%% Preallocate vectors for better memory management
+nstations = length(sx);
+dispx = zeros(1, nstations);
+dispy = zeros(1, nstations);
+dispz = zeros(1, nstations);
+D = zeros(9, nstations);
+S = zeros(6, nstations);
+
+%% Loop over a subset of stations b/c Cervelli
+%% code doesn't seem to like it if given more
+%% than ~500 stations
+for subset = 1:1:ceil(nstations/500)
+ a = (subset-1)*500 + 1;
+ b = subset*500;
+ if b > nstations
+ b = nstations;
+ end
+ X = [sx(a:b)'; sy(a:b)'; sz(a:b)'];
+
+ [Usub,Dsub,Ssub,flag,flag2] = disloc3d_mod2(M,X,mu,nu);
+ dispx(a:b) = Usub(1,:)';
+ dispy(a:b) = Usub(2,:)';
+ dispz(a:b) = Usub(3,:)';
+ D(:,a:b) = Dsub;
+ S(:,a:b) = Ssub;
+end %%subset
+clear X Usub Dsub Ssub flag flag2
+
+%% Remove the dummy displacements at the beginning
+%% of the displacement lists
+dispx(1) = [];
+dispy(1) = [];
+dispz(1) = [];
+D(:,1) = [];
+
+%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-
+function M = getM(fi,mu)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% function getM
+%% Convert fault parameters into form needed for
+%% Cervelli calculations
+%%
+%% Inputs: fi fault stuff
+%% mu Poissons ratio
+%%
+%% Outputs: M matrix with fault stuff
+%% (see below for what is on
+%% the rows, 1 column per
+%% fault patch)
+%%
+%% Notes: This could be redone to be vectorized and faster
+%% (later maybe)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M=[];
+
+%% Loop over all the fault patches to build M
+s = size(fi);
+for j = 1:1:s(1)
+ fx1 = fi(j,1);
+ fy1 = fi(j,2);
+ fx2 = fi(j,3);
+ fy2 = fi(j,4);
+ dip = fi(j,5);%%in degs here
+ D = fi(j,6);
+ Pr = fi(j,7);
+ bd = fi(j,8);
+ ss = fi(j,9);
+ ds = fi(j,10);
+ ts = fi(j,11);
+ [strike, L, W, ofx, ofy, ofxe, ofye, tfx, tfy, tfxe, tfye] =...
+ fault_params_to_okada_form(fx1,fy1,fx2,fy2,dip*pi/180.0,D,bd);
+
+ length = L;
+ width = W;
+
+ if (dip) < 0
+ east = (L/2)*cos(strike) + tfx; %% east component of top center
+ north = (L/2)*sin(strike) + tfy; %% north " " " "
+ depth = bd;
+ width = abs(width);
+ ss = -1*ss;
+ ds = -1*ds;
+ ts = -1*ts;
+ else
+ east = (L/2)*cos(strike) + ofx; %% east component of bottom center
+ north = (L/2)*sin(strike) + ofy;
+ depth = D;
+ end
+
+ str = 180.0/pi*(pi/2-strike); %%Cervelli takes strike in degs
+ %%(strike comes out of
+ %% fault_params_to_okada_form in rads)
+
+ %%nextMcol = [length;width;depth;dip;str;east;north;ss;ds;ts];
+ M(:,j) = [length;width;depth;dip;str;east;north;ss;ds;ts];
+
+end
+%%
+%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-%-
+function [strike, L, W, ofx, ofy, ofxe, ofye, tfx, tfy, tfxe, tfye] = fault_params_to_okada_form(sx1, sy1, sx2, sy2, dip, D, bd)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% %%
+%% fault_params_to_okada_form.m %%
+%% %%
+%% MATLAB script %%
+%% %%
+%% This function takes fault trace, dip, and %%
+%% locking depth information and calculates the %%
+%% anchor coordinates, length, width and strike %%
+%% of the fault plane as per Okada (1985). %%
+%% %%
+%% It may be beneficial in the future to compute %%
+%% the midpoint base anchor as well for using %%
+%% Okada (1992) Stanford binary. %%
+%% %%
+%% We should also allow for buried faults in the %%
+%% future. This can be achieved by passing a %%
+%% locking depth and a burial depth. The only %%
+%% output changed would be the width of the %%
+%% fault plane. %%
+%% %%
+%% We may need to be more formal about say %%
+%% endpoint1 < endpoint2 ... maybe %%
+%% %%
+%% Arguments: %%
+%% sx1: x coord of fault trace endpoint1 %%
+%% sy1: y coord of fault trace endpoint1 %%
+%% sx2: x coord of fault trace endpoint2 %%
+%% sy2: y coord of fault trace endpoint2 %%
+%% dip: dip of fault plane [radians] %%
+%% D : fault locking depth %%
+%% bd : burial depth (top "locking depth") %%
+%% %%
+%% Returned variables: %%
+%% strike: stike of fault plane %%
+%% L : fault length %%
+%% W : fault width %%
+%% ofx : x coord of fault anchor %%
+%% ofy : y coord of fault anchor %%
+%% ofxe : x coord of other buried corner %%
+%% ofye : y coord of other buried corner %%
+%% tfx : x coord of fault anchor %%
+%% (top relative) %%
+%% tfy : y coord of fault anchor %%
+%% (top relative) %%
+%% tfxe : x coord of other buried corner %%
+%% (top relative) %%
+%% tfye : y coord of other buried corner %%
+%% (top relative) %%
+%% %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Do some basic checks %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (D < bd)
+ disp(sprintf('Burial depth is greater than locking depth!'))
+ disp(sprintf('Fault width will be negative!'))
+end
+
+if (D == bd)
+ disp(sprintf('Burial depth is equal to locking depth!'))
+ disp(sprintf('Fault width will zero!'))
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Calculate need parameters for Okada input %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+strike = atan2(sy1 - sy2, sx1 - sx2) + pi;
+L = sqrt((sx2 - sx1)^2+(sy2 - sy1)^2);
+W = (D - bd) / sin(dip);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% For older versions without a burial depth option %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% W = D / sin(dip);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Calculate fault segment anchor and other buried point %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ofx = sx1 + D / tan(dip) * sin(strike);
+ofy = sy1 - D / tan(dip) * cos(strike);
+ofxe = sx2 + D / tan(dip) * sin(strike);
+ofye = sy2 - D / tan(dip) * cos(strike);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Calculate fault segment anchor and other buried point (top relative) %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+tfx = sx1 + bd / tan(dip) * sin(strike);
+tfy = sy1 - bd / tan(dip) * cos(strike);
+tfxe = sx2 + bd / tan(dip) * sin(strike);
+tfye = sy2 - bd / tan(dip) * cos(strike);
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/points0.h5.gz
===================================================================
(Binary files differ)
Property changes on: cs/cigma/trunk/sandbox/okada/disloc3d/points0.h5.gz
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: cs/cigma/trunk/sandbox/okada/disloc3d/points0.in.gz
===================================================================
(Binary files differ)
Property changes on: cs/cigma/trunk/sandbox/okada/disloc3d/points0.in.gz
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,257 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "txtarray.h"
+
+
+int txtarray_read1(const char *filename,
+ double **data, int *m, int *n)
+{
+ FILE *fp;
+ char linebuf[500]; // 500 > 14 * (3 + 3 + 9 + 6)
+
+ int i,j;
+ double *x;
+
+ fp = fopen(filename, "r");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ if (fscanf(fp, "# %d %d\n", m, n) != 2)
+ {
+ return -2;
+ }
+
+ fgets(linebuf, 500, fp); // read header line
+
+ *data = (double *)malloc(((*m)*(*n))*sizeof(double));
+
+ for (i = 0; i < *m; i++)
+ {
+ x = &(*data)[(*n) * i];
+ for (j = 0; j < *n; j++)
+ {
+ fscanf(fp, "%lf", &x[j]);
+ }
+ }
+
+ fclose(fp);
+
+ return 0;
+}
+
+
+int txtarray_read_stations(const char *filename,
+ double **stations,
+ int *nstations)
+{
+ int three;
+ int ret;
+
+ ret = txtarray_read1(filename, stations, nstations, &three);
+
+ if (three != 3)
+ {
+ return -1;
+ }
+
+ return ret;
+}
+
+
+/*
+ * TODO: can this family of functions be written using varargs?
+ *
+ * txtarray_write(const char *filename,
+ * const char *header,
+ * int m,
+ * double *data1, int n1,
+ * double *data2, int n2,
+ * ...);
+ *
+ */
+
+int txtarray_write1(const char *filename, const char *header,
+ double *data, int m, int n)
+{
+ FILE *fp;
+ int i,j;
+ double *x;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ fprintf(fp, "# %d %d\n", m, n);
+
+ if (header != NULL)
+ fprintf(fp, "# %s\n", header);
+ else
+ fprintf(fp, "# header\n");
+
+ for (i = 0; i < m; i++)
+ {
+ x = &data[n*i];
+ for (j = 0; j < n; j++)
+ {
+ fprintf(fp, " %0.12e", x[j]);
+ }
+ fprintf(fp, "\n");
+ }
+ fclose(fp);
+
+ return 0;
+}
+
+int txtarray_write2(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2)
+{
+ FILE *fp;
+ int i,j;
+ double *x;
+ double *y;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ fprintf(fp, "# %d %d %d\n", m, n1, n2);
+
+ if (header != NULL)
+ fprintf(fp, "# %s\n", header);
+ else
+ fprintf(fp, "# header\n");
+
+ for (i = 0; i < m; i++)
+ {
+ x = &data1[n1 * i];
+ for (j = 0; j < n1; j++)
+ {
+ fprintf(fp, " %0.12e", x[j]);
+ }
+
+ y = &data2[n2 * i];
+ for (j = 0; j < n2; j++)
+ {
+ fprintf(fp, " %0.12e", y[j]);
+ }
+
+ fprintf(fp, "\n");
+ }
+ fclose(fp);
+
+ return 0;
+}
+
+int txtarray_write4(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2,
+ double *data3, int n3,
+ double *data4, int n4)
+{
+ FILE *fp;
+ int i,j;
+ double *x, *y, *z, *w;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ fprintf(fp, "# %d %d %d %d %d\n", m, n1, n2, n3, n4);
+
+ if (header != NULL)
+ fprintf(fp, "# %s\n", header);
+ else
+ fprintf(fp, "# header\n");
+
+ for (i = 0; i < m; i++)
+ {
+ x = &data1[n1 * i];
+ for (j = 0; j < n1; j++)
+ {
+ fprintf(fp, " %0.12e", x[j]);
+ }
+
+ y = &data2[n2 * i];
+ for (j = 0; j < n2; j++)
+ {
+ fprintf(fp, " %0.12e", y[j]);
+ }
+
+ z = &data3[n3 * i];
+ for (j = 0; j < n3; j++)
+ {
+ fprintf(fp, " %0.12e", z[j]);
+ }
+
+ w = &data4[n4 * i];
+ for (j = 0; j < n4; j++)
+ {
+ fprintf(fp, " %0.12e", w[j]);
+ }
+
+ fprintf(fp, "\n");
+ }
+ fclose(fp);
+
+ return 0;
+}
+
+int txtarray_write_disp(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements)
+{
+ return txtarray_write2(filename, header, nstations,
+ stations, 3, displacements, 3);
+}
+
+int txtarray_write_deriv(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *derivatives)
+{
+ return txtarray_write2(filename, header, nstations,
+ stations, 3, derivatives, 9);
+}
+
+int txtarray_write_stress(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *stress)
+{
+ return txtarray_write2(filename, header, nstations,
+ stations, 3, stress, 6);
+}
+
+int txtarray_write_all(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements,
+ double *derivatives,
+ double *stress)
+{
+ return txtarray_write4(filename, header,
+ nstations,
+ stations, 3,
+ displacements, 3,
+ derivatives, 9,
+ stress, 6);
+}
+
Added: cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/disloc3d/txtarray.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,56 @@
+#ifndef __TXTARRAY_H__
+#define __TXTARRAY_H__
+
+
+int txtarray_read1(const char *filename,
+ double **data, int *m, int *n);
+
+int txtarray_read_stations(const char *filename,
+ double **stations,
+ int *nstations);
+
+
+int txtarray_write1(const char *filename, const char *header,
+ double *data, int m, int n);
+
+int txtarray_write2(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2);
+
+int txtarray_write4(const char *filename,
+ const char *header,
+ int m,
+ double *data1, int n1,
+ double *data2, int n2,
+ double *data3, int n3,
+ double *data4, int n4);
+
+int txtarray_write_disp(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements);
+
+int txtarray_write_deriv(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *derivatives);
+
+int txtarray_write_stress(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements);
+
+int txtarray_write_all(const char *filename,
+ const char *header,
+ int nstations,
+ double *stations,
+ double *displacements,
+ double *derivatives,
+ double *stress);
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/furls
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/furls 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/furls 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1 @@
+http://richter.colorado.edu/okada/
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.c
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.c 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.c 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,111 @@
+#include <math.h>
+#include "okada.h"
+
+/* ******************************************************************
+
+ okada.c copyright (c) 8/96 by Seth McGinnis.
+
+ This is the primary file for the code that implements Okada's
+ analytic solutions to the problem of the displacement and stress
+ produced by a rectangular displacement in a linear elastic
+ half-space. All that this file does is to include some explanatory
+ material and to #include sub-header files, which contain the actual
+ code, distributed in such a fashion as to (hopefully) be
+ comprehensible as a whole.
+
+ All the mathematical content of these files comes from Yoshimitsu
+ Okada's paper "Internal Deformation Due to Shear and Tensile Faults
+ in a Half-Space", Bulletin of the Seismological Society of America,
+ Vol. 82, No. 2, pp 1018-1040, April 1992. Okada is a dude.
+
+ Some notes on the functioning of the program: The mathematics
+ assume the fault to be a rectangle with one corner located at
+ (0,0,-c). The strike of the fault is along the x-axis, and the
+ fault has a dip angle delta. Displacement on the fault can be
+ strike-slip, dip-slip, or tensile. The convention in this program
+ is that, for a 3-d slip vector, the first component is strike-slip,
+ the second dip-slip, and the third tensile. Therefore, if you want
+ to do things in some absolute reference frame, it's your
+ responsibility to convert your observation points and slip-vectors
+ to and from the assumed reference frame of this code. Note also
+ that angles are ALWAYS assumed to be in radians.
+
+ **************************************************************** */
+
+
+static const double TWOPI=2*M_PI;
+
+/* a convenience variable. M_PI is double-precision and lives in
+ <math.h>. */
+
+
+static double lambda = 1.0;
+static double mu = 1.0;
+static double alpha = 0.666666666666666666667;
+
+#include "okada.lame.h"
+
+/* These are the lame' parameters. Alpha is defined to be
+ (lambda+mu)/(lambda +2*mu). Access functions for these variables
+ are given in the lame.h and are prototyped in the main okada.h
+ header file. */
+
+
+#include "okada.notation.h"
+
+/* Many of the quantities in Okada's equations are shorthand for more
+ complicated expressions; these have been implemented for the most
+ part with #define expressions in the notation.h file. The
+ 'chinnery' function (see p. 1026) is also in this file. */
+
+
+#include "okada.disp.h"
+
+/* The functions which are used to calculate the displacement at a
+ given observation point due to some slip on the fault. */
+
+
+#include "okada.strain.x.h"
+
+/* The functions which are used to calculate the x-derivatives of
+ displacement at a given observation point due to some slip on the
+ fault. Strains are used to calculate stress; details are in
+ response.h. */
+
+
+#include "okada.strain.y.h"
+
+/* The functions which are used to calculate the y-derivatives of
+ displacement at a given observation point due to some slip on the
+ fault. */
+
+
+#include "okada.strain.z.h"
+
+/* The functions which are used to calculate the z-derivatives of
+ displacement at a given observation point due to some slip on the
+ fault. */
+
+
+#include "okada.response.h"
+
+/* Top-level functions that will calculate displacement and/or
+ stresses at an observation point due to slip on a fault. These
+ functions are prototyped in okada.h; this is the header file to
+ look at for finding out how the Okada code tells you what your
+ displacements/strains are. */
+
+
+/* #include "okada.self.h" */
+
+/* Function for calculating the Kself matrix terms, which tell you
+ about the shear stress at the center of a fault due to slip on the
+ fault itself. This is very useful in slider-block type code. Also
+ prototyped in okada.h */
+
+
+#include "okada.cleanup.h"
+
+/* #undefs all the things #defined in notation.h; this is necessary
+ because a lot of the #defined variables have names like 'h' and 'q'
+ and things that you might expect to show up elsewhere. */
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.cleanup.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.cleanup.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.cleanup.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,76 @@
+#ifndef OKADA_CLEANUP_H
+#define OKADA_CLEANUP_H
+
+/* This is the header file that cleans up the #defines, to prevent
+ name clashing. */
+
+#undef d
+#undef p
+#undef q
+#undef ybar
+#undef dbar
+#undef cbar
+
+
+#undef X11
+#undef X32
+#undef X53
+#undef Y11
+#undef Y32
+#undef Y53
+
+#undef h
+#undef Z32
+#undef Z53
+
+#undef Y0
+#undef Z0
+
+
+
+#undef lnRxi
+#undef lnReta
+
+#undef Theta
+#undef Xsq
+#undef X
+
+#undef I3
+#undef I4
+#undef I1
+#undef I2
+
+
+
+#undef D11
+#undef K1
+#undef K3
+#undef K2
+#undef K4
+#undef J2
+#undef J3
+#undef J4
+#undef J5
+#undef J6
+#undef J1
+
+
+
+#undef E
+#undef F
+#undef G
+#undef H
+#undef P
+#undef Q
+
+
+
+#undef Eprime
+#undef Fprime
+#undef Gprime
+#undef Hprime
+#undef Pprime
+#undef Qprime
+
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.disp.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.disp.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.disp.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,404 @@
+#ifndef OKADA_DISP_H
+#define OKADA_DISP_H
+
+/* This file contains code for calculating displacements */
+
+
+
+/*************************************************************
+ A/B/C
+this part of the file is the f s, the entries in Table 6.
+ i
+
+(notation: f1As, for instance, is the 1st (x) component of fA
+for the strike-slip case.)
+**************************************************************/
+
+
+
+
+
+double f1As(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( Theta/2 + alpha/2 * xi * q * Y11 );
+}
+
+double f1Bs(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return ( -xi*q*Y11 - Theta - (1-alpha)/alpha * I1 * sindel);
+}
+
+double f1Cs(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)*xi*Y11*cosdel - alpha * xi * q * Z32);
+}
+
+
+
+double f2As(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( alpha/2 * q/R);
+}
+
+double f2Bs(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -q/R + (1-alpha)/alpha * ybar/(R+dbar) * sindel);
+}
+
+double f2Cs(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)*(cosdel/R + 2*q*Y11*sindel) - alpha*cbar*q/(R*Rsq));
+}
+
+
+
+double f3As(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)/2 * lnReta - alpha/2 * q*q*Y11);
+}
+
+double f3Bs(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return ( q*q*Y11 - (1-alpha)/alpha * I2 * sindel);
+}
+
+double f3Cs(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)*q*Y11*cosdel - alpha*(cbar*eta/(Rsq*R) - z*Y11 + xi*xi*Z32));
+}
+
+
+
+
+
+double f1Ad(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( alpha/2 * q/R);
+}
+
+double f1Bd(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return ( -q/R + (1-alpha)/alpha * I3 * sindel * cosdel);
+}
+
+double f1Cd(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)*cosdel/R - q*Y11*sindel - alpha*cbar*q/(R*Rsq));
+}
+
+
+
+double f2Ad(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( Theta/2 + alpha/2 * eta*q*X11);
+}
+
+double f2Bd(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -eta*q*X11 - Theta - (1-alpha)/alpha * xi/(R+dbar) * sindel*cosdel);
+}
+
+double f2Cd(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)*ybar*X11 - alpha*cbar*eta*q*X32);
+}
+
+
+
+double f3Ad(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)/2 * lnRxi - alpha/2 * q*q*X11);
+}
+
+double f3Bd(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( q*q*X11 + (1-alpha)/alpha * I4 * sindel * cosdel);
+}
+
+double f3Cd(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -dbar*X11 - xi*Y11*sindel - alpha*cbar * (X11 - q*q*X32));
+}
+
+
+
+
+
+double f1At(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -(1-alpha)/2 * lnReta - alpha/2 * q*q*Y11);
+}
+
+double f1Bt(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( q*q*Y11 - (1-alpha)/alpha * I3 * sindel * sindel);
+}
+
+double f1Ct(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -(1-alpha)*(sindel/R + q*Y11*cosdel) - alpha*(z*Y11 - q*q*Z32));
+}
+
+
+
+double f2At(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -(1-alpha)/2 * lnReta - alpha/2 * q*q*X11);
+}
+
+double f2Bt(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( q*q*X11 + (1-alpha)/alpha * xi/(R+dbar) * sindel * sindel);
+}
+
+double f2Ct(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha)*2*xi*Y11*sindel + dbar*X11 - alpha*cbar*(X11 - q*q*X32));
+}
+
+
+
+double f3At(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( Theta/2 - alpha/2 * q * (eta*X11 + xi*Y11));
+}
+
+double f3Bt(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( q*(eta*X11+xi*Y11) - Theta - (1-alpha)/alpha * I4 * sindel * sindel);
+}
+
+double f3Ct(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( (1-alpha) * (ybar*X11 + xi*Y11*cosdel)
+ + alpha * q * (cbar*eta*X32 + xi*Z32));
+}
+
+
+
+/********************************************************************/
+/* These displacement-component functions are at the top of Table 6 */
+
+
+double ux(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ chinnery(f1As,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f1As,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f1Bs,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f1Cs,x,y,z,c,L,W,cosdel,sindel))
+ + Ud*(
+ chinnery(f1Ad,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f1Ad,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f1Bd,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f1Cd,x,y,z,c,L,W,cosdel,sindel))
+ + Ut*(
+ chinnery(f1At,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f1At,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f1Bt,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f1Ct,x,y,z,c,L,W,cosdel,sindel));
+
+ return(ans/TWOPI);
+}
+
+
+
+double uy(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ cosdel*( chinnery(f2As,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2As,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f2Bs,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f2Cs,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(f3As,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3As,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f3Bs,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f3Cs,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ cosdel*( chinnery(f2Ad,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2Ad,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f2Bd,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f2Cd,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(f3Ad,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3Ad,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f3Bd,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f3Cd,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ cosdel*( chinnery(f2At,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2At,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f2Bt,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f2Ct,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(f3At,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3At,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f3Bt,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(f3Ct,x,y,z,c,L,W,cosdel,sindel)));
+
+ return(ans/TWOPI);
+}
+
+
+
+double uz(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ sindel*( chinnery(f2As,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2As,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f2Bs,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(f2Cs,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(f3As,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3As,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f3Bs,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(f3Cs,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ sindel*( chinnery(f2Ad,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2Ad,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f2Bd,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(f2Cd,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(f3Ad,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3Ad,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f3Bd,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(f3Cd,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ sindel*( chinnery(f2At,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2At,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f2Bt,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(f2Ct,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(f3At,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3At,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(f3Bt,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(f3Ct,x,y,z,c,L,W,cosdel,sindel)));
+
+ return(ans/TWOPI);
+}
+
+
+
+
+#endif
+
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,107 @@
+#ifndef OKADA_H
+#define OKADA_H
+
+/*******************************************************************
+
+ okada.h copyright (c) 8/96 by Seth McGinnis.
+
+ This is the header file for the okada library, which implements
+ Okada's analytic solutions to the problem of the displacement,
+ strain, and/or stress produced by a rectangular displacement in a
+ linear elastic half-space.
+
+ All the mathematical content of these files comes from Yoshimitsu
+ Okada's paper "Internal Deformation Due to Shear and Tensile Faults
+ in a Half-Space", Bulletin of the Seismological Society of America,
+ Vol. 82, No. 2, pp 1018-1040, April 1992. Okada is a dude.
+
+ Some notes on the functioning of the code: The mathematics assume
+ the fault to be a rectangle with one corner located at (0,0,-c).
+ The strike of the fault is along the x-axis, and the fault has a dip
+ angle delta. Displacement on the fault can be strike-slip,
+ dip-slip, or tensile. The convention in this program is that, for a
+ 3-d slip vector, the first component is strike-slip, the second
+ dip-slip, and the third tensile. Therefore, if you want to do
+ things in some absolute reference frame, it's your responsibility to
+ convert your observation points and slip-vectors to and from the
+ assumed reference frame of this code. Angles are ALWAYS assumed to
+ be in radians.
+
+ ******************************************************************/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+ /* okada.c uses the lame' parameters, lambda and mu. It also uses
+ alpha, which is defined as alpha = (lambda+mu)/(lambda+2*mu). By
+ default, lambda and mu are 1.0, which makes alpha = 2/3. These
+ access functions allow other pieces of code to inspect and change
+ the values. If you change lambda or mu, alpha will be adjusted
+ so that the three values remain consistent. The "SetLame"
+ function sets both lame parameters simultaneously. Note,
+ however, that no checking is done to ensure that values are
+ realistic (i.e., nonnegative, etc.) Consider yourself warned. */
+
+ double OkadaLambda(void);
+ double OkadaMu(void);
+ double OkadaAlpha(void);
+
+ void OkadaSetLambda(double newval);
+ void OkadaSetMu(double newval);
+ void OkadaSetLame(double lambda, double mu);
+
+
+ /*****************************************************************/
+
+
+ /* Functions that calculate displacement and/or stress at an
+ observation point due to slip on a fault due to a given fault
+ displacement. DispResponse() returns the displacements;
+ StrainResponse() returns the strains, StressResponse() returns
+ the stresses as calculated from the strains, and just plain
+ response() returns the displacements and the stresses. Stress is
+ always in the x=strike, z=up frame. */
+
+ /* slip, obspoint, and disp are always x[3]; stress and strain are
+ tensors and should always be x[9]. (recording convention:
+
+ s = 0 1 2
+ 3 4 5
+ 6 7 8 ) */
+
+
+ void response(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *disp, double *stress);
+
+ void DispResponse(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *disp);
+
+ void StrainResponse(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *strain);
+
+ void StressResponse(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *stress);
+
+
+ /* function for calculating the self terms-- shear stress at the
+ center of the fault due to displacement of the fault. Strike-slip
+ generates only shear stress in the strike direction; dip-slip
+ generates only dip-wise shear stress. Selfstrike and selfdip are
+ the strike- and dip- shear stress (components 31 and 32 of stress
+ in the frame of the fault) due to a unit displacement. */
+/*
+ void SelfMatrix(double c, double dip, double length, double width,
+ double *selfstrike, double *selfdip);
+*/
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.lame.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.lame.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.lame.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,48 @@
+#ifndef OKADA_LAME_H
+#define OKADA_LAME_H
+
+/* This header file contains the access functions for the
+ (library-file scope) global variables alpha, lambda, and mu, which
+ determine the material properties of the elastic medium. Names
+ should be self-explanatory. The relationship amongst the variables
+ is: alpha = (lambda+mu)/(lambda+2*mu). Default values for lambda,
+ mu, and alpha are 1, 1, and 2/3, respectively. */
+
+double OkadaLambda(void)
+{
+ return(lambda);
+}
+
+double OkadaMu(void)
+{
+ return(mu);
+}
+
+double OkadaAlpha(void)
+{
+ return (alpha);
+}
+
+void OkadaSetLambda(double newval)
+{
+ lambda = newval;
+ alpha = (lambda+mu)/(lambda+2*mu);
+ return;
+}
+
+void OkadaSetMu(double newval)
+{
+ mu = newval;
+ alpha = (lambda+mu)/(lambda+2*mu);
+ return;
+}
+
+void OkadaSetLame(double new_lambda, double new_mu)
+{
+ lambda = new_lambda;
+ mu = new_mu;
+ alpha = (lambda+mu)/(lambda+2*mu);
+ return;
+}
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.notation.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.notation.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.notation.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,178 @@
+#ifndef OKADA_NOTATION_H
+#define OKADA_NOTATION_H
+
+/* This is the header file that does most of the "notational" parts of
+ the Okada math. */
+
+
+
+/********************/
+/* general things */
+
+
+
+/* for numerics and if(x==0) statements */
+
+static const double EPS = 0.00001;
+
+
+/* from tables 6-9 in Okada */
+#define d (c - z)
+#define p (y*cosdel + d*sindel)
+#define q (y*sindel - d*cosdel)
+#define ybar (eta*cosdel + q*sindel)
+#define dbar (eta*sindel - q*cosdel)
+#define cbar (dbar + z)
+
+
+
+/********************************************************/
+
+/* Chinnery's notation (p1026),
+ ||f(xi,eta) = f(x,p) - f(x,p-W) - f(x-L,p) + f(x-L,p-W)
+ is implemented by this function:
+ */
+
+inline double chinnery(double (*f)(double,double,double,double,
+ double,double,double),
+ double x, double y, double z,
+ double c, double L, double W,
+ double cosdel, double sindel)
+{
+ return( (*f)( x , p , y, z, c, cosdel, sindel)
+ - (*f)( x , p-W, y, z, c, cosdel, sindel)
+ - (*f)(x-L, p , y, z, c, cosdel, sindel)
+ + (*f)(x-L, p-W, y, z, c, cosdel, sindel));
+}
+
+
+
+/*************************/
+/* more general things */
+
+/* from p. 1026 of Okada */
+
+#define X11 ((fabs(R+xi)<EPS)?0: \
+ 1/(R*(R+xi)))
+#define X32 ((fabs(R+xi)<EPS)?0: \
+ (2*R+xi)/(Rsq*R*(R+xi)*(R+xi)))
+#define X53 ((fabs(R+xi)<EPS)?0: \
+ (8*Rsq+9*R*xi+3*xi*xi)/(Rsq*Rsq*R*(R+xi)*(R+xi)*(R+xi)))
+
+#define Y11 ((fabs(R+eta)<EPS)?0: \
+ 1/(R*(R+eta)))
+#define Y32 ((fabs(R+eta)<EPS)?0: \
+ (2*R+eta)/(Rsq*R*(R+eta)*(R+eta)))
+#define Y53 ((fabs(R+eta)<EPS)?0: \
+ (8*Rsq+9*R*eta+3*eta*eta)/(Rsq*Rsq*R*(R+eta)*(R+eta)*(R+eta)))
+
+#define h (q*cosdel - z)
+#define Z32 (sindel/(Rsq*R) - h*Y32)
+#define Z53 (3*sindel/(Rsq*Rsq*R) - h*Y53)
+
+#define Y0 (Y11 - xi*xi*Y32)
+#define Z0 (Z32 - xi*xi*Z53)
+
+
+
+/***************************************************/
+/* Displacement-specific notation (from Table 6) */
+
+#define lnRxi ((fabs(R+xi)<EPS) ? -log(R-xi) : log(R+xi))
+#define lnReta ((fabs(R+eta)<EPS) ? -log(R-eta) : log(R+eta))
+
+#define Theta ((fabs(q)<EPS) ? 0 : atan(xi*eta/(q*R)))
+#define Xsq (xi*xi+q*q)
+#define X (sqrt(Xsq))
+
+
+#define I3 ((fabs(cosdel)<EPS) ? \
+ ((eta/(R+dbar)+ybar*q/((R+dbar)*(R+dbar))-lnReta)/2.0) : \
+ ((ybar/(R+dbar)-(lnReta-sindel*log(R+dbar))/cosdel)/cosdel))
+
+#define I4 ((fabs(cosdel)<EPS) ? \
+ (xi*ybar/(2.0*(R+dbar)*(R+dbar))) : \
+ ( (sindel*xi/(R+dbar) \
+ +2/cosdel*atan( \
+ (eta*(X+q*cosdel)+X*(R+X)*sindel) \
+ /(xi*(R+X)*cosdel))) \
+ /cosdel))
+
+#define I1 (-xi*cosdel/(R+dbar) - I4*sindel)
+
+#define I2 (log(R+dbar)+I3*sindel)
+
+
+
+
+/******************************************/
+/* X-strain specific notation (Table 7) */
+
+#define D11 (1/(R*(R+dbar)))
+
+#define K1 ((fabs(cosdel)<EPS) ? \
+ (xi*q/(R+dbar)*D11) : \
+ (xi/cosdel*(D11-Y11*sindel)))
+
+#define K3 ((fabs(cosdel)<EPS) ? \
+ (sindel/(R+dbar)*(xi*xi*D11-1)) : \
+ ((q*Y11-ybar*D11)/cosdel))
+
+#define K2 (1/R+K3*sindel)
+
+#define K4 (xi*Y11*cosdel - K1*sindel)
+
+#define J2 (xi*ybar*D11/(R+dbar))
+
+#define J3 ((fabs(cosdel)<EPS) ? \
+ (-xi/((R+dbar)*(R+dbar))*(q*q*D11-1/2.0)) : \
+ ((K1-J2*sindel)/cosdel))
+
+#define J4 (-xi*Y11-J2*cosdel+J3*sindel)
+
+#define J5 (-(dbar+ybar*ybar/(R+dbar))*D11)
+
+#define J6 ((fabs(cosdel)<EPS) ? \
+ (-ybar/((R+dbar)*(R+dbar))*(xi*xi*D11-1/2.0)) : \
+ (1/cosdel*(K3-J5*sindel)))
+
+#define J1 (J5*cosdel-J6*sindel)
+
+
+
+/******************************************/
+/* Y-strain specific notation (Table 8) */
+
+#define E (sindel/R - ybar*q/(R*Rsq))
+
+#define F (dbar/(R*Rsq) + xi*xi*Y32*sindel)
+
+#define G (2*X11*sindel - ybar*q*X32)
+
+#define H (dbar*q*X32 + xi*q*Y32*sindel)
+
+#define P (cosdel/(R*Rsq) + q*Y32*sindel)
+
+#define Q (3*cbar*dbar/(R*Rsq*Rsq) - (z*Y32 + Z32 +Z0)*sindel)
+
+
+
+/******************************************/
+/* Z-strain specific notation (Table 9) */
+
+#define Eprime (cosdel/R + dbar*q/(R*Rsq))
+
+#define Fprime (ybar/(R*Rsq) + xi*xi*Y32*cosdel)
+
+#define Gprime (2*X11*cosdel + dbar*q*X32)
+
+#define Hprime (ybar*q*X32 + xi*q*Y32*cosdel)
+
+#define Pprime (sindel/(R*Rsq) - q*Y32*cosdel)
+
+#define Qprime (3*cbar*ybar/(R*Rsq*Rsq) + q*Y32 - (z*Y32 + Z32 + Z0)*cosdel)
+
+
+
+#endif
+
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.response.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.response.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.response.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,219 @@
+#ifndef OKADA_RESPONSE_H
+#define OKADA_RESPONSE_H
+
+/* This header file contains the functions that call all the functions
+ in the other files-- you use these functions to actually find out
+ what the stress or displacement is at an observation point due to a
+ given fault displacement. disp_response() returns the
+ displacements; strain_response() returns the strains,
+ stress_response() returns the stresses as calculated from the
+ strains, and just plain response() returns the displacements and
+ the stresses.*/
+
+
+void response(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *disp, double *stress)
+{
+ double strain[3][3];
+
+
+ /* calculate displacements */
+
+ disp[0] = ux(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ disp[1] = uy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ disp[2] = uz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+
+ /* calculate strains */
+
+ strain[0][0] = duxdx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[0][1] = duydx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[0][2] = duzdx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ strain[1][0] = duxdy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[1][1] = duydy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[1][2] = duzdy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ strain[2][0] = duxdz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[2][1] = duydz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[2][2] = duzdz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ /* stress comes from strain, as follows: (Okada p. 1024)
+ e(i,j) = 1/2(dui/dxj + duj/dxi)
+ sigma(i,j) = lambda e(k,k) delta(i,j) + 2 mu e(i,j)
+
+ Note that e(i,j) == e(j,i), and similarly with sigma.
+ (However, dui/dxj != duj/dxi, so we really did have to do all
+ those strain function calls up there.) We can now calculate and
+ return the stresses...
+ */
+
+ /* sigma11 = lambda(du1/dux1+du2/dx2+du3/dx3) + 2 mu du1/dx1 */
+ /* sigma22 = lambda(du1/dux1+du2/dx2+du3/dx3) + 2 mu du2/dx2 */
+ /* sigma33 = lambda(du1/dux1+du2/dx2+du3/dx3) + 2 mu du3/dx3 */
+
+ /* sigma12 = sigma21 = mu(du1/dx2 + du2/dx1) */
+ /* sigma23 = sigma32 = mu(du2/dx3 + du3/dx2) */
+ /* sigma31 = sigma13 = mu(du3/dx1 + du1/dx3) */
+
+ /* (stress[0...8] = s11, s12, s13, s21, s22, s23, s31, s32, s33) */
+ /* stress[0] serves as a temporary holder for epsilon(k,k) */
+
+
+ stress[0] = lambda*(strain[0][0]+strain[1][1]+strain[2][2]);
+ stress[1] = mu*(strain[0][1] + strain[1][0]);
+ stress[2] = mu*(strain[2][0] + strain[0][2]);
+ stress[3] = stress[1];
+ stress[4] = stress[0] + 2*mu*strain[1][1];
+ stress[5] = mu*(strain[1][2] + strain[2][1]);
+ stress[6] = stress[2];
+ stress[7] = stress[5];
+ stress[8] = stress[0] + 2*mu*strain[2][2];
+
+ stress[0] += 2*mu*strain[0][0];
+
+ return;
+}
+
+
+
+
+void DispResponse(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *disp)
+
+{
+
+ /* calculate displacements */
+
+ disp[0] = ux(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ disp[1] = uy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ disp[2] = uz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ return;
+}
+
+
+
+
+void StrainResponse(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *strain)
+{
+ /* calculate strains */
+
+ strain[0] = duxdx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[1] = duydx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[2] = duzdx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ strain[3] = duxdy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[4] = duydy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[5] = duzdy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ strain[6] = duxdz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[7] = duydz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[8] = duzdz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ return;
+}
+
+
+
+
+
+void StressResponse(double c, double dip, double length, double width,
+ const double *slip, const double *obspoint,
+ double *stress)
+{
+ double strain[3][3];
+
+
+ /* calculate strains */
+
+ strain[0][0] = duxdx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[0][1] = duydx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[0][2] = duzdx(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ strain[1][0] = duxdy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[1][1] = duydy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[1][2] = duzdy(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ strain[2][0] = duxdz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[2][1] = duydz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+ strain[2][2] = duzdz(obspoint[0], obspoint[1], obspoint[2], slip[0],
+ slip[1], slip[2], c, dip, length, width);
+
+ /* stress comes from strain, as follows: (Okada p. 1024)
+ epsilon(i,j) = 1/2(dui/dxj + duj/dxi)
+ sigma(i,j) = lambda epsilon(k,k) delta(i,j) + 2 mu epsilon(i,j)
+
+ Note that epsilon(i,j) == epsilon(j,i), and similarly with sigma.
+ (However, dui/dxj != duj/dxi, so we really did have to do all
+ those strain function calls up there.) We can now calculate and
+ return the stresses...
+ */
+
+ /* sigma11 = lambda(du1/dux1+du2/dx2+du3/dx3) + 2 mu du1/dx1 */
+ /* sigma22 = lambda(du1/dux1+du2/dx2+du3/dx3) + 2 mu du2/dx2 */
+ /* sigma33 = lambda(du1/dux1+du2/dx2+du3/dx3) + 2 mu du3/dx3 */
+
+ /* sigma12 = sigma21 = mu(du1/dx2 + du2/dx1) */
+ /* sigma23 = sigma32 = mu(du2/dx3 + du3/dx2) */
+ /* sigma31 = sigma13 = mu(du3/dx1 + du1/dx3) */
+
+ /* (stress[0...8] = s11, s12, s13, s21, s22, s23, s31, s32, s33) */
+ /* stress[0] serves as a temporary holder for epsilon(k,k) */
+
+ stress[0] = lambda*(strain[0][0]+strain[1][1]+strain[2][2]);
+ stress[1] = mu*(strain[0][1] + strain[1][0]);
+ stress[2] = mu*(strain[2][0] + strain[0][2]);
+ stress[3] = stress[1];
+ stress[4] = stress[0] + 2*mu*strain[1][1];
+ stress[5] = mu*(strain[1][2] + strain[2][1]);
+ stress[6] = stress[2];
+ stress[7] = stress[5];
+ stress[8] = stress[0] + 2*mu*strain[2][2];
+
+ stress[0] += 2*mu*strain[0][0];
+
+ return;
+}
+
+
+
+
+#endif
+
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.self.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.self.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.self.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,133 @@
+#ifndef OKADA_SELF_H
+#define OKADA_SELF_H
+
+/* This header file contains the functions used to calculate the
+ "self" term; the matrix that describes the (SHEAR) stresses induced
+ at the center of a fault patch by the motion of that fault patch.
+ Delightfully, physical intution suggests (and a little math and
+ experimentation verifies) that if we are only concerned with shear
+ stresses (which we are for modeling fault patches), the self matrix
+ is diagonal! Slip in the 1 (strike) direction will generate zero
+ stress in the sigma32 orientation; it will only generate stress in
+ the sigma31 orientation. Likewise for dip-slip. Therefore, there
+ are only two numbers which determine the shear stress produced by a
+ given (non-tensile) slip:
+
+ sigma31 = selfstrike * Ustrike;
+ sigma32 = selfstrike * Udip;
+
+ By generating these two numbers, the amount of slip needed to
+ relieve a given state of shear stress can be found.
+
+ The way to find these two numbers is to calculate the strain for a
+ unit displacement in the slip/dip direction, calculate the stress
+ resulting from that strain, rotate that stress into the frame of
+ the fault, and then consider the 31 and 32 elements of the stress.
+ Thus, we have the linear relation between displacement in the
+ slip/dip direction and stress on the fault plane in one of those
+ two directions. We can do this because the response of a given
+ stress component to a given displacement component is independent
+ of the other displacement components. And since the matrix turns
+ out to be diagonal (after rotation), we need only extract two
+ numbers.
+
+ */
+
+
+/* *************************************************************
+
+ Initially, okada.h #included "matrixrot.h" for this headerfile to
+ use the DipMatrixRotation function. However, since that function
+ and only that function is used here, and it is used here and
+ nowhere else, it has been duplicated here, to keep you from needing
+ to link to libmatrixrot.a when you #include "okada.h". Will you
+ ever use the okada code without using the matrixrot code also? Who
+ knows. Perhaps...
+
+ */
+
+#include <math.h> /* just to be sure */
+
+void OkadaDipMatrixRotation(double* in, double* out, double xangle)
+{
+ /* Rotation of a matrix around the xaxis only.
+
+ [1 0 0 ]
+ out = Rx.in.Rx-transpose, Rx = [0 cx sx]
+ [0 -sx cx]
+
+ */
+
+ double s, c, csq, ssq;
+
+ s = sin(xangle);
+ c = cos(xangle);
+ ssq = s*s;
+ csq = c*c;
+
+ out[0] = in[0];
+ out[1] = in[1]*c + in[2]*s;
+ out[2] = in[2]*c - in[1]*s;
+ out[3] = in[3]*c + in[6]*s;
+ out[4] = in[4]*csq + (in[5] + in[7])*c*s + in[8]*ssq;
+ out[5] = in[5]*csq + (in[8] - in[4])*c*s - in[7]*ssq;
+ out[6] = in[6]*c - in[3]*s;
+ out[7] = in[7]*csq + (in[8] - in[4])*c*s - in[5]*ssq;
+ out[8] = in[8]*csq - (in[5] + in[7])*c*s + in[4]*ssq;
+
+ return;
+}
+
+/***************************************************************/
+
+
+
+void SelfMatrix(double c, double dip, double length, double width,
+ double *selfstrike, double *selfdip)
+{
+ double strikestress[9], dipstress[9], temp1[9], temp2[9];
+
+ double center[3];
+ double unitdisp[3];
+
+ center[0] = length/2;
+ center[1] = width/2*cos(dip);
+ center[2] = width/2*sin(dip) - c;
+
+
+ /* get stress due to unit strike displacement... */
+
+ unitdisp[0] = 1;
+ unitdisp[1] = 0;
+ unitdisp[2] = 0;
+
+ StressResponse(c, dip, length, width, unitdisp, center, temp1);
+
+
+ /* ...and unit dip displacement */
+
+ unitdisp[0] = 0;
+ unitdisp[1] = 1;
+ unitdisp[2] = 0;
+
+ StressResponse(c, dip, length, width, unitdisp, center, temp2);
+
+ /* Now we rotate those stresses from strike-horiz-vert frame back
+ into the faultpatch frame (strike-dip-normal) */
+
+ OkadaDipMatrixRotation(temp1, strikestress, dip);
+ OkadaDipMatrixRotation(temp2, dipstress, dip);
+
+ /* and now we pull out selfstrike and selfdip. */
+
+ *selfstrike = strikestress[6];
+ *selfdip = dipstress[7];
+
+ /* All done! */
+
+ return;
+}
+
+
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.x.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.x.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.x.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,400 @@
+#ifndef OKADA_STRAIN_X_H
+#define OKADA_STRAIN_X_H
+
+/* This header file contains the functions used to calculate the
+ x-derivatives of displacement */
+
+
+
+/*********************************************************/
+/** And here we put the actual functions from Table 7.***/
+
+double df1Asdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return( -(1-alpha)/2.0*q*Y11 - alpha/2.0*xi*xi*q*Y32 );
+}
+
+double df1Bsdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(xi*xi*q*Y32 - (1-alpha)/alpha*J1*sindel);
+}
+
+double df1Csdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*Y0*cosdel - alpha*q*Z0);
+}
+
+
+
+
+double df2Asdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-alpha/2.0*xi*q/(Rsq*R));
+}
+
+double df2Bsdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(xi*q/(R*Rsq) - (1-alpha)/alpha*J2*sindel);
+}
+
+double df2Csdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*xi*(cosdel/(R*Rsq) + 2*q*Y32*sindel) + alpha*3*cbar*xi*q/(Rsq*Rsq*R));
+}
+
+
+
+
+double df3Asdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*xi*Y11 + alpha/2.0*xi*q*q*Y32);
+}
+
+double df3Bsdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-xi*q*q*Y32 - (1-alpha)/alpha*J3*sindel);
+}
+
+double df3Csdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*xi*q*Y32*cosdel + alpha*xi*(3*cbar*eta/(Rsq*Rsq*R)-z*Y32-Z32-Z0));
+}
+
+
+
+
+
+
+
+double df1Addx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-alpha/2.0*xi*q/(Rsq*R));
+}
+
+double df1Bddx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(xi*q/(R*Rsq) +(1-alpha)/alpha*J4*sindel*cosdel);
+}
+
+double df1Cddx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*xi*cosdel/(R*Rsq) + xi*q*Y32*sindel + alpha*3*cbar*xi*q/(Rsq*Rsq*R));
+}
+
+
+
+
+double df2Addx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-q/2.0*Y11 - alpha/2.0*eta*q/(R*Rsq));
+}
+
+double df2Bddx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(eta*q/(R*Rsq) + q*Y11 + (1-alpha)/alpha*J5*sindel*cosdel);
+}
+
+double df2Cddx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*ybar/(R*Rsq) + alpha*3*cbar*eta*q/(R*Rsq*Rsq));
+}
+
+
+
+
+double df3Addx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/(2.0*R) + alpha/2.0*q*q/(R*Rsq));
+}
+
+double df3Bddx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-q*q/(R*Rsq) + (1-alpha)/alpha*J6*sindel*cosdel);
+}
+
+double df3Cddx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(dbar/(R*Rsq) - Y0*sindel + alpha*cbar/(R*Rsq)*(1-3*q*q/Rsq));
+}
+
+
+
+
+
+
+
+double df1Atdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/2.0*xi*Y11 + alpha/2*xi*q*q*Y32);
+}
+
+double df1Btdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-xi*q*q*Y32 - (1-alpha)/alpha*J4*sindel*sindel);
+}
+
+double df1Ctdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*xi*sindel/(R*Rsq) + xi*q*Y32*cosdel + alpha*xi*(3*cbar*eta/(Rsq*Rsq*R)-2*Z32-Z0));
+}
+
+
+
+
+double df2Atdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/(2.0*R) + alpha*q*q/(2.0*Rsq*R));
+}
+
+double df2Btdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-q*q/(R*Rsq) - (1-alpha)/alpha*J5*sindel*sindel);
+}
+
+double df2Ctdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*2*Y0*sindel - dbar/(R*Rsq) + alpha*cbar/(R*Rsq)*(1-3*q*q/Rsq));
+}
+
+
+
+
+double df3Atdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/2.0*q*Y11 - alpha/2.0*q*q*q*Y32);
+}
+
+double df3Btdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*q*q*Y32 - (1-alpha)/alpha*J6*sindel*sindel);
+}
+
+double df3Ctdx(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*(ybar/(R*Rsq)-Y0*cosdel) - alpha*(3*cbar*eta*q/(Rsq*Rsq*R)-q*Z0));
+}
+
+
+/**********************************************************/
+/* x strain-component functions. (At the top of Table 7) */
+
+double duxdx(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ chinnery(df1Asdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df1Asdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Bsdx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Csdx,x,y,z,c,L,W,cosdel,sindel))
+
+ + Ud*(
+ chinnery(df1Addx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df1Addx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Bddx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Cddx,x,y,z,c,L,W,cosdel,sindel))
+
+ + Ut*(
+ chinnery(df1Atdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df1Atdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Btdx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Ctdx,x,y,z,c,L,W,cosdel,sindel));
+
+ return(ans/TWOPI);
+}
+
+double duydx(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ cosdel*( chinnery(df2Asdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Asdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bsdx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Csdx,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Asdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Asdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bsdx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Csdx,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ cosdel*( chinnery(df2Addx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Addx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bddx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Cddx,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Addx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Addx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bddx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Cddx,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ cosdel*( chinnery(df2Atdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Atdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Btdx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Ctdx,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Atdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Atdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Btdx,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Ctdx,x,y,z,c,L,W,cosdel,sindel)));
+ return(ans/TWOPI);
+}
+
+
+
+
+double duzdx(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ sindel*( chinnery(df2Asdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Asdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bsdx,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Csdx,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Asdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Asdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bsdx,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Csdx,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ sindel*( chinnery(df2Addx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Addx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bddx,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Cddx,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Addx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Addx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bddx,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Cddx,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ sindel*( chinnery(df2Atdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Atdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Btdx,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Ctdx,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Atdx,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Atdx,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Btdx,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Ctdx,x,y,z,c,L,W,cosdel,sindel)));
+
+ return(ans/TWOPI);
+}
+
+
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.y.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.y.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.y.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,431 @@
+#ifndef OKADA_STRAIN_Y_H
+#define OKADA_STRAIN_Y_H
+
+/* This header file contains the functions used to calculate the
+ y-derivatives of displacement.*/
+
+
+
+
+/********************************************************/
+/*** Here we put the actual functions from Table 8. ***/
+
+double df1Asdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*xi*Y11*sindel + dbar*X11/2.0 + alpha/2.0*xi*F);
+}
+
+
+double df1Bsdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-xi*F - dbar*X11 + (1-alpha)/alpha*(xi*Y11+J4)*sindel);
+}
+
+
+double df1Csdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*xi*P*cosdel - alpha*xi*Q);
+}
+
+
+
+
+
+double df2Asdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(alpha/2.0*E);
+}
+
+
+double df2Bsdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-E + (1-alpha)/alpha*(1/R+J5)*sindel);
+}
+
+
+double df2Csdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(2*(1-alpha)*(dbar/(Rsq*R)-Y0*sindel)*sindel - ybar/(R*Rsq)*cosdel - alpha*((cbar+dbar)/(R*Rsq)*sindel - eta/(R*Rsq) - 3*cbar*ybar*q/(Rsq*Rsq*R)));
+}
+
+
+
+
+double df3Asdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*(cosdel/R+q*Y11*sindel) - alpha/2.0*q*F);
+}
+
+
+double df3Bsdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*F - (1-alpha)/alpha*(q*Y11-J6)*sindel);
+}
+
+
+double df3Csdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*q/(R*Rsq) + (ybar/(R*Rsq)-Y0*cosdel)*sindel + alpha*((cbar+dbar)/(R*Rsq)*cosdel + 3*cbar*dbar*q/(Rsq*Rsq*R) - (Y0*cosdel+q*Z0)*sindel));
+}
+
+
+
+
+
+
+
+
+
+double df1Addy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(alpha/2*E);
+}
+
+
+double df1Bddy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-E + (1-alpha)/alpha*J1*sindel*cosdel);
+}
+
+
+double df1Cddy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*eta/(R*Rsq) + Y0*sindel*sindel - alpha*((cbar+dbar)/(R*Rsq)*sindel-3*cbar*ybar*q/(Rsq*Rsq*R)));
+}
+
+
+
+
+
+double df2Addy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*dbar*X11 + xi/2.0*Y11*sindel + alpha/2.0*eta*G);
+}
+
+
+double df2Bddy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-eta*G - xi*Y11*sindel + (1-alpha)/alpha*J2*sindel*cosdel);
+}
+
+
+double df2Cddy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*(X11-ybar*ybar*X32) - alpha*cbar*((dbar+2*q*cosdel)*X32-ybar*eta*q*X53));
+}
+
+
+
+
+double df3Addy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*ybar*X11 - alpha/2.0*q*G);
+}
+
+
+double df3Bddy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*G + (1-alpha)/alpha*J3*sindel*cosdel);
+}
+
+
+double df3Cddy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(xi*P*sindel + ybar*dbar*X32 + alpha*cbar*((ybar+2*q*sindel)*X32-ybar*q*q*X53));
+}
+
+
+
+
+
+
+
+
+
+double df1Atdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/2.0*(cosdel/R+q*Y11*sindel) - alpha/2.0*q*F);
+}
+
+
+double df1Btdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*F - (1-alpha)/alpha*J1*sindel*sindel);
+}
+
+
+double df1Ctdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*(q/(R*Rsq)+Y0*sindel*cosdel) + alpha*(z/(R*Rsq)*cosdel + 3*cbar*dbar*q/(Rsq*Rsq*R) - q*Z0*sindel));
+}
+
+
+
+
+
+double df2Atdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/2.0*ybar*X11 - alpha/2.0*q*G);
+}
+
+
+double df2Btdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*G - (1-alpha)/alpha*J2*sindel*sindel);
+}
+
+
+double df2Ctdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*2*xi*P*sindel - ybar*dbar*X32 + alpha*cbar*((ybar+2*q*sindel)*X32 - ybar*q*q*X53));
+}
+
+
+
+
+double df3Atdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*(dbar*X11+xi*Y11*sindel) + alpha/2.0*q*H);
+}
+
+
+double df3Btdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-q*H - (1-alpha)/alpha*J3*sindel*sindel);
+}
+
+
+double df3Ctdy(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)*(xi*P*cosdel - X11 + ybar*ybar*X32) + alpha*cbar*((dbar+2*q*cosdel)*X32 - ybar*eta*q*X53) + alpha*xi*Q);
+}
+
+
+
+
+
+
+
+/***********************************************************/
+/* y strain-component functions. (At the top of Table 8) */
+
+double duxdy(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ chinnery(df1Asdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df1Asdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Bsdy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Csdy,x,y,z,c,L,W,cosdel,sindel))
+
+ + Ud*(
+ chinnery(df1Addy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df1Addy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Bddy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Cddy,x,y,z,c,L,W,cosdel,sindel))
+
+ + Ut*(
+ chinnery(df1Atdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df1Atdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Btdy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Ctdy,x,y,z,c,L,W,cosdel,sindel));
+
+ return(ans/TWOPI);
+}
+
+double duydy(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ cosdel*( chinnery(df2Asdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Asdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bsdy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Csdy,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Asdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Asdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bsdy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Csdy,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ cosdel*( chinnery(df2Addy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Addy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bddy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Cddy,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Addy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Addy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bddy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Cddy,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ cosdel*( chinnery(df2Atdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Atdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Btdy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Ctdy,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Atdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Atdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Btdy,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Ctdy,x,y,z,c,L,W,cosdel,sindel)));
+ return(ans/TWOPI);
+}
+
+
+
+
+double duzdy(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ sindel*( chinnery(df2Asdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Asdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bsdy,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Csdy,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Asdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Asdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bsdy,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Csdy,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ sindel*( chinnery(df2Addy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Addy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bddy,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Cddy,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Addy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Addy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bddy,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Cddy,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ sindel*( chinnery(df2Atdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df2Atdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Btdy,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Ctdy,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Atdy,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(df3Atdy,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Btdy,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Ctdy,x,y,z,c,L,W,cosdel,sindel)));
+
+ return(ans/TWOPI);
+}
+
+
+
+#endif
Added: cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.z.h
===================================================================
--- cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.z.h 2007-05-03 23:30:26 UTC (rev 6773)
+++ cs/cigma/trunk/sandbox/okada/mcginnis/okada.strain.z.h 2007-05-03 23:33:53 UTC (rev 6774)
@@ -0,0 +1,447 @@
+#ifndef OKADA_STRAIN_Z_H
+#define OKADA_STRAIN_Z_H
+
+/* This header file contains the functions used to calculate the
+ z-derivatives of displacement.*/
+
+
+
+
+/********************************************************/
+/*** Here we put the actual functions from Table 9. ***/
+
+double df1Asdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*xi*Y11*cosdel + ybar/2.0*X11 + alpha/2.0*xi*Fprime);
+}
+
+
+double df1Bsdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-xi*Fprime - ybar*X11 + (1-alpha)/alpha*K1*sindel);
+}
+
+
+double df1Csdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*xi*Pprime*cosdel - alpha*xi*Qprime);
+}
+
+
+
+
+
+double df2Asdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(alpha/2.0*Eprime);
+}
+
+
+double df2Bsdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-Eprime + (1-alpha)/alpha*ybar*D11*sindel);
+}
+
+
+double df2Csdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(2*(1-alpha)*(ybar/(R*Rsq)-Y0*cosdel)*sindel + dbar/(R*Rsq)*cosdel - alpha*((cbar+dbar)/(R*Rsq)*cosdel + 3*cbar*dbar*q/(Rsq*Rsq*R)));
+}
+
+
+
+
+double df3Asdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/2.0*(sindel/R - q*Y11*cosdel) - alpha/2.0*q*Fprime);
+}
+
+
+double df3Bsdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*Fprime + (1-alpha)/alpha*K2*sindel);
+}
+
+
+double df3Csdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((ybar/(R*Rsq)-Y0*cosdel)*cosdel - alpha*((cbar+dbar)/(R*Rsq)*sindel - 3*cbar*ybar*q/(Rsq*Rsq*R) - Y0*sindel*sindel + q*Z0*cosdel));
+}
+
+
+
+
+
+
+
+
+
+double df1Addz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(alpha/2.0*Eprime);
+}
+
+
+double df1Bddz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-Eprime - (1-alpha)/alpha*K3*sindel*cosdel);
+}
+
+
+double df1Cddz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-q/(R*Rsq) + Y0*sindel*cosdel - alpha*((cbar+dbar)/(R*Rsq)*cosdel + 3*cbar*dbar*q/(R*Rsq*Rsq)));
+}
+
+
+
+
+
+double df2Addz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*ybar*X11 + xi/2.0*Y11*cosdel + alpha/2.0*eta*Gprime);
+}
+
+
+double df2Bddz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-eta*Gprime - xi*Y11*cosdel - (1-alpha)/alpha*xi*D11*sindel*cosdel);
+}
+
+
+double df2Cddz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*ybar*dbar*X32 - alpha*cbar*((ybar-2*q*sindel)*X32 + dbar*eta*q*X53));
+}
+
+
+
+
+double df3Addz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-(1-alpha)/2.0*dbar*X11 - alpha/2.0*q*Gprime);
+}
+
+
+double df3Bddz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*Gprime - (1-alpha)/alpha*K4*sindel*cosdel);
+}
+
+
+double df3Cddz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-xi*Pprime*sindel + X11 - dbar*dbar*X32 - alpha*cbar*((dbar-2*q*cosdel)*X32 - dbar*q*q*X53));
+}
+
+
+
+
+
+
+
+
+
+double df1Atdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*(sindel/R - q*Y11*cosdel) - alpha/2.0*q*Fprime);
+}
+
+
+double df1Btdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*Fprime + (1-alpha)/alpha*K3*sindel*sindel);
+}
+
+
+double df1Ctdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-eta/(R*Rsq) + Y0*cosdel*cosdel - alpha*(z/(R*Rsq)*sindel - 3*cbar*ybar*q/(Rsq*Rsq*R) - Y0*sindel*sindel + q*Z0*cosdel));
+}
+
+
+
+
+
+double df2Atdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*dbar*X11 - alpha/2.0*q*Gprime);
+}
+
+
+double df2Btdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(q*Gprime + (1-alpha)/alpha*xi*D11*sindel*sindel);
+}
+
+
+double df2Ctdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*2*xi*Pprime*sindel - X11 + dbar*dbar*X32 - alpha*cbar*((dbar-2*q*cosdel)*X32 - dbar*q*q*X53));
+}
+
+
+
+
+double df3Atdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)/2.0*(ybar*X11 + xi*Y11*cosdel) + alpha/2.0*q*Hprime);
+}
+
+
+double df3Btdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return(-q*Hprime +(1-alpha)/alpha*K4*sindel*sindel);
+}
+
+
+double df3Ctdz(double xi, double eta, double y, double z, double c, double cosdel, double sindel)
+{
+ double Rsq, R;
+ Rsq = (xi*xi + eta*eta + q*q);
+ R = sqrt(Rsq);
+
+ return((1-alpha)*(xi*Pprime*cosdel + ybar*dbar*X32) + alpha*cbar*((ybar-2*q*sindel)*X32 + dbar*eta*q*X53) + alpha*xi*Qprime);
+}
+
+
+
+
+
+
+
+/***********************************************************/
+/* z strain-component functions. (At the top of Table 9) */
+
+double duxdz(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ chinnery(df1Asdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df1Asdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Bsdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f1Cs,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Csdz,x,y,z,c,L,W,cosdel,sindel))
+
+ + Ud*(
+ chinnery(df1Addz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df1Addz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Bddz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f1Cd,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Cddz,x,y,z,c,L,W,cosdel,sindel))
+
+ + Ut*(
+ chinnery(df1Atdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df1Atdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df1Btdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f1Ct,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df1Ctdz,x,y,z,c,L,W,cosdel,sindel));
+
+ return(ans/TWOPI);
+}
+
+double duydz(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ cosdel*( chinnery(df2Asdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df2Asdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bsdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f2Cs,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Csdz,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Asdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df3Asdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bsdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f3Cs,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Csdz,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ cosdel*( chinnery(df2Addz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df2Addz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bddz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f2Cd,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Cddz,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Addz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df3Addz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bddz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f3Cd,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Cddz,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ cosdel*( chinnery(df2Atdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df2Atdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Btdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f2Ct,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df2Ctdz,x,y,z,c,L,W,cosdel,sindel))
+ -sindel*( chinnery(df3Atdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df3Atdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Btdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(f3Ct,x,y,z,c,L,W,cosdel,sindel)
+ + z * chinnery(df3Ctdz,x,y,z,c,L,W,cosdel,sindel)));
+
+ return(ans/TWOPI);
+}
+
+
+
+
+double duzdz(double x, double y, double z,
+ double Us, double Ud, double Ut,
+ double c, double delta, double L, double W)
+{
+ double cosdel, sindel, ans;
+ cosdel = cos(delta);
+ sindel = sin(delta);
+
+ ans = Us*(
+ sindel*( chinnery(df2Asdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df2Asdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bsdz,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2Cs,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Csdz,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Asdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df3Asdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bsdz,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3Cs,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Csdz,x,y,z,c,L,W,cosdel,sindel)))
+ + Ud*(
+ sindel*( chinnery(df2Addz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df2Addz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Bddz,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2Cd,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Cddz,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Addz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df3Addz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Bddz,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3Cd,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Cddz,x,y,z,c,L,W,cosdel,sindel)))
+ + Ut*(
+ sindel*( chinnery(df2Atdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df2Atdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df2Btdz,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f2Ct,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df2Ctdz,x,y,z,c,L,W,cosdel,sindel))
+ +cosdel*( chinnery(df3Atdz,x,y,z,c,L,W,cosdel,sindel)
+ + chinnery(df3Atdz,x,y,-z,c,L,W,cosdel,sindel)
+ + chinnery(df3Btdz,x,y,z,c,L,W,cosdel,sindel)
+ - chinnery(f3Ct,x,y,z,c,L,W,cosdel,sindel)
+ - z * chinnery(df3Ctdz,x,y,z,c,L,W,cosdel,sindel)));
+
+ return(ans/TWOPI);
+}
+
+
+
+#endif
More information about the cig-commits
mailing list