[cig-commits] r7505 - in cs/cigma/branches/cigma-0.9/sandbox: .
dlopen numpy okada-disloc3d python python/_cigmalib python/cigma
luis at geodynamics.org
luis at geodynamics.org
Tue Jun 26 09:53:45 PDT 2007
Author: luis
Date: 2007-06-26 09:53:40 -0700 (Tue, 26 Jun 2007)
New Revision: 7505
Added:
cs/cigma/branches/cigma-0.9/sandbox/dlopen/
cs/cigma/branches/cigma-0.9/sandbox/dlopen/Makefile
cs/cigma/branches/cigma-0.9/sandbox/dlopen/cos2.c
cs/cigma/branches/cigma-0.9/sandbox/dlopen/furls
cs/cigma/branches/cigma-0.9/sandbox/dlopen/hello.cpp
cs/cigma/branches/cigma-0.9/sandbox/dlopen/main1.cpp
cs/cigma/branches/cigma-0.9/sandbox/dlopen/main2.cpp
cs/cigma/branches/cigma-0.9/sandbox/dlopen/polygon.hpp
cs/cigma/branches/cigma-0.9/sandbox/dlopen/triangle.cpp
cs/cigma/branches/cigma-0.9/sandbox/numpy/
cs/cigma/branches/cigma-0.9/sandbox/numpy/Makefile
cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.c
cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.f
cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.c
cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.h
cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3_wrap.c
cs/cigma/branches/cigma-0.9/sandbox/numpy/langtangen.py
cs/cigma/branches/cigma-0.9/sandbox/numpy/numpy_macros.h
cs/cigma/branches/cigma-0.9/sandbox/numpy/tests.py
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/Makefile
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/ReverseSlip_ng.m
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/StrikeSlip_ng.m
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/bm5_analytic_README_2.txt
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.c
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.h
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/dc3d.h
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.c
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.h
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d_mod2.f
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.c
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.h
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5points.c
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/mexfunction.f
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.c
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.m
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/points0.h5.gz
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/points0.in.gz
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.c
cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.h
cs/cigma/branches/cigma-0.9/sandbox/python/
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/Makefile.am
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.c
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.h
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.c
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.h
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.c
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.h
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.c
cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.h
cs/cigma/branches/cigma-0.9/sandbox/python/cigma/
cs/cigma/branches/cigma-0.9/sandbox/python/cigma/Makefile.am
cs/cigma/branches/cigma-0.9/sandbox/python/cigma/__init__.py
cs/cigma/branches/cigma-0.9/sandbox/python/setup.py
Log:
darcs patches:
* Added disloc3d solution to sandbox
* Added directory for cigma python modules to the sandbox
* Removed some kwargs from setup.py (adding them back later)
* Added example numpy extension module to the sandbox
* Added dlopen examples to the sandbox
* Added reference urls to setup.py
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/Makefile
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/Makefile 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/Makefile 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,30 @@
+TARGETS = cos2 main1 main2
+
+all: $(TARGETS)
+
+cos2: cos2.c
+ gcc -rdynamic $< -o $@ -ldl
+
+hello.o: hello.cpp
+ g++ -c $< -o $@
+
+hello.so: hello.o
+ g++ -shared $< -o $@
+
+main1: main1.cpp hello.so
+ g++ $^ -o $@ -ldl
+
+triangle.o: triangle.cpp
+ g++ -c $< -o $@
+
+triangle.so: triangle.o
+ g++ -shared $< -o $@
+
+main2: main2.cpp triangle.so
+ g++ $^ -o $@ -ldl
+
+clean:
+ rm -f *~ *.o *.so
+ rm -f $(TARGETS)
+
+.PHONY: clean
\ No newline at end of file
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/cos2.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/cos2.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/cos2.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,31 @@
+/* http://linux.about.com/library/cmd/blcmdl3_dlopen.htm
+ *
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <dlfcn.h>
+
+int main(int argc, char **argv)
+{
+ void *handle;
+ double (*cosine)(double);
+ char *error;
+
+ handle = dlopen("libm.so", RTLD_LAZY);
+ if (!handle)
+ {
+ fprintf(stderr, "%s\n", dlerror());
+ exit(1);
+ }
+
+ cosine = dlsym(handle, "cos");
+ if ((error = dlerror()) != NULL)
+ {
+ fprintf(stderr, "%s\n", error);
+ exit(1);
+ }
+
+ printf("%lf\n", (*cosine)(2.0));
+ dlclose(handle);
+ return 0;
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/furls
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/furls 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/furls 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,6 @@
+# google: dlopen
+http://www.die.net/doc/linux/man/man3/dlopen.3.html
+http://www.isotton.com/howtos/C++-dlopen-mini-HOWTO/C++-dlopen-mini-HOWTO.html
+
+# google: gcc shared
+http://www.adp-gmbh.ch/cpp/gcc/create_lib.html
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/hello.cpp
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/hello.cpp 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/hello.cpp 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,7 @@
+#include <iostream>
+
+extern "C"
+void hello()
+{
+ std::cout << "hello!" << '\n';
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/main1.cpp
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/main1.cpp 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/main1.cpp 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,44 @@
+#include <iostream>
+#include <dlfcn.h>
+
+int main()
+{
+ using std::cout;
+ using std::cerr;
+
+ cout << "C++ dlopen demo\n\n";
+
+ // open the library
+ cout << "Opening hello.so...\n";
+ void *handle = dlopen("./hello.so", RTLD_LAZY);
+ if (!handle)
+ {
+ cerr << "Cannot open library: " << dlerror() << '\n';
+ return 1;
+ }
+
+ // load the symbol
+ cout << "Loading symbol hello...\n";
+ typedef void (*hello_t)();
+
+ // reset errors
+ dlerror();
+ hello_t hello = (hello_t)dlsym(handle, "hello");
+ const char *dlsym_error = dlerror();
+ if (dlsym_error)
+ {
+ cerr << "Cannot load symbol 'hello': " << dlsym_error << '\n';
+ dlclose(handle);
+ return 1;
+ }
+
+ // use it to do the calculation
+ cout << "Calling hello...\n";
+ hello();
+
+ // close the library
+ cout << "Closing library...\n";
+ dlclose(handle);
+
+ return 0;
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/main2.cpp
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/main2.cpp 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/main2.cpp 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,52 @@
+#include "polygon.hpp"
+#include <iostream>
+#include <dlfcn.h>
+
+int main()
+{
+ using std::cout;
+ using std::cerr;
+
+ // load the triangle library
+ void *triangle = dlopen("./triangle.so", RTLD_LAZY);
+ if (!triangle)
+ {
+ cerr << "Cannot load library: " << dlerror() << '\n';
+ return 1;
+ }
+
+ // reset errors
+ dlerror();
+
+ // load the symbols
+ create_t *create_triangle = (create_t *)dlsym(triangle, "create");
+ const char *dlsym_error = dlerror();
+ if (dlsym_error)
+ {
+ cerr << "Cannot load symbol create: " << dlsym_error << '\n';
+ return 1;
+ }
+
+ destroy_t *destroy_triangle = (destroy_t *)dlsym(triangle, "destroy");
+ dlsym_error = dlerror();
+ if (dlsym_error)
+ {
+ cerr << "Cannot load symbol destroy: " << dlsym_error << '\n';
+ return 1;
+ }
+
+ // create an instance of the class
+ polygon *poly = create_triangle();
+
+ // use the class
+ poly->set_side_length(7);
+ cout << "The area is: " << poly->area() << '\n';
+
+ // destroy the class
+ destroy_triangle(poly);
+
+ // unload the triangle library
+ dlclose(triangle);
+
+ return 0;
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/polygon.hpp
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/polygon.hpp 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/polygon.hpp 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,27 @@
+#ifndef POLYGON_HPP
+#define POLYGON_HPP
+
+class polygon
+{
+protected:
+ double _side_length;
+
+public:
+
+ polygon() : _side_length(0) {}
+
+ virtual ~polygon() {}
+
+ void set_side_length(double side_length)
+ {
+ _side_length = side_length;
+ }
+
+ virtual double area() const = 0;
+};
+
+// the types of the class factories
+typedef polygon *create_t();
+typedef void destroy_t(polygon *);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/dlopen/triangle.cpp
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/dlopen/triangle.cpp 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/dlopen/triangle.cpp 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,23 @@
+#include "polygon.hpp"
+#include <cmath>
+
+class triangle : public polygon
+{
+public:
+ virtual double area() const
+ {
+ return _side_length * _side_length * sqrt(3) / 2;
+ }
+};
+
+// the class factories
+
+extern "C" polygon *create()
+{
+ return new triangle;
+}
+
+extern "C" void destroy(polygon *p)
+{
+ delete p;
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/Makefile
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/Makefile 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/Makefile 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,29 @@
+CC = gcc
+CFLAGS = -O3 -g -Wall
+INCLUDES = -I/usr/include/python2.4 -I/usr/lib/python2.4/site-packages/numpy/core/include
+
+all: ext_gridloop.so ext_gridloop_c.so ext_gridloop_f.so
+
+ext_gridloop_f.so: gridloop.f
+ f2py -m ext_gridloop_f -c gridloop.f
+
+gridloop_c.o: gridloop.c
+ $(CC) $(CFLAGS) $(INCLUDES) -c $< -o $@
+
+gridloop3.o: gridloop3.c gridloop3.h
+ $(CC) $(CFLAGS) $(INCLUDES) -c $< -o $@
+
+gridloop3_wrap.o: gridloop3_wrap.c
+ $(CC) $(CFLAGS) $(INCLUDES) -c $< -o $@
+
+
+ext_gridloop_c.so: gridloop_c.o
+ gcc -shared $^ -o $@
+
+ext_gridloop.so: gridloop3.o gridloop3_wrap.o
+ gcc -shared $^ -o $@
+
+clean:
+ rm -f *.so *.o *.pyc
+
+.PHONY: clean
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,213 @@
+#include <Python.h> /* Python as seen from C */
+#include <numpy/arrayobject.h> /* numpy as seen from C */
+
+
+static PyObject *
+gridloop1(PyObject *self, PyObject *args)
+{
+ PyArrayObject *a, *xcoor, *ycoor;
+ PyObject *func1, *arglist, *result;
+ int nx, ny, i, j;
+ double *a_ij, *x_i, *y_j;
+
+ /* arguments: a, xcoor, ycoor, func1 */
+ if (!PyArg_ParseTuple(args, "O!O!O!O:gridloop1",
+ &PyArray_Type, &a,
+ &PyArray_Type, &xcoor,
+ &PyArray_Type, &ycoor,
+ &func1))
+ {
+ return NULL; /* PyArg_ParseTuple has raised an exception */
+ }
+
+ /* alternative parsing without checking the pointer types
+ if (!PyArg_ParseTuple(args, "OOOO", &a, &xcoor, &ycoor, &func1))
+ {
+ return NULL;
+ }
+ // */
+
+ if (a->nd != 2 || a->descr->type_num != PyArray_DOUBLE)
+ {
+ PyErr_Format(PyExc_ValueError,
+ "a array is %d-dimensional or not of type float",
+ a->nd);
+ return NULL;
+ }
+
+ nx = a->dimensions[0];
+ ny = a->dimensions[1];
+
+ if (xcoor->nd != 1 || xcoor->descr->type_num != PyArray_DOUBLE ||
+ xcoor->dimensions[0] != nx)
+ {
+ PyErr_Format(PyExc_ValueError,
+ "xcoor array has wrong dimension (%d), type or length (%d)",
+ xcoor->nd, xcoor->dimensions[0]);
+ return NULL;
+ }
+
+ if (ycoor->nd != 1 || ycoor->descr->type_num != PyArray_DOUBLE ||
+ ycoor->dimensions[0] != ny)
+ {
+ PyErr_Format(PyExc_ValueError,
+ "ycoor array has wrong dimension (%d), type or length (%d)",
+ ycoor->nd, ycoor->dimensions[0]);
+ return NULL;
+ }
+
+ if (!PyCallable_Check(func1))
+ {
+ PyErr_Format(PyExc_TypeError, "func1 is not a callable function");
+ return NULL;
+ }
+
+ for (i = 0; i < nx; i++)
+ {
+ for (j = 0; j < ny; j++)
+ {
+ a_ij = (double *)(a->data + i*a->strides[0] + j*a->strides[1]);
+ x_i = (double *)(xcoor->data + i*xcoor->strides[0]);
+ y_j = (double *)(ycoor->data + j*ycoor->strides[0]);
+ // shorter: xcoor[i], ycoor[j] for one-dim arrays
+
+ /* call Python function pointed to by func1: */
+ arglist = Py_BuildValue("(dd)", *x_i, *y_j);
+ result = PyEval_CallObject(func1, arglist);
+ Py_DECREF(arglist);
+ if (result == NULL)
+ {
+ /* exception in func1 */
+ return NULL;
+ }
+
+ *a_ij = PyFloat_AS_DOUBLE(result);
+
+ Py_DECREF(result);
+
+#ifdef DEBUG
+ printf("a[%d,%d]=func1(%g,%g)=%g\n",i,j,*x_i,*y_j,*a_ij);
+#endif
+ }
+ }
+ return Py_BuildValue(""); /* return None */
+
+ /* alternative (return 0 for success):
+ return Py_BuildValue("i", 0);
+ // */
+}
+
+
+#include "numpy_macros.h"
+
+static PyObject *
+gridloop2(PyObject *self, PyObject *args)
+{
+ PyArrayObject *a, *xcoor, *ycoor;
+ int a_dims[2];
+ PyObject *func1, *arglist, *result;
+ int nx, ny, i, j;
+ double *a_ij, *x_i, *y_j;
+
+ /* arguments: xcoor, ycoor, func1 */
+ if (!PyArg_ParseTuple(args, "O!O!O:gridloop2",
+ &PyArray_Type, &xcoor,
+ &PyArray_Type, &ycoor,
+ &func1))
+ {
+ return NULL; /* PyArg_ParseTuple has raised an exception */
+ }
+
+ nx = xcoor->dimensions[0];
+ ny = ycoor->dimensions[0];
+
+ NDIM_CHECK(xcoor,1); TYPE_CHECK(xcoor, PyArray_DOUBLE);
+ NDIM_CHECK(ycoor,1); TYPE_CHECK(ycoor, PyArray_DOUBLE);
+ CALLABLE_CHECK(func1);
+
+ /* create return array */
+ a_dims[0] = nx;
+ a_dims[1] = ny;
+ a = (PyArrayObject *)PyArray_FromDims(2, a_dims, PyArray_DOUBLE);
+ if (a == NULL)
+ {
+ printf("Creating %dx%d array failed\n", a_dims[0], a_dims[1]);
+ return NULL; /* PyArray_FromDims raised an exception */
+ }
+
+ for (i = 0; i < nx; i++)
+ {
+ for (j = 0; j < ny; j++)
+ {
+ arglist = Py_BuildValue("(dd)", IND1(xcoor,i), IND1(ycoor,j));
+ result = PyEval_CallObject(func1, arglist);
+ Py_DECREF(arglist);
+ if (result == NULL)
+ {
+ return NULL; /* exception in func1 */
+ }
+
+ IND2(a,i,j) = PyFloat_AS_DOUBLE(result);
+
+ Py_DECREF(result);
+
+#ifdef DEBUG
+ printf("a[%d,%d]=func1(%g,%g)=%g\n", i, j,
+ IND1(xcoor,i), IND1(ycoor,j), IND2(a,i,j));
+#endif
+ }
+ }
+
+ return PyArray_Return(a);
+}
+
+
+
+
+/*
+ * doc strings
+ */
+
+static char gridloop1_doc[] = "gridloop1(a, xcoor, ycoor, pyfunc)";
+static char gridloop2_doc[] = "a = gridloop2(xcoor, ycoor, pyfunc)";
+static char module_doc[] = \
+ "module ext_gridloop:\n\
+ gridloop1(a, xcoor, ycoor, pyfunc)\n\
+ a = gridloop2(xcoor, ycoor, pyfunc)";
+
+/*
+ * The method table must always be present - it lists the
+ * functions that should be callable from Python:
+ */
+static PyMethodDef ext_gridloop_methods[] = {
+
+ {"gridloop1", /* name of func when called from Python */
+ gridloop1, /* corresponding C function */
+ METH_VARARGS, /* ordinary (not keyword) arguments */
+ gridloop1_doc}, /* doc string for gridloop1 function */
+
+ {"gridloop2", /* name of func when called from Python */
+ gridloop2, /* corresponding C function */
+ METH_VARARGS, /* ordinary (not keyword) arguments */
+ gridloop2_doc}, /* doc string for gridloop1 function */
+
+ {NULL, NULL} /* required ending of the method table */
+};
+
+PyMODINIT_FUNC
+initext_gridloop_c()
+{
+ /*
+ * Assign the name of the module, the name of
+ * the method table, and optionally a module
+ * doc string:
+ */
+
+ Py_InitModule3("ext_gridloop_c", ext_gridloop_methods, module_doc);
+
+ /* without module doc string:
+ Py_InitModule("ext_gridloop", ext_gridloop_methods);
+ // */
+
+ import_array(); /* required numpy initialization */
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.f
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.f 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop.f 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,20 @@
+ subroutine gridloop2(a, xcoor, ycoor, nx, ny, func1)
+ integer nx, ny
+ real*8 a(0:nx-1,0:ny-1), xcoor(0:nx-1), ycoor(0:ny-1), func1
+ external func1
+Cf2py intent(out) a
+Cf2py intent(in) xcoor
+Cf2py intent(in) ycoor
+Cf2py depend(nx,ny) a
+
+ integer i,j
+ real*8 x,y
+ do j = 0, ny-1
+ y = ycoor(j)
+ do i = 0, nx-1
+ x = xcoor(i)
+ a(i,j) = func1(x,y)
+ end do
+ end do
+ return
+ end
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,15 @@
+#include "gridloop3.h"
+
+void gridloop3(double **a, double *xcoor, double *ycoor,
+ int nx, int ny, Fxy func1)
+{
+ int i,j;
+ for (i = 0; i < nx; i++)
+ {
+ for (j = 0; j < ny; j++)
+ {
+ a[i][j] = func1(xcoor[i], ycoor[j]);
+ }
+ }
+}
+
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,9 @@
+#ifndef __GRIDLOOP3_H__
+#define __GRIDLOOP3_H__
+
+typedef double (*Fxy)(double x, double y);
+
+void gridloop3(double **a, double *xcoor, double *ycoor,
+ int nx, int ny, Fxy func1);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3_wrap.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3_wrap.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/gridloop3_wrap.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,189 @@
+#include <Python.h> /* Python as seen from C */
+#include <numpy/arrayobject.h> /* numpy as seen from C */
+#include <math.h>
+#include <stdio.h> /* for debug output */
+#include "numpy_macros.h"
+#include "gridloop3.h"
+
+PyObject* _pyfunc_ptr = NULL; /* Python function to be called */
+
+double _pycall(double x, double y)
+{
+ PyObject *arglist, *result;
+ double C_result;
+
+ if (_pyfunc_ptr == NULL)
+ {
+ printf("global pointer _pyfunc_ptr in %s is not initialized (NULL).",
+ __FILE__);
+ exit(1);
+ }
+
+ arglist = Py_BuildValue("(dd)", x, y);
+
+ /*
+ * the global variable _pyfunc_ptr points to the function provided
+ * by the calling Python code
+ */
+ result = PyEval_CallObject(_pyfunc_ptr, arglist);
+ Py_DECREF(arglist);
+ if (result == NULL)
+ {
+ /* cannot return NULL... */
+ printf("Error in callback...");
+ exit(1);
+ }
+
+ C_result = PyFloat_AS_DOUBLE(result);
+ Py_DECREF(result);
+
+ return C_result;
+}
+
+static PyObject *gridloop1(PyObject *self, PyObject *args)
+{
+ PyArrayObject *a, *xcoor, *ycoor;
+ PyObject *func1;
+ int nx, ny, i;
+ double **app;
+ double *ap, *xp, *yp;
+
+ /* arguments: a, xcoor, ycoor, func1 */
+ /* parsing without checking the function types: */
+ if (!PyArg_ParseTuple(args, "OOOO", &a, &xcoor, &ycoor, &func1))
+ {
+ return NULL;
+ }
+
+ NDIM_CHECK(a,2);
+ TYPE_CHECK(a, PyArray_DOUBLE);
+ nx = a->dimensions[0];
+ ny = a->dimensions[1];
+
+ NDIM_CHECK(xcoor, 1);
+ DIM_CHECK(xcoor, 0, nx);
+ TYPE_CHECK(xcoor, PyArray_DOUBLE);
+
+ NDIM_CHECK(ycoor, 1);
+ DIM_CHECK(ycoor, 0, ny);
+ TYPE_CHECK(ycoor, PyArray_DOUBLE);
+
+ CALLABLE_CHECK(func1);
+ _pyfunc_ptr = func1; /* store func1 for use in _pycall */
+
+ /* allocate help array for creating a double pointer: */
+ app = (double **) malloc(nx * sizeof(double *));
+ ap = (double *) a->data;
+
+ for (i = 0; i < nx; i++)
+ {
+ app[i] = &(ap[i*ny]);
+ }
+
+ xp = (double *) xcoor->data;
+ yp = (double *) ycoor->data;
+
+ gridloop3(app, xp, yp, nx, ny, _pycall);
+
+ free(app);
+ return Py_BuildValue(""); /* return None */
+}
+
+
+static PyObject *gridloop2(PyObject *self, PyObject *args)
+{
+ PyArrayObject *a, *xcoor, *ycoor;
+ int a_dims[2];
+ PyObject *func1;
+ int nx, ny, i;
+ double **app;
+ double *ap, *xp, *yp;
+
+ /* arguments: xcoor, ycoor, func1 */
+ if (!PyArg_ParseTuple(args, "OOO", &xcoor, &ycoor, &func1))
+ {
+ return NULL;
+ }
+
+ nx = xcoor->dimensions[0];
+ ny = ycoor->dimensions[0];
+
+ NDIM_CHECK(xcoor, 1);
+ TYPE_CHECK(xcoor, PyArray_DOUBLE);
+
+ NDIM_CHECK(ycoor, 1);
+ TYPE_CHECK(ycoor, PyArray_DOUBLE);
+
+ CALLABLE_CHECK(func1);
+ _pyfunc_ptr = func1;
+
+ /* create return array: */
+ a_dims[0] = nx;
+ a_dims[1] = ny;
+ a = (PyArrayObject *)PyArray_FromDims(2, a_dims, PyArray_DOUBLE);
+ if (a == NULL)
+ {
+ printf("creating %dx%d array failed\n", a_dims[0], a_dims[1]);
+ return NULL; /* PyArray_FromDims raises an exception */
+ }
+
+ /* allocate help array for creating a double pointer: */
+ app = (double **)malloc(nx * sizeof(double *));
+ ap = (double *)a->data;
+
+ for (i = 0; i < nx; i++)
+ {
+ app[i] = &(ap[i*ny]);
+ }
+
+ xp = (double *)xcoor->data;
+ yp = (double *)ycoor->data;
+
+ gridloop3(app, xp, yp, nx, ny, _pycall);
+
+ free(app);
+ return PyArray_Return(a);
+}
+
+
+/* doc strings: */
+static char gridloop1_doc[] = "gridloop1(a, xcoor, ycoor, pyfunc)";
+static char gridloop2_doc[] = "a = gridloop2(xcoor, ycoor, pyfunc)";
+static char module_doc[] = \
+ "module ext_gridloop:\n\
+ gridloop1(a, xcoor, ycoor, pyfunc)\n\
+ a = gridloop2(xcoor, ycoor, pyfunc)";
+
+/*
+ * The method table must always be present - it lists the
+ * functions that should be callable from Python:
+ */
+static PyMethodDef ext_gridloop_methods[] = {
+ {"gridloop1", /* name of func when called from Python */
+ gridloop1, /* corresponding C function */
+ METH_VARARGS, /* ordinary (not keyword) arguments */
+ gridloop1_doc}, /* doc string for gridloop1 function */
+
+ {"gridloop2", /* name of func when called from Python */
+ gridloop2, /* corresponding C function */
+ METH_VARARGS, /* ordinary (not keyword) arguments */
+ gridloop2_doc}, /* doc string for gridloop2 function */
+
+ {NULL, NULL} /* required ending of the method table */
+};
+
+void initext_gridloop()
+{
+ /*
+ * Assign the name of the module, the name of the
+ * method table and (optionally) a module doc string:
+ */
+ Py_InitModule3("ext_gridloop", ext_gridloop_methods, module_doc);
+
+ /* without module doc string:
+ Py_InitModule("ext_gridloop", ext_gridloop_methods);
+ // */
+
+ import_array(); /* required numpy initialization */
+}
+
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/langtangen.py
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/langtangen.py 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/langtangen.py 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,99 @@
+import numpy
+import ext_gridloop_c as ext_gridloop
+
+from numpy import sin
+from numpy import zeros, arange
+from numpy import newaxis, size
+
+def sequence(min=0.0, max=None, inc=1.0):
+ """See Langtangen (pp. 150-151)"""
+ return arange(min, max + inc/2.0, inc)
+seq = sequence
+
+class Grid2D(object):
+ """From Langtangen's Python book (pg. 377)"""
+ def __init__(self, xmin=0, xmax=1, dx=0.5,
+ ymin=0, ymax=1, dy=0.5):
+
+ self.xcoor = seq(xmin, xmax, dx)
+ self.ycoor = seq(ymin, ymax, dy)
+
+ self.xcoorv = self.xcoor[:,newaxis]
+ self.ycoorv = self.ycoor[newaxis,:]
+
+ def __call__(self, f):
+ """
+ Treat f either as an expression containing x and y
+ or as a standard Python function f(x,y). Evaluate
+ the formula or function for all grid points and
+ return the corresponding two-dimensional array.
+ """
+ if isinstance(f, basestring):
+ # introduce the names x and y such that
+ # a simple eval(f) will work for the arrays directly
+ x = self.xcoorv
+ y = self.ycoorv
+ a = eval(f)
+ else:
+ a = f(self.xcoorv, self.ycoorv)
+ return a
+
+ def gridloop(self, f):
+ f_is_str = isinstance(f, basestring)
+ if f_is_str:
+ f_compiled = compile(f, '<string>', 'eval')
+ lx,ly = size(self.xcoor), size(self.ycoor)
+ _a = zeros((lx,ly), float)
+ for i in xrange(lx):
+ x = self.xcoor[i]
+ for j in xrange(ly):
+ y = self.ycoor[j]
+ if f_is_str:
+ _a[i,j] = eval(f_compiled)
+ else:
+ _a[i,j] = f(x,y)
+ return _a
+
+class Grid2Deff(Grid2D):
+ """Langtangen (pg. 433)"""
+
+
+ def ext_gridloop1(self, f):
+ """Compute a[i,j] = f(xi,yj) in an external routine."""
+
+ # a is made here, sent to the routine, and then returned
+ lx,ly = size(self.xcoor), size(self.ycoor)
+ a = zeros((lx,ly), float)
+
+ # <Fortran-specific initialization...>
+
+ if isinstance(f, basestring):
+ f_compiled = compile(f, '<string>', 'eval')
+ def f_str(x,y):
+ return eval(f_compiled)
+ func = f_str
+ else:
+ func = f
+
+ ext_gridloop.gridloop1(a, self.xcoor, self.ycoor, func)
+
+ return a
+
+
+ def ext_gridloop2(self, f):
+ """Compute a[i,j] = f(xi,yj) in an external routine."""
+
+ # a is made in the external routine
+ if isinstance(f, basestring):
+ f_compiled = compile(f, '<string>', 'eval')
+ def f_str(x,y):
+ return eval(f_compiled)
+ func = f_str
+ else:
+ func = f
+
+ a = ext_gridloop.gridloop2(self.xcoor, self.ycoor, func)
+
+ return a
+
+
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/numpy_macros.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/numpy_macros.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/numpy_macros.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,57 @@
+#ifndef __NUMPY_MACROS_H__
+#define __NUMPY_MACROS_H__
+
+/*
+ * This file defines some macros for programming with
+ * numpy arrays in C extension modules.
+ */
+
+/* define some macros for array dimension/type check and subscripting */
+
+#define QUOTE(s) # s /* turn s into string "s" */
+
+#define NDIM_CHECK(a, expected_ndim) \
+ if (a->nd != expected_ndim) { \
+ PyErr_Format(PyExc_ValueError, \
+ "%s array is %d-dimensional, but expected to be %d-dimensional",\
+ QUOTE(a), a->nd, expected_ndim); \
+ return NULL; \
+ }
+
+#define DIM_CHECK(a, dim, expected_length) \
+ if (dim > a->nd) { \
+ PyErr_Format(PyExc_ValueError, \
+ "%s array has no %d dimension (max dim is %d)", \
+ QUOTE(a), dim, a->nd); \
+ return NULL; \
+ } \
+ if (a->dimensions[dim] != expected_length) { \
+ PyErr_Format(PyExc_ValueError, \
+ "%s array has wrong %d-dimension=%d (expected %d)", \
+ QUOTE(a), dim, a->dimensions[dim], expected_length); \
+ return NULL; \
+ }
+
+#define TYPE_CHECK(a, tp) \
+ if (a->descr->type_num != tp) { \
+ PyErr_Format(PyExc_TypeError, \
+ "%s array is not of correct type (%d)", \
+ QUOTE(a), tp); \
+ return NULL; \
+ }
+
+#define CALLABLE_CHECK(func) \
+ if (!PyCallable_Check(func)) { \
+ PyErr_Format(PyExc_TypeError, \
+ "%s is not a callable function", \
+ QUOTE(func)); \
+ return NULL; \
+ }
+
+
+#define IND1(a,i) *((double *)(a->data + i*a->strides[0]))
+#define IND2(a,i,j) *((double *)(a->data + i*a->strides[0] + j*a->strides[1]))
+#define IND3(a,i,j,k) *((double *)(a->data + i*a->strides[0] + j*a->strides[1] + k*a->strides[2]))
+
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/numpy/tests.py
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/numpy/tests.py 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/numpy/tests.py 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+
+def f1(x,y):
+ print 'x+2*y =', x + 2*y
+ return x + 2*y
+
+def verify1():
+ """Basic test of the extension module."""
+ from numpy import allclose
+ from langtangen import Grid2Deff
+
+ g = Grid2Deff(dx=0.5, dy=1)
+ f_exact = g(f1) # numpy computation
+
+ f = g.ext_gridloop2(f1)
+ print 'f computed by external gridloop2 function:\n', f
+ if allclose(f, f_exact, atol=1.0e-10, rtol=1.0e-12):
+ print 'f is correct'
+
+if __name__ == '__main__':
+ verify1()
Added: cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/Makefile
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/Makefile 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/Makefile 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/ReverseSlip_ng.m
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/ReverseSlip_ng.m 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/ReverseSlip_ng.m 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/StrikeSlip_ng.m
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/StrikeSlip_ng.m 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/StrikeSlip_ng.m 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/bm5_analytic_README_2.txt
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/bm5_analytic_README_2.txt 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/bm5_analytic_README_2.txt 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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 = 0; 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/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/cervelli.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/dc3d.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/dc3d.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/dc3d.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d_mod2.f
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d_mod2.f 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/disloc3d_mod2.f 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5array.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/h5points.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5points.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/h5points.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/mexfunction.f
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/mexfunction.f 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/mexfunction.f 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.m
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.m 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/okadasoln_ng.m 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/points0.h5.gz
===================================================================
(Binary files differ)
Property changes on: cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/points0.h5.gz
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/points0.in.gz
===================================================================
(Binary files differ)
Property changes on: cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/points0.in.gz
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/okada-disloc3d/txtarray.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -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/branches/cigma-0.9/sandbox/python/_cigmalib/Makefile.am
===================================================================
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,37 @@
+#include <Python.h>
+
+#include "_cigmalibmodule.h"
+#include "exceptions.h"
+#include "bindings.h"
+
+
+char pycigma_module__doc__[] = "";
+
+/* Initialization function for the module (*must* be called initcigmalib) */
+void
+init_cigmalib()
+{
+ PyObject *m, *d;
+
+ /* create the module and add the functions */
+ m = Py_InitModule4("_cigmalib",
+ pycigma_methods,
+ pycigma_module__doc__,
+ 0,
+ PYTHON_API_VERSION);
+
+ /* get its dictionary */
+ d = PyModule_GetDict(m);
+
+ /* check for errors */
+ if (PyErr_Occurred())
+ {
+ Py_FatalError("can't initialize module _cigmalib");
+ }
+
+ /* install the module exceptions */
+ pycigma_runtimeError = PyErr_NewException("_cigmalib.runtime", 0, 0);
+ PyDict_SetItemString(d, "RuntimeException", pycigma_runtimeError);
+
+ return;
+}
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/_cigmalibmodule.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,18 @@
+#ifndef __PYCIGMA_MODULE_H__
+#define __PYCIGMA_MODULE_H__
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+void init_cigmalib();
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,15 @@
+#include <Python.h>
+#include "bindings.h"
+
+/* The method table */
+struct PyMethodDef pycigma_methods[] = {
+
+ /* Dummy entry for testing */
+ {pycigma_copyright__name__,
+ pycigma_copyright,
+ METH_VARARGS,
+ pycigma_copyright__doc__},
+
+ /* Sentinel */
+ {0, 0, 0, 0}
+};
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/bindings.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,7 @@
+#ifndef __PYCIGMA_BINDINGS_H__
+#define __PYCIGMA_BINDINGS_H__
+
+/* the method table */
+extern struct PyMethodDef pycigma_methods[];
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,4 @@
+#include <Python.h>
+
+PyObject *pycigma_runtimeError = 0;
+
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/exceptions.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,6 @@
+#ifndef __PYCIGMA_EXCEPTIONS_H__
+#define __PYCIGMA_EXCEPTIONS_H__
+
+extern PyObject *pycigma_runtimeError;
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.c
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.c 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.c 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,15 @@
+#include <Python.h>
+#include "misc.h"
+
+/* copyright */
+char pycigma_copyright__doc__[] = "";
+char pycigma_copyright__name__[] = "copyright";
+static char pycigma_copyright_note[] = \
+ "cigma python module; Copyright (c) 2007 California Institute of Technology";
+
+PyObject *
+pycigma_copyright(PyObject *self, PyObject *args)
+{
+ return PyBuildValue("s", pycigma_copyright_note);
+}
+
Added: cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.h
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.h 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/_cigmalib/misc.h 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,8 @@
+#ifndef __PYCIGMA_MISC_H__
+#define __PYCIGMA_MISC_H__
+
+extern char pycigma_copyright__name__[];
+extern char pycigma_copyright__doc__[];
+PyObject *pycigma_copyright(PyObject *, PyObject *);
+
+#endif
Added: cs/cigma/branches/cigma-0.9/sandbox/python/cigma/Makefile.am
===================================================================
Added: cs/cigma/branches/cigma-0.9/sandbox/python/cigma/__init__.py
===================================================================
Added: cs/cigma/branches/cigma-0.9/sandbox/python/setup.py
===================================================================
--- cs/cigma/branches/cigma-0.9/sandbox/python/setup.py 2007-06-26 16:51:42 UTC (rev 7504)
+++ cs/cigma/branches/cigma-0.9/sandbox/python/setup.py 2007-06-26 16:53:40 UTC (rev 7505)
@@ -0,0 +1,12 @@
+#!/usr/bin/env python
+# http://www.python.org/community/sigs/current/distutils-sig/doc/
+# http://www.python.org/files/sigs/sc_submission.html
+# http://www.python.org/~jeremy/weblog/030924.html
+
+from distutils.core import setup
+
+setup(name='Cigma',
+ version='0.9',
+ packages=['cigma','_cigmalib']
+ )
+
More information about the cig-commits
mailing list