first commit

README.txt

@@ -0,0 +1,20 @@
+These Matlab codes compute strain rates on a triangle mesh using a geodetic velocity field. The method is based on the manuscript K.Johnson (in prep). Body-forces in a thin elastic plate are used to compute a velocity and strain rate field computed at the centroids of the triangles in the mesh.  The inversion solves for the distribution of body forces (at the triangle nodes) that best-fits the geodetic velocity field and slip rates across creeping faults (if included).  You can specify a range of smoothing values (called beta) that minimize the magnitude of the body forces. 
+
+An example fault creep rate file and geodetic data set is included for testing purposes. 
+
+You need to run the scripts in the following order:
+1. setup_mesh.m
+2. build_bodyforce_Greens.m
+3. invert_bodyforce.m
+4. plot_inversion_results.m
+
+Note:  outputs from 1 and 2 and 3 are saved in mat files, so you do not need to rerun all steps every time
+
+1. setup_mesh.m creates the triangular mesh using mesh2d. You need to specify a number of parameters in the Input Section at the top of the script.  Here you specify the origin of the Cartesian coordinate system, load the GPS data file, the creeping fault data file (if you have one), and the nominal node spacing and mesh domain boundaries. You can also experiment with different levels of mesh refinement by change the value of 'refine_mesh'.
+
+2. build_bodyforce_Greens.m computes the Greens Functions which give strains and velocities at centroids of triangles due to unit body force vectors. The script also computes fault creep Greens Functions if you have specified creeping fault segments in the setup_mesh.m script.  The only input you need to specify in this script is the name of the mat file generated in step 1. 
+
+3.  invert_bodyforce.m computes the inversion of your data for the distribution of body forces (and fault creep if that applies).  You can specify a single value or a range of values for smoothing parameter, beta.  The inversion loops over all beta values and computes 'num' realizations of the strain rate and velocity fields for each beta.  You can specify the beta values, num, and the relative weight placed on fitting the creep rate data for creeping faults (if this applies).  
+
+4. plot_inversion_results.m plots the result of the inversion. There are no inputs required to run this script as long as all the parameters and results are currently loaded into the Matlab workspace.
+

build_bodyforce_Greens.m

@@ -0,0 +1,42 @@
+clear all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% 
+
+%load the mat file generated in setup_mesh
+load mesh_wus_creep
+
+%filename for saving Greens functions in mat file
+saveGreens = 'test_Greens'; 
+
+%% No need to modify below
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%assign data to triangles
+find_triangles_data
+
+%take median of data in a triangle and assign this value to the centroid
+assign_data_centroids
+
+
+%creep Greens functions (piece-wise linear slip gradients)
+refine = 1;  %creep discretization (km)
+
+if ~isempty(PatchEnds)
+    lengths = sqrt((PatchEnds(:,3)-PatchEnds(:,1)).^2 + (PatchEnds(:,4)-PatchEnds(:,2)).^2);
+    angle = atan2(PatchEnds(:,4)-PatchEnds(:,2),PatchEnds(:,3)-PatchEnds(:,1));
+    strike=90-angle*180/pi;
+    centers=[(PatchEnds(:,1)+PatchEnds(:,3))./2 (PatchEnds(:,2)+PatchEnds(:,4))./2];
+    pm = [lengths 10^6*ones(length(lengths),1) 10^6*ones(length(lengths),1) 90*ones(length(lengths),1) strike centers(:,1) centers(:,2)];
+    [G1east_creep,G2east_creep,G1north_creep,G2north_creep,...
+        G1Exx_creep,G2Exx_creep,G1Exy_creep,G2Exy_creep,...
+        G1Eyy_creep,G2Eyy_creep] = make_dispG_piecewise(pm,tri_centroids(:,1:2),refine);
+end
+
+nu = .25;  %Poisson's ratio
+[Ge_x,Ge_y,Gn_x,Gn_y,GExx_x,GExy_x,GEyy_x,GExx_y,GExy_y,GEyy_y] = buildG_PointForce(tri_centroids,nodes,nu);
+
+fprintf('\n \n %s \n \n', 'Saving Greens Functions')
+
+save(saveGreens)

creeping_faults.txt

