<reponame>Jetli100/SHRA_Chapter6_7_9
# Note: The codes were originally created by Prof. <NAME> in the MATLAB

import numpy as np
from scipy.stats import norm
from scipy.interpolate import interp1d

# CY_14 attenuation equation - Update of the Chiou and Youngs NGA Model for
# the Average Horizontal Component of Peak Ground Motion and Response
# Spectra, Chiou and Youngs (2014)


##############
### Inputs ###
##############

# M     = Moment Magnitude
# T     = Period (sec); Use Period = -1 for PGV computation
#                 Use 1000 for output the array of Sa with original period
#                 (no interpolation)
# Rup   = Closest distance (km) to the ruptured plane
# Rjb   = Joyner-Boore distance (km); closest distance (km) to surface
#       projection of rupture plane
# Rx    =    Horizontal distance from top of rupture measured perpendicular 
#       to fault strike (km). See Figures a, b and c for illustation% Ztor  = Depth(km) to the top of ruptured plane
# delta = Fault dip angle (in degree)
# lamda = Rake angle      (in degree)
# region        = 0 for global (incl. Taiwan)
#               = 1 for California
#               = 2 for Japan
#               = 3 for China 
#               = 4 for Italy 
#               = 5 for Turkey
# Z10            = Basin depth (km); depth from the groundsurface to the
#                   1km/s shear-wave horizon.
#               = 999 if unknown
# Vs30          = shear wave velocity averaged over top 30 m in m/s
#               = ref: 1130
# FVS30         = 1 for Vs30 is inferred from geology
#               = 0 for measured  Vs30
# d_DPP         = DPP centered on the site and earthquake specific average
#               DPP, = 0 for median calc
#               %d_DPP functionality added BAB 26 July 2020

###############
### Outputs ###
###############

# Sa: Median spectral acceleration prediction
# sigma: logarithmic standard deviation of spectral acceleration
#          prediction


