@@ -0,0 +1,1940 @@
+#!/usr/bin/env python%0A%0Aimport sys%0A%0Aimport numpy%0Aimport matplotlib.pyplot as plt%0A%0Aimport clawpack.geoclaw.okada2 as okada%0Aimport clawpack.geoclaw.topo as topo%0Aimport clawpack.geoclaw.dtopotools as dtopotools%0Aimport clawpack.visclaw.colormaps as colormaps%0A%0Aplot_fault = False%0Aif len(sys.argv) %3E 1:%0A    if sys.argv%5B1%5D == %22plot%22:%0A        plot_fault = True%0A%0A# Base subfault%0A# M 7.8%0A# L 90%0A# W 70%0A# depth of hypocenter  20 km (it is written in the text)%0A# coordinates of epicenter: 17.110 -99.100 (taken from: http://usuarios.geofisica.unam.mx/vladimir/sismos/100a%25F1os.html; here they give a hypocentral depth of 33, but let's try with 20, as mentioned aove)%0A# strike 296%0A# dip 25%0A# rake 90%0A# top of the fault plane at 12 km.%0A# They use a value of rigidity of 5 X 10%5E11 dyna/cm%5E2 -%3E 5e11 * 1e2 = 5e13%0Asubfaults = %5B%5D%0Atest_slips = %5B165, 200%5D%0Afor (n, slip) in enumerate(test_slips):%0A    subfaults.append(topo.SubFault(units=%7B%22slip%22:%22cm%22, %22dimensions%22:%22km%22, %22depth%22:%22km%22%7D))%0A    subfaults%5B-1%5D.coordinates = %5B-99.25, 16.6%5D%0A    subfaults%5B-1%5D.coordinate_specification = %22top center%22%0A    subfaults%5B-1%5D.slip = slip%0A    subfaults%5B-1%5D.rake = 90.0%0A    subfaults%5B-1%5D.strike = 296%0A    subfaults%5B-1%5D.dip = 25.0%0A    subfaults%5B-1%5D.depth = 12.0%0A    subfaults%5B-1%5D.dimensions = (90.0, 70.0)%0A    subfaults%5B-1%5D.write('./okada_1957_du%25s.tt3' %25 slip)%0A    %0A    print %22Subfault Characteristics:%22%0A    print %22  Mw = %25s%22 %25 str(subfaults%5B-1%5D.Mw(mu=5e11))%0A    print %22  du = %25s%22 %25 subfaults%5B-1%5D.slip%0A    print %22  Containing Rect = %25s%22 %25 subfaults%5B-1%5D.containing_rect()%0A%0Aif plot_fault:%0A    fig = plt.figure()%0A    for (n,subfault) in enumerate(subfaults):%0A        axes = fig.add_subplot(1,2,n+1)%0A        subfault.plot(axes)%0A        subfault.plot_fault_rect(axes, color='k')%0A        subfault.plot_rake(axes)%0A        axes.set_title(r%22$M_w = %25s$, $%5CDelta u = %25s$ cm%22 %25 (subfault.Mw(mu=5e11), %0A                                                            subfault.slip))%0A    plt.show()

Add simple script for creating Okada deformations