source_modelling.ccldpy

   1# Sourced from: https://github.com/tristanbuckreis/ccldpy
   2# This package is not included as a dependency but packaged
   3# directly because it cannot be installed via pip.
   4# ccldpy v2.0.0
   5# June 12, 2024
   6# Tristan E. Buckreis (tristanbuckreis@ucla.edu)
   7
   8###########################################################################################################################################
   9
  10import numpy as np
  11import pandas as pd
  12
  13###########################################################################################################################################
  14
  15
  16TECTONIC_MAPPING = {
  17    "Crustal": "crustal",
  18    "Slab": "intraslab",
  19    "Interface": "interface",
  20    "Outer-rise": "intraslab",
  21    "Undetermined": "crustal",
  22}
  23
  24
  25def WellsCoppersmith1994(magnitude, eqType):
  26    """
  27    Magnitude-Scaling relationships defined by Wells & Coppersmith (1994) for all style-of-faultings.
  28
  29    Input Arguments:
  30    - magnitude = moment magnitude (Mw)
  31    - eqType = type of earthquake (tectonic regime)
  32        "crustal" = shallow-crustal in active tectonic regimes
  33        "intraslab" = intraslab-type in subduction regimes
  34        "interface" = interface-type in subduction regimes
  35        "stable" = shallow-crustal in stable-continental regimes
  36
  37    Output Arguments:
  38    - A = rupture area (km^2)
  39    - AR = rupture aspect ratio (L/W) - median constrained to be >= 1.0 when extrapilated to small magnitudes
  40    - L = rupture along-strike length (km)
  41    - W = rupture down-dip width (km)
  42    """
  43    if eqType in ["intraslab", "interface", "stable"]:
  44        raise ValueError(
  45            "eqType: WellsCoppersmith1994 cannot be used for subduction-type or stable-continental events ('intraslab', 'interface', or 'stable')"
  46        )
  47        return (np.nan, np.nan, np.nan, np.nan)
  48
  49    elif eqType == "crustal":
  50        # Rupture Area:
  51        a1 = -3.49
  52        b1 = 0.91
  53        s1 = 0.24
  54        A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
  55        # Rupture Length:
  56        a2 = -2.44
  57        b2 = 0.59
  58        s2 = 0.16
  59        L = 10 ** (a2 + b2 * magnitude + np.random.normal(0, s2))
  60        # Aspect Ratio:
  61        AR = L**2 / A
  62        # Check if the Aspect Ratio < 1.0, if so constrain and derive L & W:
  63        if AR < 1.0:
  64            # AR = 1.0
  65            s_cy08 = 0.16
  66            AR = np.random.normal(1, s_cy08)
  67            L = np.sqrt(A * AR)
  68            W = np.sqrt(A / AR)
  69        else:
  70            W = L / AR
  71        return (A, AR, L, W)
  72
  73
  74def Leonard2014(magnitude, mechanism, eqType):
  75    """
  76    Magnitude-Scaling relationships defined by Leonard (2014) for different tectonic regimes and style-of-faultings.
  77
  78    Input Arguments:
  79    - magnitude = moment magnitude (Mw)
  80    - mechanism = known or preferred style-of-faulting; can be inferred from rake angle (e.g., Ancheta et al. 2013)
  81        "SS" = strike-slip: (-180 < rake < -150) or (-30 < rake < 30) or (150 < rake < 180)
  82        "NM" = normal: -150 < rake < -30
  83        "RV" = reverse: 30 < rake < 150
  84    - eqType = type of earthquake (tectonic regime)
  85        "crustal" = shallow-crustal in active tectonic regimes
  86        "intraslab" = intraslab-type in subduction regimes
  87        "interface" = interface-type in subduction regimes
  88        "stable" = shallow-crustal in stable-continental regimes
  89
  90    Output Arguments:
  91    - A = rupture area (km^2)
  92    - AR = rupture aspect ratio (L/W) - median constrained to be >= 1.0 when extrapilated to small magnitudes
  93    - L = rupture along-strike length (km)
  94    - W = rupture down-dip width (km)
  95    """
  96    if eqType in ["intraslab", "interface"]:
  97        raise ValueError(
  98            "eqType: Leonard2014 cannot be used for subduction-type events ('intraslab' or 'interface')"
  99        )
 100        return (np.nan, np.nan, np.nan, np.nan)
 101
 102    elif eqType == "crustal":
 103        if mechanism == "SS":
 104            # Rupture Area:
 105            a1 = 3.99
 106            b1 = 1.00
 107            s1 = 0.13
 108            A = 10 ** ((magnitude - a1 - np.random.normal(0, s1)) / b1)
 109            # Rupture Length:
 110            a2 = 4.170
 111            b2 = 1.667
 112            s2 = 0.19
 113            L = 10 ** ((magnitude - a2 - np.random.normal(0, s2)) / b2)
 114            if L > 45.0:
 115                a3 = 5.27
 116                b3 = 1.000
 117                L = 10 ** ((magnitude - a3 - np.random.normal(0, s2)) / b3)
 118        elif mechanism in ["NM", "RV"]:
 119            # Rupture Area:
 120            a1 = 4.00
 121            b1 = 1.00
 122            s1 = 0.15
 123            A = 10 ** ((magnitude - a1 - np.random.normal(0, s1)) / b1)
 124            # Rupture Length:
 125            a2 = 4.000
 126            b2 = 2.000
 127            s2 = 0.23
 128            L = 10 ** ((magnitude - a2 - np.random.normal(0, s2)) / b2)
 129            if L > 5.4:
 130                a3 = 4.240
 131                b3 = 1.667
 132                L = 10 ** ((magnitude - a3 - np.random.normal(0, s2)) / b3)
 133        # Aspect Ratio:
 134        AR = L**2 / A
 135        # Check if the Aspect Ratio < 1.0, if so constrain and derive L & W:
 136        if AR < 1.0:
 137            # AR = 1.0
 138            s_cy08 = 0.16
 139            AR = np.random.normal(1, s_cy08)
 140            L = np.sqrt(A * AR)
 141            W = np.sqrt(A / AR)
 142        else:
 143            W = L / AR
 144        return (A, AR, L, W)
 145
 146    elif eqType == "stable":
 147        if mechanism == "SS":
 148            # Rupture Area:
 149            a1 = 4.18
 150            b1 = 1.00
 151            s1 = 0.09
 152            A = 10 ** ((magnitude - a1 - np.random.normal(0, s1)) / b1)
 153            # Rupture Length:
 154            a2 = 4.250
 155            b2 = 1.667
 156            s2 = 0.18
 157            L = 10 ** ((magnitude - a2 - np.random.normal(0, s2)) / b2)
 158            if L > 60.0:
 159                a3 = 5.44
 160                b3 = 1.000
 161                L = 10 ** ((magnitude - a3 - np.random.normal(0, s2)) / b3)
 162        elif mechanism in ["NM", "RV"]:
 163            # Rupture Area:
 164            a1 = 4.19
 165            b1 = 1.00
 166            s1 = 0.10
 167            A = 10 ** ((magnitude - a1 - np.random.normal(0, s1)) / b1)
 168            # Rupture Length:
 169            a2 = 4.320
 170            b2 = 1.667
 171            s2 = 0.19
 172            L = 10 ** ((magnitude - a2 - np.random.normal(0, s2)) / b2)
 173        # Aspect Ratio:
 174        AR = L**2 / A
 175        # Check if the Aspect Ratio < 1.0, if so constrain and derive L & W:
 176        if AR < 1.0:
 177            # AR = 1.0
 178            s_cy08 = 0.16
 179            AR = np.random.normal(1, s_cy08)
 180            L = np.sqrt(A * AR)
 181            W = np.sqrt(A / AR)
 182        else:
 183            W = L / AR
 184        return (A, AR, L, W)
 185
 186
 187def ThingbaijamEtAl2017(magnitude, mechanism, eqType):
 188    """
 189    Magnitude-Scaling relationships defined by Thingbaijam et al. (2017) for different tectonic regimes and style-of-faultings.
 190
 191    Input Arguments:
 192    - magnitude = moment magnitude (Mw)
 193    - mechanism = known or preferred style-of-faulting; can be inferred from rake angle (e.g., Ancheta et al. 2013)
 194        "SS" = strike-slip: (-180 < rake < -150) or (-30 < rake < 30) or (150 < rake < 180)
 195        "NM" = normal: -150 < rake < -30
 196        "RV" = reverse: 30 < rake < 150
 197    - eqType = type of earthquake (tectonic regime)
 198        "crustal" = shallow-crustal in active tectonic regimes
 199        "intraslab" = intraslab-type in subduction regimes
 200        "interface" = interface-type in subduction regimes
 201        "stable" = shallow-crustal in stable-continental regimes
 202
 203    Output Arguments:
 204    - A = rupture area (km^2)
 205    - AR = rupture aspect ratio (L/W) - median constrained to be >= 1.0 when extrapilated to small magnitudes
 206    - L = rupture along-strike length (km)
 207    - W = rupture down-dip width (km)
 208    """
 209    if eqType in ["intraslab", "stable"]:
 210        raise ValueError(
 211            "eqType: ThingbaijamEtAl2017 cannot be used for instra-slab or stable-continental events ('intraslab' or 'stable')"
 212        )
 213        return (np.nan, np.nan, np.nan, np.nan)
 214
 215    elif eqType == "crustal":
 216        if mechanism == "SS":
 217            # Rupture Area:
 218            a1 = -3.486
 219            b1 = 0.942
 220            s1 = 0.184
 221            A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 222            # Rupture Length:
 223            a2 = -2.943
 224            b2 = 0.681
 225            s2 = 0.151
 226            L = 10 ** (a2 + b2 * magnitude + np.random.normal(0, s2))
 227        elif mechanism == "NM":
 228            # Rupture Area:
 229            a1 = -2.551
 230            b1 = 0.808
 231            s1 = 0.181
 232            A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 233            # Rupture Length:
 234            a2 = -1.722
 235            b2 = 0.485
 236            s2 = 0.128
 237            L = 10 ** (a2 + b2 * magnitude + np.random.normal(0, s2))
 238        elif mechanism == "RV":
 239            # Rupture Area:
 240            a1 = -4.362
 241            b1 = 1.049
 242            s1 = 0.121
 243            A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 244            # Rupture Length:
 245            a2 = -2.693
 246            b2 = 0.614
 247            s2 = 0.083
 248            L = 10 ** (a2 + b2 * magnitude + np.random.normal(0, s2))
 249        # Aspect Ratio:
 250        AR = L**2 / A
 251        # Check if the Aspect Ratio < 1.0, if so constrain and derive L & W:
 252        if AR < 1.0:
 253            # AR = 1.0
 254            s_cy08 = 0.16
 255            AR = np.random.normal(1, s_cy08)
 256            L = np.sqrt(A * AR)
 257            W = np.sqrt(A / AR)
 258        else:
 259            W = L / AR
 260        return (A, AR, L, W)
 261
 262    elif eqType == "interface":
 263        # Rupture Area:
 264        a1 = -3.292
 265        b1 = 0.949
 266        s1 = 0.150
 267        A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 268        # Rupture Length:
 269        a2 = -2.412
 270        b2 = 0.583
 271        s2 = 0.107
 272        L = 10 ** (a2 + b2 * magnitude + np.random.normal(0, s2))
 273        # Aspect Ratio:
 274        AR = L**2 / A
 275        # Check if the Aspect Ratio < 1.0, if so constrain and derive L & W:
 276        if AR < 1.0:
 277            # AR = 1.0
 278            s_cy08 = 0.16
 279            AR = np.random.normal(1, s_cy08)
 280            L = np.sqrt(A * AR)
 281            W = np.sqrt(A / AR)
 282        else:
 283            W = L / AR
 284        return (A, AR, L, W)
 285
 286
 287def ChiouYoungs2008(magnitude, mechanism, eqType, model):
 288    """
 289    Magnitude-Scaling aspect ratio relationship defined by Chiou & Youngs (2008) for different style-of-faultings
 290    in combination with other published magnitude-scaling area relationships.
 291
 292    Input Arguments:
 293    - magnitude = moment magnitude (Mw)
 294    - mechanism = known or preferred style-of-faulting; can be inferred from rake angle (e.g., Ancheta et al. 2013)
 295        "SS" = strike-slip: (-180 < rake < -150) or (-30 < rake < 30) or (150 < rake < 180)
 296        "NM" = normal: -150 < rake < -30
 297        "RV" = reverse: 30 < rake < 150
 298    - eqType = type of earthquake (tectonic regime)
 299        "crustal" = shallow-crustal in active tectonic regimes
 300        "intraslab" = intraslab-type in subduction regimes
 301        "interface" = interface-type in subduction regimes
 302        "stable" = shallow-crustal in stable-continental regimes
 303    - model = magntiude-scaling area relationship to use with Chiou & Youngs (2008) aspect ratio relation
 304        "WellsCoppersmith1994" = Wells & Coppersmith (1994)
 305        "Leonard2014" = Leonard (2014)
 306        "ThingbaijamEtAl2017" = Thingbaijam et al. (2017)
 307
 308    Output Arguments:
 309    - A = rupture area (km^2)
 310    - AR = rupture aspect ratio (L/W) - median constrained to be >= 1.0 when extrapilated to small magnitudes
 311    - L = rupture along-strike length (km)
 312    - W = rupture down-dip width (km)
 313    """
 314    if eqType in ["intraslab", "interface", "stable"]:
 315        raise ValueError(
 316            "eqType: ChiouYoungs2008 cannot be used for subduction-type or stable-continental events ('intraslab', 'interface', or 'stable')"
 317        )
 318        return (np.nan, np.nan, np.nan, np.nan)
 319
 320    elif eqType == "crustal":
 321        if model == "WellsCoppersmith1994":
 322            # Rupture Area:
 323            a1 = -3.49
 324            b1 = 0.91
 325            s1 = 0.24
 326            A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 327        elif model == "Leonard2014":
 328            if mechanism == "SS":
 329                # Rupture Area:
 330                a1 = 3.99
 331                b1 = 1.00
 332                s1 = 0.13
 333                A = 10 ** ((magnitude - a1 - np.random.normal(0, s1)) / b1)
 334            elif mechanism in ["NM", "RV"]:
 335                # Rupture Area:
 336                a1 = 4.00
 337                b1 = 1.00
 338                s1 = 0.15
 339                A = 10 ** ((magnitude - a1 - np.random.normal(0, s1)) / b1)
 340        elif model == "ThingbaijamEtAl2017":
 341            if mechanism == "SS":
 342                # Rupture Area:
 343                a1 = -3.486
 344                b1 = 0.942
 345                s1 = 0.184
 346                A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 347            elif mechanism == "NM":
 348                # Rupture Area:
 349                a1 = -2.551
 350                b1 = 0.808
 351                s1 = 0.181
 352                A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 353            elif mechanism == "RV":
 354                # Rupture Area:
 355                a1 = -4.362
 356                b1 = 1.049
 357                s1 = 0.121
 358                A = 10 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 359        # Aspect Ratio:
 360        s = 0.16
 361        if magnitude < 4.0:
 362            # AR = 1.0
 363            AR = np.random.normal(1, s)
 364        else:
 365            FNM = 0.0
 366            FRV = 0.0
 367            if mechanism == "NM":
 368                FNM = 1.0
 369            elif mechanism == "RV":
 370                FRV = 1.0
 371            ass = 0.01752
 372            anm = -0.00472
 373            arv = -0.01099
 374            b = 3.097
 375            AR = 10 ** (
 376                (ass + anm * FNM + arv * FRV) * (magnitude - 4.0) ** b
 377                + np.random.normal(0, s)
 378            )
 379        # Compute Length and Width:
 380        L = np.sqrt(A * AR)
 381        W = np.sqrt(A / AR)
 382        return (A, AR, L, W)
 383
 384
 385def ContrerasEtAl2022(magnitude, eqType):
 386    """
 387    Magnitude-Scaling aspect ratio relationship defined by Contreras et al. (222) for subduction-type earthquakes.
 388
 389    Input Arguments:
 390    - magnitude = moment magnitude (Mw)
 391    - eqType = type of earthquake (tectonic regime)
 392        "crustal" = shallow-crustal in active tectonic regimes
 393        "intraslab" = intraslab-type in subduction regimes
 394        "interface" = interface-type in subduction regimes
 395        "stable" = shallow-crustal in stable-continental regimes
 396
 397    Output Arguments:
 398    - A = rupture area (km^2)
 399    - AR = rupture aspect ratio (L/W) - median constrained to be >= 1.0 when extrapilated to small magnitudes
 400    - L = rupture along-strike length (km)
 401    - W = rupture down-dip width (km)
 402    """
 403    if eqType in ["crustal", "stable"]:
 404        raise ValueError(
 405            "eqType: ContrerasEtAl2022 cannot be used for non-subduction-type events ('crustal' or 'stable')"
 406        )
 407        return (np.nan, np.nan, np.nan, np.nan)
 408
 409    elif eqType == "interface":
 410        # Rupture Area:
 411        a1 = -3.8290
 412        b1 = 1.0
 413        s1 = 0.270
 414        A = 10.0 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 415        # Aspect Ratio:
 416        if magnitude > 7.25:
 417            AR = 10.0 ** (
 418                0.2759 * (magnitude - 7.25) + np.random.normal(0, 0.192)
 419            )  # std = 0.441 in ln unit
 420        else:
 421            AR = 10.00 ** (0.0 + np.random.normal(0, 0.0717))  # std = 0.165 in ln unit
 422        # Compute Length and Width:
 423        L = np.sqrt(A * AR)
 424        W = np.sqrt(A / AR)
 425        return (A, AR, L, W)
 426
 427    elif eqType == "intraslab":
 428        # Rupture Area:
 429        a1 = -3.251
 430        b1 = 0.890
 431        s1 = 0.184
 432        A = 10.0 ** (a1 + b1 * magnitude + np.random.normal(0, s1))
 433        # Aspect Ratio:
 434        if magnitude > 6.5:
 435            AR = 10.0 ** (
 436                0.0938 * (magnitude - 6.5) + np.random.normal(0, 0.164)
 437            )  # std = 0.378in ln unit
 438        else:
 439            AR = 10.00 ** (0.0 + np.random.normal(0, 0.104))  # std = 0.239 in ln unit
 440        # Compute Length and Width:
 441        L = np.sqrt(A * AR)
 442        W = np.sqrt(A / AR)
 443        return (A, AR, L, W)
 444
 445
 446def get_mechanism_based_on_rake(rake):
 447    """
 448    Infer style-of-faulting (mechanism) from rake angle (e.g., Ancheta et al. 2013)
 449
 450    Input Arguments:
 451    rake = rake angle (degrees)
 452
 453    Output Arguments:
 454    mechanism = indicator string for style-of-faulting (mechanism)
 455        "SS" = strike-slip: (-180 < rake < -150) or (-30 < rake < 30) or (150 < rake < 180)
 456        "NM" = normal: -150 < rake < -30
 457        "RV" = reverse: 30 < rake < 150
 458    """
 459    if (-180 <= rake < -150) | (-30 <= rake < 30) | (150 <= rake <= 180):
 460        return "SS"
 461    elif (-120 <= rake < -60) | (-150 <= rake < -120) | (-60 <= rake < -30):
 462        return "NM"
 463    elif (60 <= rake < 120) | (30 <= rake < 60) | (120 <= rake < 150):
 464        return "RV"
 465
 466
 467def discrete(n, x, p, tp):
 468    for i in range(10):
 469        if (tp >= p[i]) and (tp <= p[i + 1]):
 470            tx = x[i] + (x[i + 1] - x[i]) * (tp - p[i]) / (p[i + 1] - p[i])
 471            return tx
 472    return
 473
 474
 475def get_hyp_down_dip_position(eqtype, region):
 476    """
 477    Simulate hypocenter's down-dip relative position on the rupture surface based on Chiou & Youngs (2008)
 478
 479    Input arguments:
 480    - eqType = type of earthquake (tectonic regime)
 481        "crustal" = shallow-crustal in active tectonic regimes
 482        "intraslab" = intraslab-type in subduction regimes
 483        "interface" = interface-type in subduction regimes
 484        "stable" = shallow-crustal in stable-continental regimes
 485    - region = geographic region where the earthquake occured
 486        "japan"
 487        "chile"
 488        "other"
 489
 490    Output Arguments:
 491    - fd = down-dip relative position of hypocenter on rupture surface [0, 1]
 492        0 ~ along the top-depth of the rupture
 493        1 ~ along the bottom-depth of the rupture
 494    """
 495    nxf = 11
 496    xdf = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.00])
 497    ssd = np.array([0, 0.025, 0.05, 0.1, 0.175, 0.275, 0.4, 0.55, 0.7, 0.85, 1.00])
 498    iad = np.array(
 499        [0.00, 0.012, 0.051, 0.139, 0.294, 0.500, 0.706, 0.861, 0.949, 0.988, 1.00]
 500    )
 501    ied0 = np.array(
 502        [0.000, 0.024, 0.085, 0.206, 0.389, 0.599, 0.783, 0.906, 0.969, 0.993, 1.00]
 503    )
 504    ied1 = np.array(
 505        [0.0, 0.002, 0.012, 0.044, 0.121, 0.262, 0.460, 0.671, 0.843, 0.950, 1.00]
 506    )
 507    ied2 = np.array(
 508        [0.00, 0.013, 0.053, 0.143, 0.297, 0.500, 0.703, 0.857, 0.947, 0.987, 1.00]
 509    )
 510    if eqtype in ["crustal", "stable"]:
 511        fd = discrete(nxf, xdf, ssd, np.random.rand())
 512    elif eqtype == "intraslab":
 513        fd = discrete(nxf, xdf, iad, np.random.rand())
 514    elif eqtype == "interface":
 515        if region == "japan":
 516            fd = discrete(nxf, xdf, ied0, np.random.rand())
 517        elif region == "chile":
 518            fd = discrete(nxf, xdf, ied1, np.random.rand())
 519        elif region == "other":
 520            fd = discrete(nxf, xdf, ied2, np.random.rand())
 521    try:
 522        return fd
 523    except:
 524        return 0.0
 525
 526
 527def get_hyp_along_strike_position(eqtype):
 528    """
 529    Simulate hypocenter's along-strike relative position on the rupture surface based on Chiou & Youngs (2008)
 530
 531    Input arguments:
 532    - eqType = type of earthquake (tectonic regime)
 533        "crustal" = shallow-crustal in active tectonic regimes
 534        "intraslab" = intraslab-type in subduction regimes
 535        "interface" = interface-type in subduction regimes
 536        "stable" = shallow-crustal in stable-continental regimes
 537
 538    Output Arguments:
 539    - fl = along-strike relative position of hypocenter on rupture surface [0, 1]
 540        0 ~ along the "left" edge of the rupture
 541        1 ~ along the "right" edge of the rupture
 542    """
 543    nxf = 11
 544    xdf = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.00])
 545    hypx = np.array([0, 0.05, 0.125, 0.225, 0.35, 0.5, 0.65, 0.775, 0.875, 0.95, 1.00])
 546    hypxe = np.array(
 547        [0.000, 0.007, 0.034, 0.112, 0.272, 0.500, 0.728, 0.888, 0.966, 0.993, 1.00]
 548    )
 549    hypxa = np.array(
 550        [0.00, 0.015, 0.057, 0.148, 0.301, 0.500, 0.699, 0.852, 0.943, 0.985, 1.00]
 551    )
 552    if eqtype in ["crustal", "stable"]:
 553        fl = discrete(nxf, xdf, hypx, np.random.rand())
 554    elif eqtype == "intraslab":
 555        fl = discrete(nxf, xdf, hypxa, np.random.rand())
 556    elif eqtype == "interface":
 557        fl = discrete(nxf, xdf, hypxe, np.random.rand())
 558    return fl
 559
 560
 561def get_median_index(array):
 562    index = np.argpartition(array, len(array) // 2)[len(array) // 2]
 563    return index
 564
 565
 566def LL2XY(zlon, zlat, lon, lat):
 567    """
 568    Convert latitude/longitude to x/y coordinates on the basis of a spherical earth
 569    (ignoring ellipsoidal effects) using the haversine formula.
 570
 571    Reference: http://www.movable-type.co.uk/scripts/latlong.html
 572
 573    Input Arguments:
 574    - zlon = centered longitude (degrees): corresponding to x/y coordinate of (0,0)
 575    - zlat = centered latitude (degrees): corresponding to x/y coordinate of (0,0)
 576    - lon = numpy array of longitudes (degrees)
 577    - lat = numpy array of latitudes (degrees)
 578
 579    Output Arguments:
 580    -x = numpy array of x-coordinate
 581    -y = numpy array of y-coordinates
 582
 583    """
 584    Rearth = 6371.0
 585    dtr = np.arcsin(1.0) / 90
 586    lam1 = (zlon + 360) * dtr if zlon < 0 else zlon * dtr
 587    phi1 = zlat * dtr
 588
 589    # lam2 = (lon + 360)*dtr if lon < 0 else lon*dtr
 590    if type(lon) == "float":
 591        lam2 = (lon + 360) * dtr if lon < 0 else lon * dtr
 592    else:
 593        lam2 = np.where(lon < 0, (lon + 360) * dtr, lon * dtr)
 594
 595    phi2 = lat * dtr
 596    dphi = phi2 - phi1
 597    dlam = lam2 - lam1
 598    a = np.sin(dphi / 2) ** 2 + np.cos(phi1) * np.cos(phi2) * np.sin(dlam / 2) ** 2
 599    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
 600    dist = Rearth * c
 601    if dist == 0:
 602        theta = 0
 603    else:
 604        theta = np.arctan2(
 605            np.sin(dlam) * np.cos(phi2),
 606            np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2) * np.cos(dlam),
 607        )
 608    x = dist * np.sin(theta)
 609    y = dist * np.cos(theta)
 610
 611    return (x, y)
 612
 613
 614def XY2LL(zlon, zlat, x, y):
 615    """
 616    Convert x/y coordinates to latitude/longitude on the basis of a spherical earth
 617    (ignoring ellipsoidal effects) using the haversine formula.
 618
 619    Reference: http://www.movable-type.co.uk/scripts/latlong.html
 620
 621    Input Arguments:
 622    - zlon = centered longitude (degrees): corresponding to x/y coordinate of (0,0)
 623    - zlat = centered latitude (degrees): corresponding to x/y coordinate of (0,0)
 624    - x = numpy array of x-coordinate
 625    - y = numpy array of y-coordinates
 626
 627    Output Arguments:
 628    - lon = numpy array of longitudes (degrees)
 629    - lat = numpy array of latitudes (degrees)
 630
 631    """
 632    Rearth = 6371.0
 633    dtr = np.arcsin(1.0) / 90
 634    lam1 = (zlon + 360) * dtr if zlon < 0 else zlon * dtr
 635    phi1 = zlat * dtr
 636
 637    d = np.sqrt(x**2 + y**2)
 638    delta = d / Rearth
 639    theta = np.arctan2(x, y)
 640    phi2 = np.arcsin(
 641        np.sin(phi1) * np.cos(delta) + np.cos(phi1) * np.sin(delta) * np.cos(theta)
 642    )
 643    lam2 = lam1 + np.arctan2(
 644        np.sin(theta) * np.sin(delta) * np.cos(phi1),
 645        np.cos(delta) - np.sin(phi1) * np.sin(phi2),
 646    )
 647    lat = phi2 / dtr
 648    lon = lam2 / dtr
 649    if type(lon) == "float":
 650        lon = lon - 360 if lon > 180 else lon
 651    else:
 652        lon = np.where(lon > 180, lon - 360, lon)
 653
 654    return (lon, lat)
 655
 656
 657def pointTriangleDistance(TRI_xyz, P_xyz):
 658    """
 659    Compute shortest 3D distance between a triangle and point (Eberly 1999)
 660
 661    Reference: https://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
 662
 663    Coded in Python by Meera Kota (UCLA)
 664
 665    Input Arguments:
 666    - TRI_xyz = numpy arrays of 3D triangles in XY space [[x1,y1,z1],[x2,y2,z2],[x3,y3,z3]]
 667    - P_xyz = numpy array of 3D points in XY space [[x1,y1,z1],...,[xn,yn,zn]] where n = number of points
 668    """
 669    E0 = TRI_xyz[:, 1] - TRI_xyz[:, 0]
 670    E1 = TRI_xyz[:, 2] - TRI_xyz[:, 0]
 671    a = np.sum(np.multiply(E0, E0), axis=1)
 672    b = np.sum(np.multiply(E0, E1), axis=1)
 673    c = np.sum(np.multiply(E1, E1), axis=1)
 674
 675    det = a * c - b * b
 676    det[det == 0] = 1.0e-8
 677
 678    D = TRI_xyz[:, 0] - P_xyz
 679
 680    d = np.sum(np.multiply(E0, D), axis=2)
 681    e = np.sum(np.multiply(E1, D), axis=2)
 682    f = np.sum(np.multiply(D, D), axis=2)
 683
 684    s = b * e - c * d
 685    t = b * d - a * e
 686
 687    sqrdistance = np.empty((len(P_xyz), len(TRI_xyz)), dtype=float)
 688
 689    # Region 4
 690    cond = (s + t <= det) & (s < 0.0) & (t < 0.0) & (d < 0.0) & (-d >= a)
 691    sqrdistance[cond] = (a + 2.0 * d + f)[cond]
 692    cond = (s + t <= det) & (s < 0.0) & (t < 0.0) & (d < 0.0) & (-d < a)
 693    sqrdistance[cond] = (-d * d / a + f)[cond]
 694    cond = (s + t <= det) & (s < 0.0) & (t < 0.0) & (d >= 0.0) & (e >= 0.0)
 695    sqrdistance[cond] = (f)[cond]
 696    cond = (s + t <= det) & (s < 0.0) & (t < 0.0) & (d >= 0.0) & (e < 0.0) & (-e >= c)
 697    sqrdistance[cond] = (c + 2.0 * e + f)[cond]
 698    cond = (s + t <= det) & (s < 0.0) & (t < 0.0) & (d >= 0.0) & (e < 0.0) & (-e < c)
 699    sqrdistance[cond] = (-e * e / c + f)[cond]
 700
 701    # Region 3
 702    cond = (s + t <= det) & (s < 0.0) & (t >= 0.0) & (e >= 0.0)
 703    sqrdistance[cond] = (f)[cond]
 704    cond = (s + t <= det) & (s < 0.0) & (t >= 0.0) & (e < 0.0) & (-e >= c)
 705    sqrdistance[cond] = (c + 2.0 * e + f)[cond]
 706    cond = (s + t <= det) & (s < 0.0) & (t >= 0.0) & (e < 0.0) & (-e < c)
 707    sqrdistance[cond] = (-e * e / c + f)[cond]
 708
 709    # Region 5
 710    cond = (s + t <= det) & (s >= 0.0) & (t < 0.0) & (d >= 0.0)
 711    sqrdistance[cond] = (f)[cond]
 712    cond = (s + t <= det) & (s >= 0.0) & (t < 0.0) & (d < 0.0) & (-d >= a)
 713    sqrdistance[cond] = (a + 2.0 * d + f)[cond]
 714    cond = (s + t <= det) & (s >= 0.0) & (t < 0.0) & (d < 0.0) & (-d < a)
 715    sqrdistance[cond] = (-d * d / a + f)[cond]
 716
 717    # Region 0
 718    invDet = 1.0 / det
 719    stemp = s * invDet
 720    ttemp = t * invDet
 721    cond = (s + t <= det) & (s >= 0) & (t >= 0)
 722    sqrdistance[cond] = (
 723        stemp * (a * stemp + b * ttemp + 2.0 * d)
 724        + ttemp * (b * stemp + c * ttemp + 2.0 * e)
 725        + f
 726    )[cond]
 727
 728    # Region 2
 729    tmp0 = b + d
 730    tmp1 = c + e
 731    numer = tmp1 - tmp0
 732    denom = a - 2.0 * b + c
 733    denom[denom == 0] = 1.0e-8
 734    stemp = numer / denom
 735    ttemp = 1.0 - stemp
 736    cond = (s + t > det) & (s < 0.0) & (tmp1 > tmp0) & (numer >= denom)
 737    sqrdistance[cond] = (a + 2.0 * d + f)[cond]
 738    cond = (s + t > det) & (s < 0.0) & (tmp1 > tmp0) & (numer < denom)
 739    sqrdistance[cond] = (
 740        stemp * (a * stemp + b * ttemp + 2.0 * d)
 741        + ttemp * (b * stemp + c * ttemp + 2.0 * e)
 742        + f
 743    )[cond]
 744    cond = (s + t > det) & (s < 0.0) & (tmp1 <= tmp0) & (tmp1 <= 0.0)
 745    sqrdistance[cond] = (c + 2.0 * e + f)[cond]
 746    cond = (s + t > det) & (s < 0.0) & (tmp1 <= tmp0) & (tmp1 > 0.0) & (e >= 0.0)
 747    sqrdistance[cond] = (f)[cond]
 748    cond = (s + t > det) & (s < 0.0) & (tmp1 <= tmp0) & (tmp1 > 0.0) & (e < 0.0)
 749    sqrdistance[cond] = (-e * e / c + f)[cond]
 750
 751    # Region 6
 752    tmp0 = b + e
 753    tmp1 = a + d
 754    numer = tmp1 - tmp0
 755    denom = a - 2.0 * b + c
 756    denom[denom == 0] = 1.0e-8
 757    ttemp = numer / denom
 758    stemp = 1.0 - ttemp
 759    cond = (s + t > det) & (s >= 0) & (t < 0) & (tmp1 > tmp0) & (numer >= denom)
 760    sqrdistance[cond] = (c + 2.0 * e + f)[cond]
 761    cond = (s + t > det) & (s >= 0) & (t < 0) & (tmp1 > tmp0) & (numer < denom)
 762    sqrdistance[cond] = (
 763        stemp * (a * stemp + b * ttemp + 2.0 * d)
 764        + ttemp * (b * stemp + c * ttemp + 2.0 * e)
 765        + f
 766    )[cond]
 767    cond = (s + t > det) & (s >= 0) & (t < 0) & (tmp1 <= tmp0) & (tmp1 <= 0)
 768    sqrdistance[cond] = (a + 2.0 * d + f)[cond]
 769    cond = (s + t > det) & (s >= 0) & (t < 0) & (tmp1 <= tmp0) & (tmp1 > 0) & (d >= 0)
 770    sqrdistance[cond] = (f)[cond]
 771    cond = (s + t > det) & (s >= 0) & (t < 0) & (tmp1 <= tmp0) & (tmp1 > 0) & (d < 0)
 772    sqrdistance[cond] = (-d * d / a + f)[cond]
 773
 774    # Region 1
 775    numer = c + e - b - d
 776    denom = a - 2.0 * b + c
 777    stemp = numer / denom
 778    ttemp = 1.0 - stemp
 779    cond = (s + t > det) & (s >= 0) & (t >= 0) & (numer <= 0)
 780    sqrdistance[cond] = (c + 2.0 * e + f)[cond]
 781    cond = (s + t > det) & (s >= 0) & (t >= 0) & (numer > 0) & (numer >= denom)
 782    sqrdistance[cond] = (a + 2.0 * d + f)[cond]
 783    cond = (s + t > det) & (s >= 0) & (t >= 0) & (numer > 0) & (numer < denom)
 784    sqrdistance[cond] = (
 785        stemp * (a * stemp + b * ttemp + 2.0 * d)
 786        + ttemp * (b * stemp + c * ttemp + 2.0 * e)
 787        + f
 788    )[cond]
 789
 790    # account for numerical round-off error
 791    sqrdistance[sqrdistance <= 0] = 0
 792    return np.sqrt(sqrdistance).T
 793
 794
 795def check_input_arguments(
 796    eqType, method, nsims, strike, dip, rake, strike2, dip2, rake2
 797):
 798    # Check that required arguments are defined for each simulation method  -------------------------------------------------
 799    if method in ["A", "D"]:
 800        if None in np.array([strike, dip, rake]):
 801            raise ValueError(
 802                "'strike', 'dip', and 'rake' must be defined when using 'method' = '%s'"
 803                % method
 804            )
 805    if method == "B":
 806        if None in np.array([strike2, dip2, rake2]):
 807            raise ValueError(
 808                "'strike2', 'dip2', and 'rake2' must be defined when using 'method' = 'B'"
 809            )
 810    if method == "C":
 811        if None in np.array([strike, dip, rake, strike2, dip2, rake2]):
 812            raise ValueError(
 813                "'strike', 'dip', 'rake', 'strike2', 'dip2', and 'rake2' must be defined when using 'method' = 'C'"
 814            )
 815
 816    # Check that number of simulations are compatible with sepcified eqType  ------------------------------------------------
 817    if (eqType == "crustal") & (nsims[6] > 0):
 818        raise ValueError(
 819            "nsims[6] = %i: ContrerasEtAl2022 cannot be used for 'crustal' events"
 820            % nsims[6]
 821        )
 822    elif eqType == "stable":
 823        if nsims[0] > 0:
 824            raise ValueError(
 825                "nsims[0] = %i: WellsCoppersmith1994 cannot be used for 'stable' events"
 826                % nsims[0]
 827            )
 828        elif nsims[2] > 0:
 829            raise ValueError(
 830                "nsims[2] = %i: ThingbaijamEtAl2017 cannot be used for 'stable' events"
 831                % nsims[2]
 832            )
 833        elif sum(nsims[3:6]) > 0:
 834            raise ValueError(
 835                "nsims[3:6] = [%i,%i,%i]: ChiouYoungs2008 cannot be used for 'stable' events"
 836                % (nsims[3], nsims[4], nsims[5])
 837            )
 838        elif nsims[6] > 0:
 839            raise ValueError(
 840                "nsims[6] = %i: ContrerasEtAl2022 cannot be used for 'stable' events"
 841                % nsims[6]
 842            )
 843    elif eqType in ["interface", "intraslab"]:
 844        if nsims[0] > 0:
 845            raise ValueError(
 846                "nsims[0] = %i: WellsCoppersmith1994 cannot be used for '%s' events"
 847                % (nsims[0], eqType)
 848            )
 849        elif nsims[1] > 0:
 850            raise ValueError(
 851                "nsims[1] = %i: Leonard2014 cannot be used for '%s' events"
 852                % (nsims[1], eqType)
 853            )
 854        elif sum(nsims[3:6]) > 0:
 855            raise ValueError(
 856                "nsims[3:6] = [%i,%i,%i]: ChiouYoungs2008 cannot be used for '%s' events"
 857                % (nsims[3], nsims[4], nsims[5], eqType)
 858            )
 859
 860    # Check that there are only positive integers for the number of simulations  --------------------------------------------
 861    if np.any(np.array(nsims), where=np.array(nsims) < 0):
 862        raise ValueError(
 863            "nsims = %s: can not specify a negative number of simulatiuons for any scaling relationship"
 864            % str(list(nsims))
 865        )
 866
 867    # Check that the total number of simulations is odd and greater than zero  ----------------------------------------------
 868    if sum(nsims) % 2 == 0:
 869        if eqType == "crustal":
 870            if nsims[0] > 0:
 871                nsims[0] = nsims[0] + 1
 872            elif nsims[1] > 0:
 873                nsims[1] = nsims[1] + 1
 874            elif nsims[2] > 0:
 875                nsims[2] = nsims[2] + 1
 876            elif nsims[3] > 0:
 877                nsims[3] = nsims[3] + 1
 878            elif nsims[4] > 0:
 879                nsims[4] = nsims[4] + 1
 880            elif nsims[5] > 0:
 881                nsims[5] = nsims[5] + 1
 882            else:
 883                raise ValueError(
 884                    "nsims = %s: must specify the number of simulations partitioned to each scaling relationship"
 885                    % str(list(nsims))
 886                )
 887        if eqType == "stable":
 888            if nsims[1] > 0:
 889                nsims[1] = nsims[1] + 1
 890            else:
 891                raise ValueError(
 892                    "nsims = %s: must specify the number of simulations partitioned to each scaling relationship"
 893                    % str(list(nsims))
 894                )
 895        if eqType in ["interface", "intraslab"]:
 896            if nsims[6] > 0:
 897                nsims[6] = nsims[6] + 1
 898            elif nsims[2] > 0:
 899                nsims[2] = nsims[2] + 1
 900            else:
 901                raise ValueError(
 902                    "nsims = %s: must specify the number of simulations partitioned to each scaling relationship"
 903                    % str(list(nsims))
 904                )
 905
 906
 907###########################################################################################################################################
 908
 909v_WellsCoppersmith1994 = np.vectorize(WellsCoppersmith1994)
 910v_Leonard2014 = np.vectorize(Leonard2014)
 911v_ThingbaijamEtAl2017 = np.vectorize(ThingbaijamEtAl2017)
 912v_ChiouYoungs2008 = np.vectorize(ChiouYoungs2008)
 913v_ContrerasEtAl2022 = np.vectorize(ContrerasEtAl2022)
 914v_get_mechanism_based_on_rake = np.vectorize(get_mechanism_based_on_rake)
 915v_get_hyp_down_dip_position = np.vectorize(get_hyp_down_dip_position)
 916v_get_hyp_along_strike_position = np.vectorize(get_hyp_along_strike_position)
 917v_get_median_index = np.vectorize(get_median_index, signature="(n)->()")
 918v_LL2XY = np.vectorize(LL2XY, signature="(),(),(n),(n)->(n),(n)")
 919v_XY2LL = np.vectorize(XY2LL, signature="(),(),(n),(n)->(n),(n)")
 920
 921###########################################################################################################################################
 922
 923
 924def get_rupture_surface_simulations(
 925    eqid,
 926    eqType,
 927    region,
 928    elat,
 929    elon,
 930    hypd,
 931    magnitude,
 932    method,
 933    nsims,
 934    mechanism=None,
 935    strike=None,
 936    dip=None,
 937    rake=None,
 938    strike2=None,
 939    dip2=None,
 940    rake2=None,
 941):
 942    """
 943    Simulate earthquake rupture surface that minimizes the difference between the median distance of a pseudo-grid of sites
 944    and a stochastic set of possible ruptures. Randomized fault attributes are dependent on the fault category. Fault scaling
 945    models depend on earthquake type. Adapted from the CCLD routine originally coded in Frortran by Chiou and Youngs (2008).
 946
 947    Input Arguments:
 948    - eqid = unique integer identifier for the event (to set the random seed)
 949    - eqType = type of earthquake:
 950        "crustal" = shallow-crustal in active tectonic regimes
 951        "intraslab" = intraslab-type in subduction regimes
 952        "interface" = interface-type in subduction regimes
 953        "stable" = shallow-crustal in stable-continental regimes
 954    - region = geographic region where the earthquake occured"
 955        "japan"
 956        "chile"
 957        "other"
 958    - elat = hypocenter latitude (degrees)
 959    - elon = hypocenter longiude (degrees)
 960    - hypd = hypocenter depth (km); positive into the ground
 961    - magnitude = earthquake moment magnitude (Mw)
 962    - method = code for rupture simulation constraints
 963        "A" = strike, dip, and rake for first nodal plane solution preferred
 964              (optional "strike", "dip", and "rake" arguments are required)
 965              Warning: not recommended
 966        "B" = strike, dip, and rake for second nodal plane solution preferred
 967              (optional "strike2", "dip2", and "rake2" arguments are required)
 968              Warning: not recommended
 969        "C" = strike, dip, and rake are known for two nodal planes, and neither is preferred
 970              (optional "strike", "dip", "rake", "strike2", "dip2", and "rake2" arguments are required)
 971        "D" = One nodal plane solution for strike, dip, and rake; randomize the strike and dip
 972              (optional "strike", "dip", and "rake" arguments are required)
 973              Warning: not recommended
 974        "E" = No nodal plane solutions; randomize strike, dip, and rake
 975              (dip and rake are assigned based on faulting mechanism)
 976              (if optional "mechanism" argument is not speficied, simulations randomly assign one)
 977    - nsims = Number of simulations assigned to each M-scaling relationship. Total number of simulations should be odd.
 978        nsims[0] = Wells & Coppersmith (1994) - [recommended 334]
 979        nsims[1] = Leonard (2014) - [recommended 333]
 980        nsims[2] = Thingbaijam et al. (2017) [recommended 333]
 981        nsims[3] = Chiou & Youngs (2008) aspect ratio model with Wells & Coppersmith (1994) area relationship [recommended 111]
 982        nsims[4] = Chiou & Youngs (2008) aspect ratio model with Leonard (2014) area relationship [recommended 111]
 983        nsims[5] = Chiou & Youngs (2008) aspect ratio model with Thingbaijam et al. (2017) area relationship [recommended 111]
 984        nsims[6] = Contreras et al. (2022) - [recommended 333]
 985    - mechanim = known or preferred style-of-faulting [default None]
 986        "SS" = strike-slip: (-180 < rake < -150) or (-30 < rake < 30) or (150 < rake < 180)
 987        "NM" = normal: -150 < rake < -30
 988        "RV" = reverse: 30 < rake < 150
 989    - strike = strike-angle (degrees) of the first nodal plane solution [default None]
 990    - dip = dip-angle (degrees) of the first nodal plane solution [default None]
 991    - rake = rake-angle (degrees) of the first nodal plane solution [default None]
 992    - strike2 = strike-angle (degrees) of the second nodal plane solution [default None]
 993    - dip2 = dip-angle (degrees) of the second nodal plane solution [default None]
 994    - rake2 = rake-angle (degrees) of the second nodal plane solution [default None]
 995
 996    Output Arguments:
 997    - median_simulation = index of selected (median) rupture
 998    - along_strikes = numpy array of the hypocenter along-strike relative position [0, 1]
 999    - down_dips = numpy array of the hypocenter down-dip relative position [0, 1]
1000    - lon1s = numpy array of longitude (degrees) for the ULC rupture vertex
1001    - lat1s = numpy array of latitude (degrees) for the ULC rupture vertex
1002    - rpz1s = numpy array of depth (km) for the ULC rupture vertex
1003    - lon2s = numpy array of longitude (degrees) for the URC rupture vertex
1004    - lat2s = numpy array of latitude (degrees) for the URC rupture vertex
1005    - rpz2s = numpy array of depth (km) for the URC rupture vertex
1006    - lon3s = numpy array of longitude (degrees) for the LLC rupture vertex
1007    - lat3s = numpy array of latitude (degrees) for the LLC rupture vertex
1008    - rpz3s = numpy array of depth (km) for the LLC rupture vertex
1009    - lon4s = numpy array of longitude (degrees) for the LRC rupture vertex
1010    - lat4s = numpy array of latitude (degrees) for the LRC rupture vertex
1011    - rpz4s = numpy array of depth (km) for the LRC rupture vertex
1012    - strikes = numpy array of strike angle (degrees)
1013    - dips = numpy array of dip angle (degrees)
1014    - rakes = numpy array of rake angle (degrees)
1015    - models = numpy array of magnitude-scaling relationship used for each realization
1016    - areas = numpy array of area (km^2)
1017    - aspect_ratios = numpy array of aspect ratio
1018    - lengths = numpy array of rupture length (km)
1019    - widths = numpy array of rupture width (km)
1020    - top_depths = numpy array of top depth (km) of rupture
1021    - bottom_depths = numpy array of bottom depth (km) of rupture
1022    """
1023
1024    check_input_arguments(
1025        eqType, method, nsims, strike, dip, rake, strike2, dip2, rake2
1026    )
1027    try:
1028        strike = float(strike)
1029        dip = float(dip)
1030        rake = float(rake)
1031    except:
1032        pass
1033    try:
1034        strike2 = float(strike2)
1035        dip2 = float(dip2)
1036        rake2 = float(rake2)
1037    except:
1038        pass
1039
1040    # Set random seed -------------------------------------------------------------------------------------------------------
1041    np.random.seed(seed=eqid)
1042    total_sims = np.sum(nsims)
1043
1044    # Set earthquake parameters  --------------------------------------------------------------------------------------------
1045    if method == "A":
1046        # Prefer the first nodal plane solution:
1047        strikes = np.repeat([strike], total_sims)
1048        dips = np.repeat([dip], total_sims)
1049        rakes = np.repeat([rake], total_sims)
1050        mechanisms = np.repeat(get_mechanism_based_on_rake(rake), total_sims)
1051    elif method == "B":
1052        # Prefer the second nodal plane solution:
1053        strikes = np.repeat([strike2], total_sims)
1054        dips = np.repeat([dip2], total_sims)
1055        rakes = np.repeat([rake2], total_sims)
1056        mechanisms = np.repeat(get_mechanism_based_on_rake(rake2), total_sims)
1057    elif method == "C":
1058        # Randomize between the two nodal plane solutions:
1059        nodal_planes = np.random.choice([1, 2], total_sims)
1060        strikes = np.where(nodal_planes == 1, strike, strike2)
1061        dips = np.where(nodal_planes == 1, dip, dip2)
1062        rakes = np.where(nodal_planes == 1, rake, rake2)
1063        mechanisms = v_get_mechanism_based_on_rake(rakes)
1064    elif method == "D":
1065        # Randomize strike and dip with a uniform distributions centered around the perscribed values
1066        strikes = np.random.uniform(
1067            strike - 30.0, strike + 30.0, total_sims
1068        )  # +/- 30 degrees
1069        strikes = np.where(strikes < 0, strikes + 360, strikes)  # strike >= 0
1070        strikes = np.where(strikes >= 360.0, strikes - 360.0, strikes)  # strike < 360
1071        dips = np.random.uniform(dip - 10.0, dip + 10.0, total_sims)  # +/- 10 degrees
1072        dips = np.where(dips < 10.0, 10.0, dips)  # 10 <= dip
1073        dips = np.where(
1074            dips > 89.9999999, 89.9999999, dips
1075        )  # dip <= 90 # cannot divide by zero
1076        rakes = np.repeat([rake], total_sims)
1077        mechanisms = v_get_mechanism_based_on_rake(rakes)
1078    elif method == "E":
1079        # Randomize everything based on prescribed mechanism
1080        # if mechanism is not prescribed, then randomize everything
1081        if mechanism == None:
1082            mechanisms = np.random.choice(["SS", "NM", "RV"], total_sims)
1083        else:
1084            mechanisms = np.repeat(mechanism, total_sims)
1085        strikes = np.random.uniform(0.0, 360.0, total_sims)
1086        strikes = np.where(
1087            strikes >= 360.0, strikes - 360.0, strikes
1088        )  # 0 <= strike < 360
1089        rakes = np.where(mechanisms == "SS", 0.0, mechanisms)
1090        rakes = np.where(rakes == "NM", -90.0, rakes)
1091        rakes = np.where(rakes == "RV", 90.0, rakes)
1092        dips = np.where(
1093            mechanisms == "SS", 89.9999999, mechanisms
1094        )  # cannot divide by zero
1095        dips = np.where(dips == "NM", 55.0, dips)
1096        dips = np.where(dips == "RV", 40.0, dips)
1097        dips = dips.astype(float)
1098    dips = np.where(
1099        dips > 89.9999999, 89.9999999, dips
1100    )  # dip <= 90 # cannot divide by zero
1101
1102    # Convert degrees to radians  -------------------------------------------------------------------------------------------
1103    strikes *= np.pi / 180.0
1104    dips *= np.pi / 180.0
1105
1106    # Simulate hypocenter down-dip and along-strike position  ---------------------------------------------------------------
1107    down_dips = v_get_hyp_down_dip_position(np.repeat(eqType, total_sims), region)
1108    along_strikes = v_get_hyp_along_strike_position(np.repeat(eqType, total_sims))
1109
1110    # Simulate rupture geometry using magnitude-scaling relationships  ------------------------------------------------------
1111    if nsims[0] > 0:
1112        nsim_WellsCoppersmith1994 = nsims[0]
1113        relationships_WellsCoppersmith1994 = v_WellsCoppersmith1994(
1114            np.repeat(magnitude, nsim_WellsCoppersmith1994),
1115            np.repeat(eqType, nsim_WellsCoppersmith1994),
1116        )
1117        areas_WellsCoppersmith1994 = relationships_WellsCoppersmith1994[0]
1118        aspect_ratios_WellsCoppersmith1994 = relationships_WellsCoppersmith1994[1]
1119        lengths_WellsCoppersmith1994 = relationships_WellsCoppersmith1994[2]
1120        widths_WellsCoppersmith1994 = relationships_WellsCoppersmith1994[3]
1121        models = np.repeat("WellsCoppersmith1994", nsim_WellsCoppersmith1994)
1122    else:
1123        areas_WellsCoppersmith1994, aspect_ratios_WellsCoppersmith1994 = np.array(
1124            []
1125        ), np.array([])
1126        lengths_WellsCoppersmith1994, widths_WellsCoppersmith1994 = np.array(
1127            []
1128        ), np.array([])
1129        models = np.array([])
1130
1131    if nsims[1] > 0:
1132        nsim_Leonard2014 = nsims[1]
1133        mechanisms_Leonard2014 = mechanisms[np.sum(nsims[:1]) : np.sum(nsims[:2])]
1134        relationships_Leonard2014 = v_Leonard2014(
1135            np.repeat(magnitude, nsim_Leonard2014),
1136            mechanisms_Leonard2014,
1137            np.repeat(eqType, nsim_Leonard2014),
1138        )
1139        areas_Leonard2014 = relationships_Leonard2014[0]
1140        aspect_ratios_Leonard2014 = relationships_Leonard2014[1]
1141        lengths_Leonard2014 = relationships_Leonard2014[2]
1142        widths_Leonard2014 = relationships_Leonard2014[3]
1143        models = np.concatenate((models, np.repeat("Leonard2014", nsim_Leonard2014)))
1144    else:
1145        areas_Leonard2014, aspect_ratios_Leonard2014 = np.array([]), np.array([])
1146        lengths_Leonard2014, widths_Leonard2014 = np.array([]), np.array([])
1147
1148    if nsims[2] > 0:
1149        nsim_ThingbaijamEtAl2017 = nsims[2]
1150        mechanisms_ThingbaijamEtAl2017 = mechanisms[
1151            np.sum(nsims[:2]) : np.sum(nsims[:3])
1152        ]
1153        relationships_ThingbaijamEtAl2017 = v_ThingbaijamEtAl2017(
1154            np.repeat(magnitude, nsim_ThingbaijamEtAl2017),
1155            mechanisms_ThingbaijamEtAl2017,
1156            np.repeat(eqType, nsim_ThingbaijamEtAl2017),
1157        )
1158        areas_ThingbaijamEtAl2017 = relationships_ThingbaijamEtAl2017[0]
1159        aspect_ratios_ThingbaijamEtAl2017 = relationships_ThingbaijamEtAl2017[1]
1160        lengths_ThingbaijamEtAl2017 = relationships_ThingbaijamEtAl2017[2]
1161        widths_ThingbaijamEtAl2017 = relationships_ThingbaijamEtAl2017[3]
1162        models = np.concatenate(
1163            (models, np.repeat("ThingbaijamEtAl2017", nsim_ThingbaijamEtAl2017))
1164        )
1165    else:
1166        areas_ThingbaijamEtAl2017, aspect_ratios_ThingbaijamEtAl2017 = np.array(
1167            []
1168        ), np.array([])
1169        lengths_ThingbaijamEtAl2017, widths_ThingbaijamEtAl2017 = np.array(
1170            []
1171        ), np.array([])
1172
1173    if nsims[3] > 0:
1174        nsim_ChiouYoungs2008_WellsCoppersmith1994 = nsims[3]
1175        mechanisms_ChiouYoungs2008_WellsCoppersmith1994 = mechanisms[
1176            np.sum(nsims[:3]) : np.sum(nsims[:4])
1177        ]
1178        relationships_ChiouYoungs2008_WellsCoppersmith1994 = v_ChiouYoungs2008(
1179            np.repeat(magnitude, nsim_ChiouYoungs2008_WellsCoppersmith1994),
1180            mechanisms_ChiouYoungs2008_WellsCoppersmith1994,
1181            np.repeat(eqType, nsim_ChiouYoungs2008_WellsCoppersmith1994),
1182            "WellsCoppersmith1994",
1183        )
1184        areas_ChiouYoungs2008_WellsCoppersmith1994 = (
1185            relationships_ChiouYoungs2008_WellsCoppersmith1994[0]
1186        )
1187        aspect_ratios_ChiouYoungs2008_WellsCoppersmith1994 = (
1188            relationships_ChiouYoungs2008_WellsCoppersmith1994[1]
1189        )
1190        lengths_ChiouYoungs2008_WellsCoppersmith1994 = (
1191            relationships_ChiouYoungs2008_WellsCoppersmith1994[2]
1192        )
1193        widths_ChiouYoungs2008_WellsCoppersmith1994 = (
1194            relationships_ChiouYoungs2008_WellsCoppersmith1994[3]
1195        )
1196        models = np.concatenate(
1197            (
1198                models,
1199                np.repeat(
1200                    "ChiouYoungs2008_WellsCoppersmith1994",
1201                    nsim_ChiouYoungs2008_WellsCoppersmith1994,
1202                ),
1203            )
1204        )
1205    else:
1206        (
1207            areas_ChiouYoungs2008_WellsCoppersmith1994,
1208            aspect_ratios_ChiouYoungs2008_WellsCoppersmith1994,
1209        ) = np.array([]), np.array([])
1210        (
1211            lengths_ChiouYoungs2008_WellsCoppersmith1994,
1212            widths_ChiouYoungs2008_WellsCoppersmith1994,
1213        ) = np.array([]), np.array([])
1214
1215    if nsims[4] > 0:
1216        nsim_ChiouYoungs2008_Leonard2014 = nsims[4]
1217        mechanisms_ChiouYoungs2008_Leonard2014 = mechanisms[
1218            np.sum(nsims[:4]) : np.sum(nsims[:5])
1219        ]
1220        relationships_ChiouYoungs2008_Leonard2014 = v_ChiouYoungs2008(
1221            np.repeat(magnitude, nsim_ChiouYoungs2008_Leonard2014),
1222            mechanisms_ChiouYoungs2008_Leonard2014,
1223            np.repeat(eqType, nsim_ChiouYoungs2008_Leonard2014),
1224            "Leonard2014",
1225        )
1226        areas_ChiouYoungs2008_Leonard2014 = relationships_ChiouYoungs2008_Leonard2014[0]
1227        aspect_ratios_ChiouYoungs2008_Leonard2014 = (
1228            relationships_ChiouYoungs2008_Leonard2014[1]
1229        )
1230        lengths_ChiouYoungs2008_Leonard2014 = relationships_ChiouYoungs2008_Leonard2014[
1231            2
1232        ]
1233        widths_ChiouYoungs2008_Leonard2014 = relationships_ChiouYoungs2008_Leonard2014[
1234            3
1235        ]
1236        models = np.concatenate(
1237            (
1238                models,
1239                np.repeat(
1240                    "ChiouYoungs2008_Leonard2014", nsim_ChiouYoungs2008_Leonard2014
1241                ),
1242            )
1243        )
1244    else:
1245        areas_ChiouYoungs2008_Leonard2014, aspect_ratios_ChiouYoungs2008_Leonard2014 = (
1246            np.array([]),
1247            np.array([]),
1248        )
1249        lengths_ChiouYoungs2008_Leonard2014, widths_ChiouYoungs2008_Leonard2014 = (
1250            np.array([]),
1251            np.array([]),
1252        )
1253
1254    if nsims[5] > 0:
1255        nsim_ChiouYoungs2008_ThingbaijamEtAl2017 = nsims[5]
1256        mechanisms_ChiouYoungs2008_ThingbaijamEtAl2017 = mechanisms[
1257            np.sum(nsims[:5]) : np.sum(nsims[:6])
1258        ]
1259        relationships_ChiouYoungs2008_ThingbaijamEtAl2017 = v_ChiouYoungs2008(
1260            np.repeat(magnitude, nsim_ChiouYoungs2008_ThingbaijamEtAl2017),
1261            mechanisms_ChiouYoungs2008_ThingbaijamEtAl2017,
1262            np.repeat(eqType, nsim_ChiouYoungs2008_ThingbaijamEtAl2017),
1263            "ThingbaijamEtAl2017",
1264        )
1265        areas_ChiouYoungs2008_ThingbaijamEtAl2017 = (
1266            relationships_ChiouYoungs2008_ThingbaijamEtAl2017[0]
1267        )
1268        aspect_ratios_ChiouYoungs2008_ThingbaijamEtAl2017 = (
1269            relationships_ChiouYoungs2008_ThingbaijamEtAl2017[1]
1270        )
1271        lengths_ChiouYoungs2008_ThingbaijamEtAl2017 = (
1272            relationships_ChiouYoungs2008_ThingbaijamEtAl2017[2]
1273        )
1274        widths_ChiouYoungs2008_ThingbaijamEtAl2017 = (
1275            relationships_ChiouYoungs2008_ThingbaijamEtAl2017[3]
1276        )
1277        models = np.concatenate(
1278            (
1279                models,
1280                np.repeat(
1281                    "ChiouYoungs2008_ThingbaijamEtAl2017",
1282                    nsim_ChiouYoungs2008_ThingbaijamEtAl2017,
1283                ),
1284            )
1285        )
1286    else:
1287        (
1288            areas_ChiouYoungs2008_ThingbaijamEtAl2017,
1289            aspect_ratios_ChiouYoungs2008_ThingbaijamEtAl2017,
1290        ) = np.array([]), np.array([])
1291        (
1292            lengths_ChiouYoungs2008_ThingbaijamEtAl2017,
1293            widths_ChiouYoungs2008_ThingbaijamEtAl2017,
1294        ) = np.array([]), np.array([])
1295
1296    if nsims[6] > 0:
1297        nsim_ContrerasEtAl2022 = nsims[6]
1298        mechanisms_ContrerasEtAl2022 = mechanisms[np.sum(nsims[:6]) : np.sum(nsims[:7])]
1299        relationships_ContrerasEtAl2022 = v_ContrerasEtAl2022(
1300            np.repeat(magnitude, nsim_ContrerasEtAl2022),
1301            np.repeat(eqType, nsim_ContrerasEtAl2022),
1302        )
1303        areas_ContrerasEtAl2022 = relationships_ContrerasEtAl2022[0]
1304        aspect_ratios_ContrerasEtAl2022 = relationships_ContrerasEtAl2022[1]
1305        lengths_ContrerasEtAl2022 = relationships_ContrerasEtAl2022[2]
1306        widths_ContrerasEtAl2022 = relationships_ContrerasEtAl2022[3]
1307        models = np.concatenate(
1308            (models, np.repeat("ContrerasEtAl2022", nsim_ContrerasEtAl2022))
1309        )
1310    else:
1311        areas_ContrerasEtAl2022, aspect_ratios_ContrerasEtAl2022 = np.array(
1312            []
1313        ), np.array([])
1314        lengths_ContrerasEtAl2022, widths_ContrerasEtAl2022 = np.array([]), np.array([])
1315
1316    areas = np.concatenate(
1317        (
1318            areas_WellsCoppersmith1994,
1319            areas_Leonard2014,
1320            areas_ThingbaijamEtAl2017,
1321            areas_ChiouYoungs2008_WellsCoppersmith1994,
1322            areas_ChiouYoungs2008_Leonard2014,
1323            areas_ChiouYoungs2008_ThingbaijamEtAl2017,
1324            areas_ContrerasEtAl2022,
1325        )
1326    )
1327    aspect_ratios = np.concatenate(
1328        (
1329            aspect_ratios_WellsCoppersmith1994,
1330            aspect_ratios_Leonard2014,
1331            aspect_ratios_ThingbaijamEtAl2017,
1332            aspect_ratios_ChiouYoungs2008_WellsCoppersmith1994,
1333            aspect_ratios_ChiouYoungs2008_Leonard2014,
1334            aspect_ratios_ChiouYoungs2008_ThingbaijamEtAl2017,
1335            aspect_ratios_ContrerasEtAl2022,
1336        )
1337    )
1338    lengths = np.concatenate(
1339        (
1340            lengths_WellsCoppersmith1994,
1341            lengths_Leonard2014,
1342            lengths_ThingbaijamEtAl2017,
1343            lengths_ChiouYoungs2008_WellsCoppersmith1994,
1344            lengths_ChiouYoungs2008_Leonard2014,
1345            lengths_ChiouYoungs2008_ThingbaijamEtAl2017,
1346            lengths_ContrerasEtAl2022,
1347        )
1348    )
1349    widths = np.concatenate(
1350        (
1351            widths_WellsCoppersmith1994,
1352            widths_Leonard2014,
1353            widths_ThingbaijamEtAl2017,
1354            widths_ChiouYoungs2008_WellsCoppersmith1994,
1355            widths_ChiouYoungs2008_Leonard2014,
1356            widths_ChiouYoungs2008_ThingbaijamEtAl2017,
1357            widths_ContrerasEtAl2022,
1358        )
1359    )
1360
1361    # Compute geometric parameters  -----------------------------------------------------------------------------------------
1362    xf = np.sin(strikes) * lengths * along_strikes
1363    yf = np.cos(strikes) * lengths * along_strikes
1364    xb = np.sin(strikes + np.pi) * lengths * (1.0 - along_strikes)
1365    yb = np.cos(strikes + np.pi) * lengths * (1.0 - along_strikes)
1366    rwh = widths * np.cos(dips)
1367    rwv = widths * np.sin(dips)
1368    top_depths = hypd - rwv * down_dips
1369    down_dips = np.where(
1370        top_depths < 0.0, hypd / rwv, down_dips
1371    )  # make sure hypocenter is in the ground
1372    top_depths = np.where(top_depths < 0.0, 0.0, top_depths)  # top_depth > 0.0
1373    bottom_depths = top_depths + rwv
1374
1375    # Compute top points of the rupture surface 1 = ULC & 2 = URC  ----------------------------------------------------------
1376    rpx1s = xf + np.sin(strikes - np.arcsin(1.0)) * rwh * down_dips
1377    rpy1s = yf + np.cos(strikes - np.arcsin(1.0)) * rwh * down_dips
1378    rpz1s = top_depths
1379
1380    rpx2s = xb + np.sin(strikes - np.arcsin(1.0)) * rwh * down_dips
1381    rpy2s = yb + np.cos(strikes - np.arcsin(1.0)) * rwh * down_dips
1382    rpz2s = top_depths
1383
1384    # Compute bottom points of the rupture surface 3 = LRC & 4 = LLC  -------------------------------------------------------
1385    rpx3s = xf + np.sin(strikes + np.arcsin(1.0)) * rwh * (1.0 - down_dips)
1386    rpy3s = yf + np.cos(strikes + np.arcsin(1.0)) * rwh * (1.0 - down_dips)
1387    rpz3s = bottom_depths
1388
1389    rpx4s = xb + np.sin(strikes + np.arcsin(1.0)) * rwh * (1.0 - down_dips)
1390    rpy4s = yb + np.cos(strikes + np.arcsin(1.0)) * rwh * (1.0 - down_dips)
1391    rpz4s = bottom_depths
1392
1393    # Create a psudo-grid of stations  --------------------------------------------------------------------------------------
1394    r = np.array(
1395        list(np.arange(2, 20, 2))
1396        + list(np.arange(25, 55, 5))
1397        + list(np.arange(60, 110, 10))
1398        + list(np.arange(125, 325, 25))
1399    )
1400    theta = np.pi * np.linspace(0, 360, 25)[:-1] / 180
1401    psx = np.multiply(
1402        np.repeat(r.reshape(1, len(r)), len(theta), axis=0),
1403        np.sin(theta.reshape(len(theta), 1)),
1404    ).flatten()
1405    psy = np.multiply(
1406        np.repeat(r.reshape(1, len(r)), len(theta), axis=0),
1407        np.cos(theta.reshape(len(theta), 1)),
1408    ).flatten()
1409
1410    # Compute distances to each psudo-station  ------------------------------------------------------------------------------
1411    SITES = np.array([[[x, y, 0]] for x, y in zip(psx, psy)])
1412    TRI_123 = np.array(
1413        [
1414            [[x1, y1, z1], [x2, y2, z2], [x3, y3, z3]]
1415            for x1, x2, x3, y1, y2, y3, z1, z2, z3 in zip(
1416                rpx1s, rpx2s, rpx3s, rpy1s, rpy2s, rpy3s, rpz1s, rpz2s, rpz3s
1417            )
1418        ]
1419    )
1420    TRI_234 = np.array(
1421        [
1422            [[x1, y1, z1], [x2, y2, z2], [x3, y3, z3]]
1423            for x1, x2, x3, y1, y2, y3, z1, z2, z3 in zip(
1424                rpx2s, rpx4s, rpx3s, rpy2s, rpy4s, rpy3s, rpz2s, rpz4s, rpz3s
1425            )
1426        ]
1427    )
1428    rrups1 = pointTriangleDistance(TRI_123, SITES)
1429    rrups2 = pointTriangleDistance(TRI_234, SITES)
1430    rrups = np.minimum(rrups1, rrups2)
1431
1432    # Locate the median Rrup for each site  ---------------------------------------------------------------------------------
1433    median_rrups = np.median(rrups, axis=0)
1434    median_rrup_indices = v_get_median_index(np.transpose(rrups))
1435
1436    # Locate the simulation that minimizes the difference between median Rrup and Rrup across all sites  --------------------
1437    sum_of_difference_squared = np.sum((median_rrups - rrups) ** 2, axis=1)
1438    median_simulation = np.argmin(sum_of_difference_squared)
1439
1440    # Convert simulated faults xy to ll using haversine formulation  --------------------------------------------------------
1441    lon1s, lat1s = v_XY2LL(elon, elat, rpx1s, rpy1s)
1442    lon2s, lat2s = v_XY2LL(elon, elat, rpx2s, rpy2s)
1443    lon3s, lat3s = v_XY2LL(elon, elat, rpx3s, rpy3s)
1444    lon4s, lat4s = v_XY2LL(elon, elat, rpx4s, rpy4s)
1445
1446    # Convert radians to degrees  -------------------------------------------------------------------------------------------
1447    strikes = np.round(strikes * 180 / np.pi, 6)
1448    dips = np.round(dips * 180 / np.pi, 6)
1449
1450    # Return everything  ----------------------------------------------------------------------------------------------------
1451    return (
1452        median_simulation,
1453        along_strikes,
1454        down_dips,
1455        lon1s,
1456        lat1s,
1457        rpz1s,
1458        lon2s,
1459        lat2s,
1460        rpz2s,
1461        lon3s,
1462        lat3s,
1463        rpz3s,
1464        lon4s,
1465        lat4s,
1466        rpz4s,
1467        strikes,
1468        dips,
1469        rakes,
1470        models,
1471        areas,
1472        aspect_ratios,
1473        lengths,
1474        widths,
1475        top_depths,
1476        bottom_depths,
1477    )
1478
1479
1480def simulate_rupture_surface(
1481    eqid,
1482    eqType,
1483    region,
1484    elat,
1485    elon,
1486    hypd,
1487    magnitude,
1488    method,
1489    nsims,
1490    mechanism=None,
1491    strike=None,
1492    dip=None,
1493    rake=None,
1494    strike2=None,
1495    dip2=None,
1496    rake2=None,
1497):
1498    """
1499    Simulate earthquake rupture surface that minimizes the difference between the median distance of a pseudo-grid of sites
1500    and a stochastic set of possible ruptures. Randomized fault attributes are dependent on the fault category. Fault scaling
1501    models depend on earthquake type. Adapted from the CCLD routine originally coded in Frortran by Chiou and Youngs (2008).
1502
1503    Input Arguments:
1504    - eqid = unique integer identifier for the event (to set the random seed)
1505    - eqType = type of earthquake:
1506        "crustal" = shallow-crustal in active tectonic regimes
1507        "intraslab" = intraslab-type in subduction regimes
1508        "interface" = interface-type in subduction regimes
1509        "stable" = shallow-crustal in stable-continental regimes
1510    - region = geographic region where the earthquake occured"
1511        "japan"
1512        "chile"
1513        "other"
1514    - elat = hypocenter latitude (degrees)
1515    - elon = hypocenter longiude (degrees)
1516    - hypd = hypocenter depth (km); positive into the ground
1517    - magnitude = earthquake moment magnitude (Mw)
1518    - method = code for rupture simulation constraints
1519        "A" = strike, dip, and rake for first nodal plane solution preferred
1520              (optional "strike", "dip", and "rake" arguments are required)
1521              Warning: not recommended
1522        "B" = strike, dip, and rake for second nodal plane solution preferred
1523              (optional "strike2", "dip2", and "rake2" arguments are required)
1524              Warning: not recommended
1525        "C" = strike, dip, and rake are known for two nodal planes, and neither is preferred
1526              (optional "strike", "dip", "rake", "strike2", "dip2", and "rake2" arguments are required)
1527        "D" = One nodal plane solution for strike, dip, and rake; randomize the strike and dip
1528              (optional "strike", "dip", and "rake" arguments are required)
1529              Warning: not recommended
1530        "E" = No nodal plane solutions; randomize strike, dip, and rake
1531              (dip and rake are assigned based on faulting mechanism)
1532              (if optional "mechanism" argument is not speficied, simulations randomly assign one)
1533    - nsims = Number of simulations assigned to each M-scaling relationship. Total number of simulations should be odd.
1534        nsims[0] = Wells & Coppersmith (1994) - [recommended 334]
1535        nsims[1] = Leonard (2014) - [recommended 333]
1536        nsims[2] = Thingbaijam et al. (2017) [recommended 333]
1537        nsims[3] = Chiou & Youngs (2008) aspect ratio model with Wells & Coppersmith (1994) area relationship [recommended 111]
1538        nsims[4] = Chiou & Youngs (2008) aspect ratio model with Leonard (2014) area relationship [recommended 111]
1539        nsims[5] = Chiou & Youngs (2008) aspect ratio model with Thingbaijam et al. (2017) area relationship [recommended 111]
1540        nsims[6] = Contreras et al. (2022) - [recommended 333]
1541    - mechanim = known or preferred style-of-faulting [default None]
1542        "SS" = strike-slip: (-180 < rake < -150) or (-30 < rake < 30) or (150 < rake < 180)
1543        "NM" = normal: -150 < rake < -30
1544        "RV" = reverse: 30 < rake < 150
1545    - strike = strike-angle (degrees) of the first nodal plane solution [default None]
1546    - dip = dip-angle (degrees) of the first nodal plane solution [default None]
1547    - rake = rake-angle (degrees) of the first nodal plane solution [default None]
1548    - strike2 = strike-angle (degrees) of the second nodal plane solution [default None]
1549    - dip2 = dip-angle (degrees) of the second nodal plane solution [default None]
1550    - rake2 = rake-angle (degrees) of the second nodal plane solution [default None]
1551
1552    Output Objects:
1553    - SIMULATIONS = pandas DataFrame object containing all simulated rupture surfaces
1554    - SELECTED = pandas DataFrame object containing the selected rupture surface and statistics
1555    """
1556
1557    (
1558        median_simulation,
1559        along_strikes,
1560        down_dips,
1561        lon1s,
1562        lat1s,
1563        rpz1s,
1564        lon2s,
1565        lat2s,
1566        rpz2s,
1567        lon3s,
1568        lat3s,
1569        rpz3s,
1570        lon4s,
1571        lat4s,
1572        rpz4s,
1573        strikes,
1574        dips,
1575        rakes,
1576        models,
1577        areas,
1578        aspect_ratios,
1579        lengths,
1580        widths,
1581        top_depths,
1582        bottom_depths,
1583    ) = get_rupture_surface_simulations(
1584        eqid,
1585        eqType,
1586        region,
1587        elat,
1588        elon,
1589        hypd,
1590        magnitude,
1591        method,
1592        nsims,
1593        mechanism,
1594        strike,
1595        dip,
1596        rake,
1597        strike2,
1598        dip2,
1599        rake2,
1600    )
1601
1602    # Compute the statistics of the simulation results ----------------------------------------------------------------------
1603    areas_avg = 10 ** np.mean(np.log10(areas))
1604    areas_std = 10 ** np.std(np.log10(areas))
1605    aspect_ratios_avg = 10 ** np.mean(np.log10(aspect_ratios))
1606    aspect_ratios_std = 10 ** np.std(np.log10(aspect_ratios))
1607    lengths_avg = 10 ** np.mean(np.log10(lengths))
1608    lengths_std = 10 ** np.std(np.log10(lengths))
1609    widths_avg = 10 ** np.mean(np.log10(widths))
1610    widths_std = 10 ** np.std(np.log10(widths))
1611
1612    max_top_depths = np.max(top_depths)
1613    min_top_depths = np.min(top_depths)
1614    max_bottom_depths = np.max(bottom_depths)
1615    min_bottom_depths = np.min(bottom_depths)
1616
1617    # Construct pandas DataFrame objects to package the results -------------------------------------------------------------
1618    total_sims = np.sum(nsims)
1619    SIMULATIONS = pd.DataFrame(
1620        {
1621            "Simulation": np.array(range(total_sims)) + 1,
1622            "EQID": np.repeat(eqid, total_sims),
1623            "Magnitude": np.repeat(magnitude, total_sims),
1624            "Hypocenter Longitude": np.repeat(elon, total_sims),
1625            "Hypocenter Latitude": np.repeat(elat, total_sims),
1626            "Hypocenter Depth (km)": np.repeat(hypd, total_sims),
1627            "Hypocenter Along-Strike Position": along_strikes,
1628            "Hypocenter Down-Dip Position": down_dips,
1629            "ULC Longitude": lon2s,
1630            "ULC Latitude": lat2s,
1631            "ULC Depth (km)": rpz2s,
1632            "URC Longitude": lon1s,
1633            "URC Latitude": lat1s,
1634            "URC Depth (km)": rpz1s,
1635            "LRC Longitude": lon3s,
1636            "LRC Latitude": lat3s,
1637            "LRC Depth (km)": rpz3s,
1638            "LLC Longitude": lon4s,
1639            "LLC Latitude": lat4s,
1640            "LLC Depth (km)": rpz4s,
1641            "Strike": strikes,
1642            "Dip": dips,
1643            "Rake": rakes,
1644            "Scaling Relation": models,
1645            "Area (km^2)": areas,
1646            "Aspect Ratio": aspect_ratios,
1647            "Rupture Length (km)": lengths,
1648            "Rupture Width (km)": widths,
1649            "Rupture Top Depth (km)": top_depths,
1650            "Rupture Bottom Depth (km)": bottom_depths,
1651        }
1652    )
1653
1654    SELECTED = pd.DataFrame(
1655        {
1656            "Median Simulation": [median_simulation + 1],
1657            "EQID": [eqid],
1658            "Magnitude": [magnitude],
1659            "Hypocenter Longitude": [elon],
1660            "Hypocenter Latitude": [elat],
1661            "Hypocenter Depth (km)": [hypd],
1662            "Hypocenter Along-Strike Position": [along_strikes[median_simulation]],
1663            "Hypocenter Down-Dip Position": [down_dips[median_simulation]],
1664            "ULC Longitude": [lon2s[median_simulation]],
1665            "ULC Latitude": [lat2s[median_simulation]],
1666            "ULC Depth (km)": [rpz2s[median_simulation]],
1667            "URC Longitude": [lon1s[median_simulation]],
1668            "URC Latitude": [lat1s[median_simulation]],
1669            "URC Depth (km)": [rpz1s[median_simulation]],
1670            "LRC Longitude": [lon3s[median_simulation]],
1671            "LRC Latitude": [lat3s[median_simulation]],
1672            "LRC Depth (km)": [rpz3s[median_simulation]],
1673            "LLC Longitude": [lon4s[median_simulation]],
1674            "LLC Latitude": [lat4s[median_simulation]],
1675            "LLC Depth (km)": [rpz4s[median_simulation]],
1676            "Strike": [strikes[median_simulation]],
1677            "Dip": [dips[median_simulation]],
1678            "Rake": [rakes[median_simulation]],
1679            "Scaling Relation": [models[median_simulation]],
1680            "Area (km^2)": [areas[median_simulation]],
1681            "Average Area (km^2)": [areas_avg],
1682            "Area Standard Deviation (km^2)": [areas_std],
1683            "Aspect Ratio": [aspect_ratios[median_simulation]],
1684            "Average Aspect Ratio": [aspect_ratios_avg],
1685            "Aspect Ratio Standard Deviation": [aspect_ratios_std],
1686            "Length (km)": [lengths[median_simulation]],
1687            "Average Length (km)": [lengths_avg],
1688            "Length Standard Deviation (km)": [lengths_std],
1689            "Width (km)": [widths[median_simulation]],
1690            "Average Width (km)": [widths_avg],
1691            "Width Standard Deviation (km)": [widths_std],
1692            "Rupture Top Depth (km)": [top_depths[median_simulation]],
1693            "Rupture Bottom Depth (km)": [bottom_depths[median_simulation]],
1694        }
1695    )
1696
1697    return (SIMULATIONS, SELECTED)