def CY_2014_sub(M, ip, R_RUP, R_JB, Rx, Ztor, delta, F_RV, F_NM, HW, Z10, Vs30, FVS30, region, d_DPP):
    
    dict = {}
    dict["c2"] = [1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06]
    dict["c4"] = [-2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1]
    dict["c4_a"] = [-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5]
    dict["c_RB"] = [50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50]
    dict["c8"] = [0.2154, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.0991, 0.1982, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154]
    dict["c8_a"] = [0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695]
    dict["c1"] = [ 2.3549, -1.5065, -1.5065, -1.4798, -1.2972, -1.1007, -0.9292, -0.658 , -0.5613, -0.5342, -0.5462, -0.5858, -0.6798, -0.8663, -1.0514, -1.3794, -1.6508, -2.1511, -2.5365, -3.0686, -3.4148, -3.9013, -4.2466, -4.5143, -5.0009, -5.3461]
    dict["c1_a"] = [0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165, 0.165, 0.165, 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.1645, 0.1168, 0.0732, 0.0484, 0.022 , 0.0124]
    dict["c1_b"] = [-0.0626, -0.255 , -0.255 , -0.255 , -0.255 , -0.255 , -0.255 , -0.254 , -0.253 , -0.252 , -0.25  , -0.248 , -0.2449, -0.2382, -0.2313, -0.2146, -0.1972, -0.162 , -0.14  , -0.1184, -0.11  , -0.104 , -0.102 , -0.101 , -0.101 , -0.1   ]
    dict["c1_c"] = [-0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.1645, -0.1168, -0.0732, -0.0484, -0.022 , -0.0124]
    dict["c1_d"] = [0.0626, 0.255 , 0.255 , 0.255 , 0.255 , 0.255 , 0.255 , 0.254 , 0.253 , 0.252 , 0.25  , 0.248 , 0.2449, 0.2382, 0.2313, 0.2146, 0.1972, 0.162 , 0.14  , 0.1184, 0.11  , 0.104 , 0.102 , 0.101 , 0.101 , 0.1   ]
    dict["c_n"] = [ 3.3024, 16.0875, 16.0875, 15.7118, 15.8819, 16.4556, 17.6453, 20.1772, 19.9992, 18.7106, 16.6246, 15.3709, 13.7012, 11.2667, 9.1908,  6.5459,  5.2305,  3.7896,  3.3024,  2.8498,  2.5417, 2.1488,  1.8957,  1.7228,  1.5737,  1.5265]
    dict["c_m"] = [5.423 , 4.9993, 4.9993, 4.9993, 4.9993, 4.9993, 4.9993, 5.0031, 5.0172, 5.0315, 5.0547, 5.0704, 5.0939, 5.1315, 5.167 , 5.2317, 5.2893, 5.4109, 5.5106, 5.6705, 5.7981, 5.9983, 6.1552, 6.2856, 6.5428, 6.7415]
    dict["c3"] = [2.3152, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9795, 2.0362, 2.0823, 2.1521, 2.2574, 2.344 , 2.4709, 2.5567, 2.6812, 2.7474, 2.8161, 2.8514, 2.8875, 2.9058, 2.9169, 2.932 , 2.9396]
    dict["c5"] = [5.8096, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.8305, 7.1333, 7.3621, 7.4365, 7.4972, 7.5416, 7.56  , 7.5735, 7.5778, 7.5808, 7.5814, 7.5817, 7.5818, 7.5818, 7.5818, 7.5818, 7.5818, 7.5818]
    dict["c_HM"] = [3.0514, 3.0956, 3.0956, 3.0963, 3.0974, 3.0988, 3.1011, 3.1094, 3.2381, 3.3407, 3.43  , 3.4688, 3.5146, 3.5746, 3.6232, 3.6945, 3.7401, 3.7941, 3.8144, 3.8284, 3.833 , 3.8361, 3.8369, 3.8376, 3.838 , 3.838 ]
    dict["c6"] = [0.4407, 0.4908, 0.4908, 0.4925, 0.4992, 0.5037, 0.5048, 0.5048, 0.5048, 0.5048, 0.5045, 0.5036, 0.5016, 0.4971, 0.4919, 0.4807, 0.4707, 0.4575, 0.4522, 0.4501, 0.45  , 0.45  , 0.45  , 0.45  , 0.45  , 0.45  ]
    dict["c7"] = [0.0324, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.016 , 0.0062, 0.0029, 0.0007, 0.0003]
    dict["c7_b"] = [ 0.0097,  0.0462,  0.0462,  0.0472,  0.0533,  0.0596,  0.0639, 0.063 ,  0.0532,  0.0452,  0.0345,  0.0283, 0.0202,  0.009 , -0.0004, -0.0155, -0.0278, -0.0477, -0.0559, -0.063 , -0.0665, -0.0516, -0.0448, -0.0424, -0.0348, -0.0253]
    dict["c8_b"] = [5.    , 0.4833, 0.4833, 1.2144, 1.6421, 1.9456, 2.181 , 2.6087, 2.9122, 3.1045, 3.3399, 3.4719, 3.6434, 3.8787, 4.0711, 4.3745, 4.6099, 5.0376, 5.3411, 5.7688, 6.0723, 6.5   , 6.8035, 7.0389, 7.4666, 7.77  ]
    dict["c9"] = [0.3079, 0.9228, 0.9228, 0.9296, 0.9396, 0.9661, 0.9794, 1.026 , 1.0177, 1.0008, 0.9801, 0.9652, 0.9459, 0.9196, 0.8829, 0.8302, 0.7884, 0.6754, 0.6196, 0.5101, 0.3917, 0.1244, 0.0086, 0.    , 0.    , 0.    ]
    dict["c9_a"] = [0.1   , 0.1202, 0.1202, 0.1217, 0.1194, 0.1166, 0.1176, 0.1171, 0.1146, 0.1128, 0.1106, 0.115 , 0.1208, 0.1208, 0.1175, 0.106 , 0.1061, 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   ]
    dict["c9_b"] = [6.5   , 6.8607, 6.8607, 6.8697, 6.9113, 7.0271, 7.0959, 7.3298, 7.2588, 7.2372, 7.2109, 7.2491, 7.2988, 7.3691, 6.8789, 6.5334, 6.526 , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   ]
    dict["c11"] = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
    dict["c11_b"] = [-0.3834, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.444 , -0.3539, -0.2688, -0.1793, -0.1428, -0.1138, -0.1062, -0.102 , -0.1009, -0.1003, -0.1001, -0.1001, -0.1, -0.1]
    dict["c_g1"] = [-0.001852, -0.007146, -0.007146, -0.007249, -0.007869, -0.008316, -0.008743, -0.009537, -0.00983 , -0.009913, -0.009896, -0.009787, -0.009505, -0.008918, -0.008251, -0.007267, -0.006492, -0.005147, -0.004277, -0.002979, -0.002301, -0.001344, -0.001084, -0.00101 , -0.000964, -0.00095 ]
    dict["c_g2"] = [-0.007403, -0.006758, -0.006758, -0.006758, -0.006758, -0.006758, -0.006758, -0.00619 , -0.005332, -0.004732, -0.003806, -0.00328 , -0.00269 , -0.002128, -0.001812, -0.001274, -0.001074, -0.001115, -0.001197, -0.001675, -0.002349, -0.003306, -0.003566, -0.00364 , -0.003686, -0.0037  ]
    dict["c_g3"] = [4.3439, 4.2542, 4.2542, 4.2386, 4.2519, 4.296 , 4.3578, 4.5455, 4.7603, 4.8963, 5.0644, 5.1371, 5.188 , 5.2164, 5.1954, 5.0899, 4.7854, 4.3304, 4.1667, 4.0029, 3.8949, 3.7928, 3.7443, 3.709 , 3.6632, 3.623 ]
    dict["phi1"] = [-0.7936, -0.521 , -0.521 , -0.5055, -0.4368, -0.3752, -0.3469, -0.3747, -0.444 , -0.4895, -0.5477, -0.5922, -0.6693, -0.7766, -0.8501, -0.9431, -1.0044, -1.0602, -1.0941, -1.1142, -1.1154, -1.1081, -1.0603, -0.9872, -0.8274, -0.7053]
    dict["phi2"] = [-0.0699, -0.1417, -0.1417, -0.1364, -0.1403, -0.1591, -0.1862, -0.2538, -0.2943, -0.3077, -0.3113, -0.3062, -0.2927, -0.2662, -0.2405, -0.1975, -0.1633, -0.1028, -0.0699, -0.0425, -0.0302, -0.0129, -0.0016,  0.    ,  0.    ,  0.    ]
    dict["phi3"] = [-0.008444, -0.00701 , -0.00701 , -0.007279, -0.007354, -0.006977, -0.006467, -0.005734, -0.005604, -0.005696, -0.005845, -0.005959, -0.006141, -0.006439, -0.006704, -0.007125, -0.007435, -0.00812 , -0.008444, -0.007707, -0.004792, -0.001828, -0.001523, -0.00144 , -0.001369, -0.001361]
    dict["phi4"] = [5.41000e+00, 1.02151e-01, 1.02151e-01, 1.08360e-01, 1.19888e-01, 1.33641e-01, 1.48927e-01, 1.90596e-01, 2.30662e-01, 2.53169e-01, 2.66468e-01, 2.65060e-01, 2.55253e-01, 2.31541e-01, 2.07277e-01, 1.65464e-01, 1.33828e-01, 8.51530e-02, 5.85950e-02, 3.17870e-02, 1.97160e-02, 9.64300e-03, 5.37900e-03, 3.22300e-03, 1.13400e-03, 5.15000e-04]
    dict["phi5"] = [0.0202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.001 , 0.004 , 0.01  , 0.034 , 0.067 , 0.143 , 0.203, 0.277, 0.309 , 0.321 , 0.329 , 0.33  ]
    dict["phi6"] = [300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300]
    dict["tau1"] = [0.3894, 0.4   , 0.4   , 0.4026, 0.4063, 0.4095, 0.4124, 0.4179, 0.4219, 0.4244, 0.4275, 0.4292, 0.4313, 0.4341, 0.4363, 0.4396, 0.4419, 0.4459, 0.4484, 0.4515, 0.4534, 0.4558, 0.4574, 0.4584, 0.4601, 0.4612]
    dict["tau2"] = [0.2578, 0.26  , 0.26  , 0.2637, 0.2689, 0.2736, 0.2777, 0.2855, 0.2913, 0.2949, 0.2993, 0.3017, 0.3047, 0.3087, 0.3119, 0.3165, 0.3199, 0.3255, 0.3291, 0.3335, 0.3363, 0.3398, 0.3419, 0.3435, 0.3459, 0.3474]
    dict["sigma1"] = [0.4785, 0.4912, 0.4912, 0.4904, 0.4988, 0.5049, 0.5096, 0.5179, 0.5236, 0.527 , 0.5308, 0.5328, 0.5351, 0.5377, 0.5395, 0.5422, 0.5433, 0.5294, 0.5105, 0.4783, 0.4681, 0.4617, 0.4571, 0.4535, 0.4471, 0.4426]
    dict["sigma2"] = [0.3629, 0.3762, 0.3762, 0.3762, 0.3849, 0.391 , 0.3957, 0.4043, 0.4104, 0.4143, 0.4191, 0.4217, 0.4252, 0.4299, 0.4338, 0.4399, 0.4446, 0.4533, 0.4594, 0.468 , 0.4681, 0.4617, 0.4571, 0.4535, 0.4471, 0.4426]
    dict["sigma3"] = [0.7504, 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.7999, 0.7997, 0.7988, 0.7966, 0.7792, 0.7504, 0.7136, 0.7035, 0.7006, 0.7001, 0.7   , 0.7   , 0.7   ]
    dict["sigma2_JP"] = [0.3918, 0.4528, 0.4528, 0.4551, 0.4571, 0.4642, 0.4716, 0.5022, 0.523 , 0.5278, 0.5304, 0.531 , 0.5312, 0.5309, 0.5307, 0.531 , 0.5313, 0.5309, 0.5302, 0.5276, 0.5167, 0.4917, 0.4682, 0.4517, 0.4167, 0.3755]
    dict["gamma_JP_IT"] = [2.2306, 1.5817, 1.5817, 1.574 , 1.5544, 1.5502, 1.5391, 1.4804, 1.4094, 1.3682, 1.3241, 1.3071, 1.2931, 1.315 , 1.3514, 1.4051, 1.4402, 1.528 , 1.6523, 1.8872, 2.1348, 3.5752, 3.8646, 3.7292, 2.3763, 1.7679]
    dict["gamma_Wn"] = [0.335 , 0.7594, 0.7594, 0.7606, 0.7642, 0.7676, 0.7739, 0.7956, 0.7932, 0.7768, 0.7437, 0.7219, 0.6922, 0.6579, 0.6362, 0.6049, 0.5507, 0.3582, 0.2003, 0.0356, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ]
    dict["phi1_JP"] = [-0.7966, -0.6846, -0.6846, -0.6681, -0.6314, -0.5855, -0.5457, -0.4685, -0.4985, -0.5603, -0.6451, -0.6981, -0.7653, -0.8469, -0.8999, -0.9618, -0.9945, -1.0225, -1.0002, -0.9245, -0.8626, -0.7882, -0.7195, -0.656 , -0.5202, -0.4068]
    dict["phi5_JP"] = [0.9488, 0.459 , 0.459 , 0.458 , 0.462 , 0.453 , 0.436 , 0.383 , 0.375 , 0.377 , 0.379 , 0.38  , 0.384 , 0.393 , 0.408 , 0.462 , 0.524 , 0.658 , 0.78  , 0.96  , 1.11  , 1.291 , 1.387 , 1.433 , 1.46  , 1.464 ]
    dict["phi6_JP"] = [800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800]
    
    
    if region == 2:
        dict["sigma2"] = dict["sigma2_JP"]
        dict["phi1"] = dict["phi1_JP"]
        dict["phi5"] = dict["phi5_JP"]
        phi6 = phi6_JP
    