@@ -0,0 +1,183 @@
+% fault id, lon endpoint 1, lat endpoint 1, lon endpoint 2, lat endpoint 2, creep rate (mm/yr)
+32,-120.3,35.752,-120.56,36.003,15
+49,-118.92,34.824,-118.87,34.828,0.1
+49,-118.87,34.828,-118.82,34.849,0.1
+49,-118.82,34.849,-118.77,34.882,0.4
+49,-118.77,34.882,-118.67,34.931,1.8
+49,-118.67,34.931,-118.47,34.996,1.8
+49,-118.47,34.996,-118.38,35.048,1
+49,-118.38,35.048,-118.33,35.079,0.9
+49,-118.33,35.079,-118.12,35.186,0.7
+49,-118.12,35.186,-118.06,35.223,0.4
+49,-118.06,35.223,-118.01,35.271,0.2
+97,-115.17,32.47,-115.5,32.836,9
+97,-115.5,32.836,-115.54,32.88,9
+97,-115.54,32.88,-115.55,32.926,9
+98,-115.89,33.034,-115.8,33.008,3.3
+98,-115.8,33.008,-115.66,32.9,2
+98,-115.66,32.9,-115.61,32.83,1
+99,-116.2,33.207,-116.16,33.182,1.4
+99,-116.16,33.182,-116.13,33.143,2
+99,-116.13,33.143,-116.06,33.098,4
+99,-116.06,33.098,-116.01,33.03,5.5
+99,-116.01,33.03,-115.94,33.001,5.5
+101,-116.28,33.264,-116.2,33.207,0.6
+103,-116.01,32.779,-116.1,32.814,1
+103,-116.1,32.814,-116.18,32.838,1
+103,-116.18,32.838,-116.23,32.889,1
+103,-116.23,32.889,-116.26,32.931,1
+103,-116.26,32.931,-116.34,32.975,1
+104,-115.1,32.161,-115.2,32.201,2
+104,-115.2,32.201,-115.24,32.222,2
+104,-115.24,32.222,-115.29,32.243,2
+104,-115.29,32.243,-115.33,32.263,2
+104,-115.33,32.263,-115.4,32.288,2
+104,-115.4,32.288,-115.72,32.545,2
+104,-115.72,32.545,-115.81,32.592,2
+104,-115.81,32.592,-115.84,32.663,2
+104,-115.84,32.663,-115.88,32.727,2
+119,-117.5,34.263,-117.39,34.192,1
+119,-117.39,34.192,-117.33,34.101,2.5
+119,-117.33,34.101,-117.24,34.017,1
+172,-115.61,32.651,-115.38,32.467,5
+172,-115.38,32.467,-115.34,32.44,5
+172,-115.34,32.44,-115.29,32.405,5
+172,-115.29,32.405,-115.24,32.372,5
+172,-115.24,32.372,-115.16,32.314,5
+172,-115.16,32.314,-115.08,32.262,5
+172,-115.08,32.262,-115.01,32.209,5
+172,-115.01,32.209,-114.93,32.152,5
+283,-117.22,34.15,-117.07,34.093,0.2
+283,-117.07,34.093,-117.01,34.074,0.5
+283,-117.01,34.074,-116.9,34.034,0.6
+283,-116.9,34.034,-116.87,34.011,0.7
+283,-116.87,34.011,-116.82,33.959,0.8
+284,-116.43,33.848,-116.52,33.885,0.7
+284,-116.52,33.885,-116.58,33.907,0.75
+284,-116.58,33.907,-116.62,33.918,0.8
+284,-116.62,33.918,-116.69,33.944,0.85
+284,-116.69,33.944,-116.78,33.937,0.85
+284,-116.78,33.937,-116.8,33.953,0.9
+285,-120.09,35.533,-120.3,35.752,2
+287,-119.21,34.864,-119.03,34.829,0.5
+289,-117.24,34.017,-117.22,34.007,0.1
+289,-117.22,34.007,-117.09,33.905,0.1
+293,-116.91,33.736,-116.84,33.697,1.6
+293,-116.84,33.697,-116.55,33.511,0.8
+293,-116.55,33.511,-116.51,33.49,0.2
+294,-117.14,34.134,-117.08,34.117,0.4
+295,-116.25,33.788,-115.71,33.35,2.5
+341,-118.02,35.283,-117.75,35.423,0.1
+341,-117.75,35.423,-117.49,35.497,0.4
+341,-117.49,35.497,-117.16,35.568,0.5
+341,-117.16,35.568,-117.01,35.602,0.2
+401,-117.09,33.905,-116.91,33.736,1.2
+601,-122.04,37.844,-121.96,37.755,3
+601,-121.96,37.755,-121.93,37.693,6
+601,-121.93,37.693,-121.89,37.62,4
+601,-121.89,37.62,-121.83,37.51,3
+601,-121.83,37.51,-121.82,37.457,2
+602,-121.52,37.057,-121.56,37.132,13
+602,-121.56,37.132,-121.61,37.174,9
+602,-121.61,37.174,-121.67,37.252,8
+602,-121.67,37.252,-121.69,37.3,7
+602,-121.69,37.3,-121.77,37.39,5
+602,-121.77,37.39,-121.82,37.457,3
+603,-121.4,36.839,-121.42,36.894,8
+603,-121.42,36.894,-121.46,36.966,8
+603,-121.46,36.966,-121.52,37.056,13
+621,-121.4,36.839,-121.35,36.785,5
+621,-121.35,36.785,-121.3,36.739,0.5
+621,-121.3,36.739,-121.24,36.683,3
+621,-121.24,36.683,-121.13,36.578,8
+621,-121.13,36.578,-121.06,36.522,10
+621,-121.06,36.522,-120.97,36.426,10
+622,-122.09,38.044,-122.06,37.998,3.4
+622,-122.06,37.998,-122.01,37.944,3
+622,-122.01,37.944,-121.99,37.9,1.2
+623,-122.22,38.415,-122.16,38.233,3.6
+623,-122.16,38.233,-122.13,38.154,4.4
+623,-122.13,38.154,-122.11,38.109,4.4
+623,-122.11,38.109,-122.09,38.044,3.8
+636,-121.54,37.507,-121.62,37.622,2.1
+636,-121.62,37.622,-121.66,37.671,2.1
+636,-121.66,37.671,-121.77,37.788,2.1
+636,-121.77,37.788,-121.8,37.827,2.1
+636,-121.8,37.827,-121.88,37.871,2.1
+638,-121.8,37.397,-121.91,37.483,1.6
+638,-121.91,37.483,-122,37.595,4.6
+638,-122,37.595,-122.05,37.64,5.2
+638,-122.05,37.64,-122.09,37.681,4.7
+638,-122.09,37.681,-122.18,37.782,4
+639,-122.18,37.782,-122.22,37.832,3.9
+639,-122.22,37.832,-122.3,37.929,4.3
+639,-122.3,37.929,-122.41,38.057,4.2
+639,-122.41,38.057,-122.44,38.106,2.5
+639,-122.44,38.106,-122.47,38.196,1.4
+640,-122.22,38.415,-122.25,38.466,3.1
+640,-122.25,38.466,-122.24,38.555,2.8
+640,-122.24,38.555,-122.3,38.623,2.6
+640,-122.3,38.623,-122.34,38.658,2.5
+640,-122.34,38.658,-122.45,38.751,2.3
+644,-122.8,38.681,-122.96,38.816,2
+644,-122.96,38.816,-123.08,39.018,5.6
+644,-123.08,39.018,-123.2,39.2,3.5
+644,-123.2,39.2,-123.23,39.251,3
+644,-123.23,39.251,-123.26,39.294,2.2
+644,-123.26,39.294,-123.34,39.395,1.4
+644,-123.34,39.395,-123.38,39.436,2.2
+644,-123.38,39.436,-123.55,39.723,1.5
+651,-122.44,38.167,-122.49,38.207,3.5
+651,-122.49,38.207,-122.54,38.27,3
+651,-122.54,38.27,-122.64,38.363,2.2
+651,-122.64,38.363,-122.67,38.409,1.8
+651,-122.67,38.409,-122.71,38.468,1.9
+651,-122.71,38.468,-122.76,38.523,5
+651,-122.76,38.523,-122.87,38.638,4.4
+651,-122.87,38.638,-123,38.763,1.6
+657,-121.73,36.987,-121.66,36.933,2
+657,-121.66,36.933,-121.62,36.905,2
+657,-121.62,36.905,-121.52,36.835,8
+657,-121.52,36.835,-121.48,36.804,8
+658,-120.56,36.003,-120.67,36.105,25
+658,-120.67,36.105,-120.74,36.163,25
+658,-120.74,36.163,-120.83,36.243,25
+658,-120.83,36.243,-120.89,36.305,22
+658,-120.89,36.305,-121,36.416,24
+658,-121,36.416,-121.09,36.501,24
+658,-121.09,36.501,-121.17,36.574,21
+658,-121.17,36.574,-121.21,36.611,21
+658,-121.21,36.611,-121.31,36.693,16
+658,-121.31,36.693,-121.37,36.739,14
+658,-121.37,36.739,-121.48,36.804,13
+660,-122.42,37.378,-122.37,37.255,0.5
+660,-122.37,37.255,-122.32,37.164,1
+660,-122.32,37.164,-122.27,37.088,0.5
+662,-121.95,37.139,-121.85,37.108,0.35
+662,-121.85,37.108,-121.78,37.08,1.1
+662,-121.78,37.08,-121.68,37.01,1.9
+662,-121.68,37.01,-121.63,36.98,2.7
+662,-121.63,36.98,-121.52,36.913,2.9
+662,-121.52,36.913,-121.41,36.871,2.9
+665,-122.38,38.372,-122.32,38.308,0.1
+668,-122.5,38.997,-122.6,39.09,0.5
+668,-122.6,39.09,-122.65,39.152,0.2
+668,-122.76,39.216,-122.87,39.35,0.2
+668,-122.87,39.35,-122.92,39.398,0.9
+668,-122.92,39.398,-122.96,39.466,3
+668,-122.96,39.466,-123.05,39.583,6
+668,-123.05,39.583,-123.12,39.65,5
+668,-123.12,39.65,-123.22,39.703,3
+668,-123.22,39.703,-123.27,39.748,1.5
+668,-123.27,39.748,-123.29,39.805,0.5
+668,-123.29,39.805,-123.34,39.9,0.3
+668,-123.34,39.9,-123.4,39.963,0.2
+668,-123.4,39.963,-123.44,40.03,0.2
+668,-123.44,40.03,-123.42,40.077,0.1
+668,-123.42,40.077,-123.45,40.171,0.1
+677,-122.36,38.681,-122.39,38.783,2
+677,-122.39,38.783,-122.38,38.836,2
+677,-122.38,38.836,-122.41,38.884,1.8
+677,-122.41,38.884,-122.46,38.941,1.3
+677,-122.46,38.941,-122.5,38.997,0.8
+999,-124.7,40.37,-126.8,40.412,43