#     fmag
    term6 = dict["c2"][ip]*(M-6)
    term7 = (dict["c2"][ip]-dict["c3"][ip])/dict["c_n"][ip]*np.log(1+np.exp(dict["c_n"][ip]*(dict["c_m"][ip]-M)))

#     Distance Scaling and attenuation term
    term8 = dict["c4"][ip] * np.log(R_RUP + dict["c5"][ip] * np.cosh(dict["c6"][ip] * max(M-dict["c_HM"][ip],0)))
    term9 = (dict["c4_a"][ip] - dict["c4"][ip]) * np.log(np.sqrt(R_RUP**2 + dict["c_RB"][ip]**2))
    term10 = (dict["c_g1"][ip] + dict["c_g2"][ip] / (np.cosh(max(M-dict["c_g3"][ip],0)))) * R_RUP

    if (region == 2 or region == 4):
        if (M>6 and M<6.9):
            term10 = dict["gamma_JP_IT"][ip] * term10

    if region == 3:
        term10 = dict["gamma_Wn"][ip] * term10

#     Style of faulting term
    term2=(dict["c1_a"][ip] + dict["c1_c"][ip] / np.cosh(2 * max(M-4.5,0))) * F_RV
    term3=(dict["c1_b"][ip] + dict["c1_d"][ip] / np.cosh(2 * max(M-4.5,0))) * F_NM

#   Ztor term
    if F_RV == 1:
        E_Ztor = (max(2.704 - 1.226*max(M-5.849,0), 0)) ** 2
    else:
        E_Ztor = (max(2.673 - 1.136*max(M-4.970,0), 0)) ** 2

    if Ztor == 999:
        Ztor = E_Ztor

    delta_ZTOR = Ztor - E_Ztor
    
    term4 = (dict["c7"][ip] + dict["c7_b"][ip] / np.cosh(2 * max(M-4.5,0))) * delta_ZTOR
    
#     Hanging wall term
    term12 = dict["c9"][ip] * HW * np.cos(delta) * (dict["c9_a"][ip] + (1 - dict["c9_a"][ip]) * np.tanh(Rx/dict["c9_b"][ip])) * (1 - np.sqrt(R_JB**2 + Ztor**2) / (R_RUP + 1))
    
#     Basin Depth term
#   Z1.0 (m) ~ Vs30 (m/s) relationship

    if region != 2: # in California and non-Japan region
        z_1 = np.exp(-7.15 / 4 * np.log((Vs30**4 + 570.94**4) / (1360**4 + 570.94**4)))
    else:
        z_1 = np.exp(-5.23 / 2 * np.log((Vs30**2 + 412.39**2) / (1360**2 + 412.39**2)))

    if Z10 ==999:
        d_Z1 = 0
    else:
        d_Z1 = Z10 * 1000 - z_1
    
#     Dip term
    term5 = (dict["c11"][ip] + dict["c11_b"][ip] / np.cosh(2 * max(M-4.5,0))) * (np.cos(delta)**2)
    
#     Directivity
    term11 = dict["c8"][ip] * max(1 - max(R_RUP-40,0)/30,0) * min(max(M-5.5,0)/0.8,1) * np.exp(-dict["c8_a"][ip] * (M-dict["c8_b"][ip])**2) * d_DPP
    term1 = dict["c1"][ip]
    ln_yrefij = term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9 + term10 + term11 + term12
    yrefij = np.exp(ln_yrefij)
    
#     Site response
    term14 = dict["phi1"][ip] * min(np.log(Vs30/1130),0)
    term15 = dict["phi2"][ip] * (np.exp(dict["phi3"][ip]*(min(Vs30,1130)-360)) - np.exp(dict["phi3"][ip] * (1130-360))) * np.log((yrefij + dict["phi4"][ip]) / dict["phi4"][ip])
    term16 = dict["phi5"][ip] * (1 - np.exp(-d_Z1 / dict["phi6"][ip]))

    Sa = yrefij * np.exp(term14 + term15 + term16)
    
#     Compute standard deviation
    Finferred = (FVS30==0) # equal 1 if Vs30 is inferred from geology.
    Fmeasured = (FVS30==1) # equal 1 if Vs30 is measured.


    NL0 = dict["phi2"][ip] * (np.exp(dict["phi3"][ip] * (min(Vs30,1130)-360)) - np.exp(dict["phi3"][ip] * (1130-360))) * (yrefij / (yrefij + dict["phi4"][ip]))
    
    sigmaNL0 = (dict["sigma1"][ip] + (dict["sigma2"][ip] - dict["sigma1"][ip]) / 1.5 * (min(max(M,5),6.5) - 5)) * np.sqrt((dict["sigma3"][ip] * Finferred + 0.7 * Fmeasured) + (1 + NL0)**2)

    tau = dict["tau1"][ip] + (dict["tau2"][ip] - dict["tau1"][ip]) / 1.5 * (min(max(M,5),6.5) - 5)
    sigma = np.sqrt((1 + NL0)**2 * tau**2 + sigmaNL0**2)
    
    return Sa, sigma


def gmpe_CY_2014(M, T, Rup, Rjb, Rx, Ztor, delta, lambda0, Z10, Vs30, Fhw, FVS30, region, *d_DPP):
    period = [-1, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.12, 0.15, 0.17, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 7.5, 10]

    delta = delta * np.pi / 180.0
    frv = (lambda0 >= 30 and lambda0 <= 150) # frv: 1 for lambda between 30 and 150, 0 otherwise
    fnm = (lambda0 >= -120 and lambda0 <= -60)  # fnm: 1 for lambda between -120 and -60, 0 otherwise

    if Fhw == 1:
        HW = 1
    elif Fhw == 0:
        HW = 0
    else:
        HW = Rx>=0
    
    arg = list(d_DPP)
    if len(arg) < 1: # for backward compatibility with prior versions that omit this term
        d_DPP = 0 # for median calculatio, d_DPP=0.

    if type(T) is not list:
        T = [T]
        
    if (len(T) == 1 and T == 1000): # Compute Sa and sigma with pre-defined period
        Sa = np.zeros((1,len(period)-2))
        sigma = np.zeros((1,len(period)-2))
        period1 = period[2:]
        for ip in range(2, len(period)): # 3:length(period)
            Sa[ip-2],sigma[ip-2] = CY_2014_sub(M, ip, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30,region, d_DPP)
            [PGA,sigma_PGA] = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30,region, d_DPP)
            if (Sa[ip-2]<PGA and period1[ip-2]<=0.3):
                Sa[ip-2]= PGA
    else:                            # Compute Sa and sigma with user-defined period
        Sa = np.zeros((1, len(T)))
        sigma = np.zeros((1, len(T)))
        period1 = T
        for i in range(len(T)):
            Ti = T[i]
            if len(np.where(abs(np.array(period)-Ti) < 0.0001)[0]) == 0: # The user defined period requires interpolation
                T_low = max(np.array(period)[np.array(period) < Ti])
                T_high = min(np.array(period)[np.array(period) > Ti])
                ip_low  = period.index(T_low)
                ip_high = period.index(T_high)

                Sa_low, sigma_low = CY_2014_sub(M, ip_low, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                Sa_high, sigma_high = CY_2014_sub(M, ip_high, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                PGA, sigma_PGA = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)

                x = [np.log(T_low), np.log(T_high)]
                Y_sa = [np.log(Sa_low), np.log(Sa_high)]
                Y_sigma = [sigma_low, sigma_high]
                interpF1 = interp1d(x, Y_sa)
                Sa[i] = np.exp(interpF1(np.log(Ti)))
                interpF2 = interp1d(x, Y_sigma)
                sigma[i] = interpF2(np.log(Ti))
                if (Sa[i] < PGA and period1[i] <= 0.3):
                    Sa[i]= PGA
            else:
                ip_T = np.where(abs(np.array(period) - Ti) < 0.0001)[0][0]
                Sa[i], sigma[i] = CY_2014_sub(M, ip_T, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                PGA, sigma_PGA = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                if (Sa[i]<PGA and period1[i] <= 0.3):
                    Sa[i] = PGA
    
    return Sa[0][0], sigma[0][0], period1
# Note: The codes were originally created by Prof. <NAME> in the MATLAB

import numpy as np
from scipy.stats import norm
from scipy.interpolate import interp1d

# CY_14 attenuation equation - Update of the Chiou and Youngs NGA Model for
# the Average Horizontal Component of Peak Ground Motion and Response
# Spectra, Chiou and Youngs (2014)


##############
### Inputs ###
##############

# M     = Moment Magnitude
# T     = Period (sec); Use Period = -1 for PGV computation
#                 Use 1000 for output the array of Sa with original period
#                 (no interpolation)
# Rup   = Closest distance (km) to the ruptured plane
# Rjb   = Joyner-Boore distance (km); closest distance (km) to surface
#       projection of rupture plane
# Rx    =    Horizontal distance from top of rupture measured perpendicular 
#       to fault strike (km). See Figures a, b and c for illustation% Ztor  = Depth(km) to the top of ruptured plane
# delta = Fault dip angle (in degree)
# lamda = Rake angle      (in degree)
# region        = 0 for global (incl. Taiwan)
#               = 1 for California
#               = 2 for Japan
#               = 3 for China 
#               = 4 for Italy 
#               = 5 for Turkey
# Z10            = Basin depth (km); depth from the groundsurface to the
#                   1km/s shear-wave horizon.
#               = 999 if unknown
# Vs30          = shear wave velocity averaged over top 30 m in m/s
#               = ref: 1130
# FVS30         = 1 for Vs30 is inferred from geology
#               = 0 for measured  Vs30
# d_DPP         = DPP centered on the site and earthquake specific average
#               DPP, = 0 for median calc
#               %d_DPP functionality added BAB 26 July 2020

###############
### Outputs ###
###############

# Sa: Median spectral acceleration prediction
# sigma: logarithmic standard deviation of spectral acceleration
#          prediction


def CY_2014_sub(M, ip, R_RUP, R_JB, Rx, Ztor, delta, F_RV, F_NM, HW, Z10, Vs30, FVS30, region, d_DPP):
    
    dict = {}
    dict["c2"] = [1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06]
    dict["c4"] = [-2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1]
    dict["c4_a"] = [-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5]
    dict["c_RB"] = [50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50]
    dict["c8"] = [0.2154, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.0991, 0.1982, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154]
    dict["c8_a"] = [0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695]
    dict["c1"] = [ 2.3549, -1.5065, -1.5065, -1.4798, -1.2972, -1.1007, -0.9292, -0.658 , -0.5613, -0.5342, -0.5462, -0.5858, -0.6798, -0.8663, -1.0514, -1.3794, -1.6508, -2.1511, -2.5365, -3.0686, -3.4148, -3.9013, -4.2466, -4.5143, -5.0009, -5.3461]
    dict["c1_a"] = [0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165, 0.165, 0.165, 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.1645, 0.1168, 0.0732, 0.0484, 0.022 , 0.0124]
    dict["c1_b"] = [-0.0626, -0.255 , -0.255 , -0.255 , -0.255 , -0.255 , -0.255 , -0.254 , -0.253 , -0.252 , -0.25  , -0.248 , -0.2449, -0.2382, -0.2313, -0.2146, -0.1972, -0.162 , -0.14  , -0.1184, -0.11  , -0.104 , -0.102 , -0.101 , -0.101 , -0.1   ]
    dict["c1_c"] = [-0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.1645, -0.1168, -0.0732, -0.0484, -0.022 , -0.0124]
    dict["c1_d"] = [0.0626, 0.255 , 0.255 , 0.255 , 0.255 , 0.255 , 0.255 , 0.254 , 0.253 , 0.252 , 0.25  , 0.248 , 0.2449, 0.2382, 0.2313, 0.2146, 0.1972, 0.162 , 0.14  , 0.1184, 0.11  , 0.104 , 0.102 , 0.101 , 0.101 , 0.1   ]
    dict["c_n"] = [ 3.3024, 16.0875, 16.0875, 15.7118, 15.8819, 16.4556, 17.6453, 20.1772, 19.9992, 18.7106, 16.6246, 15.3709, 13.7012, 11.2667, 9.1908,  6.5459,  5.2305,  3.7896,  3.3024,  2.8498,  2.5417, 2.1488,  1.8957,  1.7228,  1.5737,  1.5265]
    dict["c_m"] = [5.423 , 4.9993, 4.9993, 4.9993, 4.9993, 4.9993, 4.9993, 5.0031, 5.0172, 5.0315, 5.0547, 5.0704, 5.0939, 5.1315, 5.167 , 5.2317, 5.2893, 5.4109, 5.5106, 5.6705, 5.7981, 5.9983, 6.1552, 6.2856, 6.5428, 6.7415]
    dict["c3"] = [2.3152, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9795, 2.0362, 2.0823, 2.1521, 2.2574, 2.344 , 2.4709, 2.5567, 2.6812, 2.7474, 2.8161, 2.8514, 2.8875, 2.9058, 2.9169, 2.932 , 2.9396]
    dict["c5"] = [5.8096, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.8305, 7.1333, 7.3621, 7.4365, 7.4972, 7.5416, 7.56  , 7.5735, 7.5778, 7.5808, 7.5814, 7.5817, 7.5818, 7.5818, 7.5818, 7.5818, 7.5818, 7.5818]
    dict["c_HM"] = [3.0514, 3.0956, 3.0956, 3.0963, 3.0974, 3.0988, 3.1011, 3.1094, 3.2381, 3.3407, 3.43  , 3.4688, 3.5146, 3.5746, 3.6232, 3.6945, 3.7401, 3.7941, 3.8144, 3.8284, 3.833 , 3.8361, 3.8369, 3.8376, 3.838 , 3.838 ]
    dict["c6"] = [0.4407, 0.4908, 0.4908, 0.4925, 0.4992, 0.5037, 0.5048, 0.5048, 0.5048, 0.5048, 0.5045, 0.5036, 0.5016, 0.4971, 0.4919, 0.4807, 0.4707, 0.4575, 0.4522, 0.4501, 0.45  , 0.45  , 0.45  , 0.45  , 0.45  , 0.45  ]
    dict["c7"] = [0.0324, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.016 , 0.0062, 0.0029, 0.0007, 0.0003]
    dict["c7_b"] = [ 0.0097,  0.0462,  0.0462,  0.0472,  0.0533,  0.0596,  0.0639, 0.063 ,  0.0532,  0.0452,  0.0345,  0.0283, 0.0202,  0.009 , -0.0004, -0.0155, -0.0278, -0.0477, -0.0559, -0.063 , -0.0665, -0.0516, -0.0448, -0.0424, -0.0348, -0.0253]
    dict["c8_b"] = [5.    , 0.4833, 0.4833, 1.2144, 1.6421, 1.9456, 2.181 , 2.6087, 2.9122, 3.1045, 3.3399, 3.4719, 3.6434, 3.8787, 4.0711, 4.3745, 4.6099, 5.0376, 5.3411, 5.7688, 6.0723, 6.5   , 6.8035, 7.0389, 7.4666, 7.77  ]
    dict["c9"] = [0.3079, 0.9228, 0.9228, 0.9296, 0.9396, 0.9661, 0.9794, 1.026 , 1.0177, 1.0008, 0.9801, 0.9652, 0.9459, 0.9196, 0.8829, 0.8302, 0.7884, 0.6754, 0.6196, 0.5101, 0.3917, 0.1244, 0.0086, 0.    , 0.    , 0.    ]
    dict["c9_a"] = [0.1   , 0.1202, 0.1202, 0.1217, 0.1194, 0.1166, 0.1176, 0.1171, 0.1146, 0.1128, 0.1106, 0.115 , 0.1208, 0.1208, 0.1175, 0.106 , 0.1061, 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   ]
    dict["c9_b"] = [6.5   , 6.8607, 6.8607, 6.8697, 6.9113, 7.0271, 7.0959, 7.3298, 7.2588, 7.2372, 7.2109, 7.2491, 7.2988, 7.3691, 6.8789, 6.5334, 6.526 , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   ]
    dict["c11"] = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
    dict["c11_b"] = [-0.3834, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.444 , -0.3539, -0.2688, -0.1793, -0.1428, -0.1138, -0.1062, -0.102 , -0.1009, -0.1003, -0.1001, -0.1001, -0.1, -0.1]
    dict["c_g1"] = [-0.001852, -0.007146, -0.007146, -0.007249, -0.007869, -0.008316, -0.008743, -0.009537, -0.00983 , -0.009913, -0.009896, -0.009787, -0.009505, -0.008918, -0.008251, -0.007267, -0.006492, -0.005147, -0.004277, -0.002979, -0.002301, -0.001344, -0.001084, -0.00101 , -0.000964, -0.00095 ]
    dict["c_g2"] = [-0.007403, -0.006758, -0.006758, -0.006758, -0.006758, -0.006758, -0.006758, -0.00619 , -0.005332, -0.004732, -0.003806, -0.00328 , -0.00269 , -0.002128, -0.001812, -0.001274, -0.001074, -0.001115, -0.001197, -0.001675, -0.002349, -0.003306, -0.003566, -0.00364 , -0.003686, -0.0037  ]
    dict["c_g3"] = [4.3439, 4.2542, 4.2542, 4.2386, 4.2519, 4.296 , 4.3578, 4.5455, 4.7603, 4.8963, 5.0644, 5.1371, 5.188 , 5.2164, 5.1954, 5.0899, 4.7854, 4.3304, 4.1667, 4.0029, 3.8949, 3.7928, 3.7443, 3.709 , 3.6632, 3.623 ]
    dict["phi1"] = [-0.7936, -0.521 , -0.521 , -0.5055, -0.4368, -0.3752, -0.3469, -0.3747, -0.444 , -0.4895, -0.5477, -0.5922, -0.6693, -0.7766, -0.8501, -0.9431, -1.0044, -1.0602, -1.0941, -1.1142, -1.1154, -1.1081, -1.0603, -0.9872, -0.8274, -0.7053]
    dict["phi2"] = [-0.0699, -0.1417, -0.1417, -0.1364, -0.1403, -0.1591, -0.1862, -0.2538, -0.2943, -0.3077, -0.3113, -0.3062, -0.2927, -0.2662, -0.2405, -0.1975, -0.1633, -0.1028, -0.0699, -0.0425, -0.0302, -0.0129, -0.0016,  0.    ,  0.    ,  0.    ]
    dict["phi3"] = [-0.008444, -0.00701 , -0.00701 , -0.007279, -0.007354, -0.006977, -0.006467, -0.005734, -0.005604, -0.005696, -0.005845, -0.005959, -0.006141, -0.006439, -0.006704, -0.007125, -0.007435, -0.00812 , -0.008444, -0.007707, -0.004792, -0.001828, -0.001523, -0.00144 , -0.001369, -0.001361]
    dict["phi4"] = [5.41000e+00, 1.02151e-01, 1.02151e-01, 1.08360e-01, 1.19888e-01, 1.33641e-01, 1.48927e-01, 1.90596e-01, 2.30662e-01, 2.53169e-01, 2.66468e-01, 2.65060e-01, 2.55253e-01, 2.31541e-01, 2.07277e-01, 1.65464e-01, 1.33828e-01, 8.51530e-02, 5.85950e-02, 3.17870e-02, 1.97160e-02, 9.64300e-03, 5.37900e-03, 3.22300e-03, 1.13400e-03, 5.15000e-04]
    dict["phi5"] = [0.0202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.001 , 0.004 , 0.01  , 0.034 , 0.067 , 0.143 , 0.203, 0.277, 0.309 , 0.321 , 0.329 , 0.33  ]
    dict["phi6"] = [300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300]
    dict["tau1"] = [0.3894, 0.4   , 0.4   , 0.4026, 0.4063, 0.4095, 0.4124, 0.4179, 0.4219, 0.4244, 0.4275, 0.4292, 0.4313, 0.4341, 0.4363, 0.4396, 0.4419, 0.4459, 0.4484, 0.4515, 0.4534, 0.4558, 0.4574, 0.4584, 0.4601, 0.4612]
    dict["tau2"] = [0.2578, 0.26  , 0.26  , 0.2637, 0.2689, 0.2736, 0.2777, 0.2855, 0.2913, 0.2949, 0.2993, 0.3017, 0.3047, 0.3087, 0.3119, 0.3165, 0.3199, 0.3255, 0.3291, 0.3335, 0.3363, 0.3398, 0.3419, 0.3435, 0.3459, 0.3474]
    dict["sigma1"] = [0.4785, 0.4912, 0.4912, 0.4904, 0.4988, 0.5049, 0.5096, 0.5179, 0.5236, 0.527 , 0.5308, 0.5328, 0.5351, 0.5377, 0.5395, 0.5422, 0.5433, 0.5294, 0.5105, 0.4783, 0.4681, 0.4617, 0.4571, 0.4535, 0.4471, 0.4426]
    dict["sigma2"] = [0.3629, 0.3762, 0.3762, 0.3762, 0.3849, 0.391 , 0.3957, 0.4043, 0.4104, 0.4143, 0.4191, 0.4217, 0.4252, 0.4299, 0.4338, 0.4399, 0.4446, 0.4533, 0.4594, 0.468 , 0.4681, 0.4617, 0.4571, 0.4535, 0.4471, 0.4426]
    dict["sigma3"] = [0.7504, 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.7999, 0.7997, 0.7988, 0.7966, 0.7792, 0.7504, 0.7136, 0.7035, 0.7006, 0.7001, 0.7   , 0.7   , 0.7   ]
    dict["sigma2_JP"] = [0.3918, 0.4528, 0.4528, 0.4551, 0.4571, 0.4642, 0.4716, 0.5022, 0.523 , 0.5278, 0.5304, 0.531 , 0.5312, 0.5309, 0.5307, 0.531 , 0.5313, 0.5309, 0.5302, 0.5276, 0.5167, 0.4917, 0.4682, 0.4517, 0.4167, 0.3755]
    dict["gamma_JP_IT"] = [2.2306, 1.5817, 1.5817, 1.574 , 1.5544, 1.5502, 1.5391, 1.4804, 1.4094, 1.3682, 1.3241, 1.3071, 1.2931, 1.315 , 1.3514, 1.4051, 1.4402, 1.528 , 1.6523, 1.8872, 2.1348, 3.5752, 3.8646, 3.7292, 2.3763, 1.7679]
    dict["gamma_Wn"] = [0.335 , 0.7594, 0.7594, 0.7606, 0.7642, 0.7676, 0.7739, 0.7956, 0.7932, 0.7768, 0.7437, 0.7219, 0.6922, 0.6579, 0.6362, 0.6049, 0.5507, 0.3582, 0.2003, 0.0356, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ]
    dict["phi1_JP"] = [-0.7966, -0.6846, -0.6846, -0.6681, -0.6314, -0.5855, -0.5457, -0.4685, -0.4985, -0.5603, -0.6451, -0.6981, -0.7653, -0.8469, -0.8999, -0.9618, -0.9945, -1.0225, -1.0002, -0.9245, -0.8626, -0.7882, -0.7195, -0.656 , -0.5202, -0.4068]
    dict["phi5_JP"] = [0.9488, 0.459 , 0.459 , 0.458 , 0.462 , 0.453 , 0.436 , 0.383 , 0.375 , 0.377 , 0.379 , 0.38  , 0.384 , 0.393 , 0.408 , 0.462 , 0.524 , 0.658 , 0.78  , 0.96  , 1.11  , 1.291 , 1.387 , 1.433 , 1.46  , 1.464 ]
    dict["phi6_JP"] = [800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800]
    
    
    if region == 2:
        dict["sigma2"] = dict["sigma2_JP"]
        dict["phi1"] = dict["phi1_JP"]
        dict["phi5"] = dict["phi5_JP"]
        phi6 = phi6_JP
    
#     fmag
    term6 = dict["c2"][ip]*(M-6)
    term7 = (dict["c2"][ip]-dict["c3"][ip])/dict["c_n"][ip]*np.log(1+np.exp(dict["c_n"][ip]*(dict["c_m"][ip]-M)))

#     Distance Scaling and attenuation term
    term8 = dict["c4"][ip] * np.log(R_RUP + dict["c5"][ip] * np.cosh(dict["c6"][ip] * max(M-dict["c_HM"][ip],0)))
    term9 = (dict["c4_a"][ip] - dict["c4"][ip]) * np.log(np.sqrt(R_RUP**2 + dict["c_RB"][ip]**2))
    term10 = (dict["c_g1"][ip] + dict["c_g2"][ip] / (np.cosh(max(M-dict["c_g3"][ip],0)))) * R_RUP

    if (region == 2 or region == 4):
        if (M>6 and M<6.9):
            term10 = dict["gamma_JP_IT"][ip] * term10

    if region == 3:
        term10 = dict["gamma_Wn"][ip] * term10

#     Style of faulting term
    term2=(dict["c1_a"][ip] + dict["c1_c"][ip] / np.cosh(2 * max(M-4.5,0))) * F_RV
    term3=(dict["c1_b"][ip] + dict["c1_d"][ip] / np.cosh(2 * max(M-4.5,0))) * F_NM

#   Ztor term
    if F_RV == 1:
        E_Ztor = (max(2.704 - 1.226*max(M-5.849,0), 0)) ** 2
    else:
        E_Ztor = (max(2.673 - 1.136*max(M-4.970,0), 0)) ** 2

    if Ztor == 999:
        Ztor = E_Ztor

    delta_ZTOR = Ztor - E_Ztor
    
    term4 = (dict["c7"][ip] + dict["c7_b"][ip] / np.cosh(2 * max(M-4.5,0))) * delta_ZTOR
    
#     Hanging wall term
    term12 = dict["c9"][ip] * HW * np.cos(delta) * (dict["c9_a"][ip] + (1 - dict["c9_a"][ip]) * np.tanh(Rx/dict["c9_b"][ip])) * (1 - np.sqrt(R_JB**2 + Ztor**2) / (R_RUP + 1))
    
#     Basin Depth term
#   Z1.0 (m) ~ Vs30 (m/s) relationship

    if region != 2: # in California and non-Japan region
        z_1 = np.exp(-7.15 / 4 * np.log((Vs30**4 + 570.94**4) / (1360**4 + 570.94**4)))
    else:
        z_1 = np.exp(-5.23 / 2 * np.log((Vs30**2 + 412.39**2) / (1360**2 + 412.39**2)))

    if Z10 ==999:
        d_Z1 = 0
    else:
        d_Z1 = Z10 * 1000 - z_1
    
#     Dip term
    term5 = (dict["c11"][ip] + dict["c11_b"][ip] / np.cosh(2 * max(M-4.5,0))) * (np.cos(delta)**2)
    
#     Directivity
    term11 = dict["c8"][ip] * max(1 - max(R_RUP-40,0)/30,0) * min(max(M-5.5,0)/0.8,1) * np.exp(-dict["c8_a"][ip] * (M-dict["c8_b"][ip])**2) * d_DPP
    term1 = dict["c1"][ip]
    ln_yrefij = term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9 + term10 + term11 + term12
    yrefij = np.exp(ln_yrefij)
    
#     Site response
    term14 = dict["phi1"][ip] * min(np.log(Vs30/1130),0)
    term15 = dict["phi2"][ip] * (np.exp(dict["phi3"][ip]*(min(Vs30,1130)-360)) - np.exp(dict["phi3"][ip] * (1130-360))) * np.log((yrefij + dict["phi4"][ip]) / dict["phi4"][ip])
    term16 = dict["phi5"][ip] * (1 - np.exp(-d_Z1 / dict["phi6"][ip]))

    Sa = yrefij * np.exp(term14 + term15 + term16)
    
#     Compute standard deviation
    Finferred = (FVS30==0) # equal 1 if Vs30 is inferred from geology.
    Fmeasured = (FVS30==1) # equal 1 if Vs30 is measured.


    NL0 = dict["phi2"][ip] * (np.exp(dict["phi3"][ip] * (min(Vs30,1130)-360)) - np.exp(dict["phi3"][ip] * (1130-360))) * (yrefij / (yrefij + dict["phi4"][ip]))
    
    sigmaNL0 = (dict["sigma1"][ip] + (dict["sigma2"][ip] - dict["sigma1"][ip]) / 1.5 * (min(max(M,5),6.5) - 5)) * np.sqrt((dict["sigma3"][ip] * Finferred + 0.7 * Fmeasured) + (1 + NL0)**2)

    tau = dict["tau1"][ip] + (dict["tau2"][ip] - dict["tau1"][ip]) / 1.5 * (min(max(M,5),6.5) - 5)
    sigma = np.sqrt((1 + NL0)**2 * tau**2 + sigmaNL0**2)
    
    return Sa, sigma


def gmpe_CY_2014(M, T, Rup, Rjb, Rx, Ztor, delta, lambda0, Z10, Vs30, Fhw, FVS30, region, *d_DPP):
    period = [-1, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.12, 0.15, 0.17, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 7.5, 10]

    delta = delta * np.pi / 180.0
    frv = (lambda0 >= 30 and lambda0 <= 150) # frv: 1 for lambda between 30 and 150, 0 otherwise
    fnm = (lambda0 >= -120 and lambda0 <= -60)  # fnm: 1 for lambda between -120 and -60, 0 otherwise

    if Fhw == 1:
        HW = 1
    elif Fhw == 0:
        HW = 0
    else:
        HW = Rx>=0
    
    arg = list(d_DPP)
    if len(arg) < 1: # for backward compatibility with prior versions that omit this term
        d_DPP = 0 # for median calculatio, d_DPP=0.

    if type(T) is not list:
        T = [T]
        
    if (len(T) == 1 and T == 1000): # Compute Sa and sigma with pre-defined period
        Sa = np.zeros((1,len(period)-2))
        sigma = np.zeros((1,len(period)-2))
        period1 = period[2:]
        for ip in range(2, len(period)): # 3:length(period)
            Sa[ip-2],sigma[ip-2] = CY_2014_sub(M, ip, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30,region, d_DPP)
            [PGA,sigma_PGA] = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30,region, d_DPP)
            if (Sa[ip-2]<PGA and period1[ip-2]<=0.3):
                Sa[ip-2]= PGA
    else:                            # Compute Sa and sigma with user-defined period
        Sa = np.zeros((1, len(T)))
        sigma = np.zeros((1, len(T)))
        period1 = T
        for i in range(len(T)):
            Ti = T[i]
            if len(np.where(abs(np.array(period)-Ti) < 0.0001)[0]) == 0: # The user defined period requires interpolation
                T_low = max(np.array(period)[np.array(period) < Ti])
                T_high = min(np.array(period)[np.array(period) > Ti])
                ip_low  = period.index(T_low)
                ip_high = period.index(T_high)

                Sa_low, sigma_low = CY_2014_sub(M, ip_low, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                Sa_high, sigma_high = CY_2014_sub(M, ip_high, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                PGA, sigma_PGA = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)

                x = [np.log(T_low), np.log(T_high)]
                Y_sa = [np.log(Sa_low), np.log(Sa_high)]
                Y_sigma = [sigma_low, sigma_high]
                interpF1 = interp1d(x, Y_sa)
                Sa[i] = np.exp(interpF1(np.log(Ti)))
                interpF2 = interp1d(x, Y_sigma)
                sigma[i] = interpF2(np.log(Ti))
                if (Sa[i] < PGA and period1[i] <= 0.3):
                    Sa[i]= PGA
            else:
                ip_T = np.where(abs(np.array(period) - Ti) < 0.0001)[0][0]
                Sa[i], sigma[i] = CY_2014_sub(M, ip_T, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                PGA, sigma_PGA = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                if (Sa[i]<PGA and period1[i] <= 0.3):
                    Sa[i] = PGA
    
    return Sa[0][0], sigma[0][0], period1
# Note: The codes were originally created by Prof. <NAME> in the MATLAB # CY_14 attenuation equation - Update of the Chiou and Youngs NGA Model for # the Average Horizontal Component of Peak Ground Motion and Response # Spectra, Chiou and Youngs (2014) ############## ### Inputs ### ############## # M     = Moment Magnitude # T     = Period (sec); Use Period = -1 for PGV computation #                 Use 1000 for output the array of Sa with original period #                 (no interpolation) # Rup   = Closest distance (km) to the ruptured plane # Rjb   = Joyner-Boore distance (km); closest distance (km) to surface #       projection of rupture plane # Rx    =    Horizontal distance from top of rupture measured perpendicular #       to fault strike (km). See Figures a, b and c for illustation% Ztor  = Depth(km) to the top of ruptured plane # delta = Fault dip angle (in degree) # lamda = Rake angle      (in degree) # region        = 0 for global (incl. Taiwan) #               = 1 for California #               = 2 for Japan #               = 3 for China #               = 4 for Italy #               = 5 for Turkey # Z10            = Basin depth (km); depth from the groundsurface to the #                   1km/s shear-wave horizon. #               = 999 if unknown # Vs30          = shear wave velocity averaged over top 30 m in m/s #               = ref: 1130 # FVS30         = 1 for Vs30 is inferred from geology #               = 0 for measured  Vs30 # d_DPP         = DPP centered on the site and earthquake specific average #               DPP, = 0 for median calc #               %d_DPP functionality added BAB 26 July 2020 ############### ### Outputs ### ############### # Sa: Median spectral acceleration prediction # sigma: logarithmic standard deviation of spectral acceleration #          prediction #     fmag #     Distance Scaling and attenuation term #     Style of faulting term #   Ztor term #     Hanging wall term #     Basin Depth term #   Z1.0 (m) ~ Vs30 (m/s) relationship # in California and non-Japan region #     Dip term #     Directivity #     Site response #     Compute standard deviation # equal 1 if Vs30 is inferred from geology. # equal 1 if Vs30 is measured. # frv: 1 for lambda between 30 and 150, 0 otherwise # fnm: 1 for lambda between -120 and -60, 0 otherwise # for backward compatibility with prior versions that omit this term # for median calculatio, d_DPP=0. # Compute Sa and sigma with pre-defined period # 3:length(period) # Compute Sa and sigma with user-defined period # The user defined period requires interpolation