import numpy as np

PERIODS = np.array([0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.250, 0.3, 0.4, 0.5, 0.6, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.5, 10.0])
VLIN = [865.1, 865.1, 1053.5, 1085.7, 1032.5, 877.6, 748.2, 654.3, 587.1, 503.0, 456.6, 430.3, 410.5, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0]
B = [-1.186, -1.186, -1.346, -1.471, -1.624, -1.931, -2.188, -2.381, -2.518, -2.657, -2.699, -2.599, -2.401, -1.955, -1.025, -0.299, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
T1 = [4.2203, 4.2203, 4.5371, 5.0733, 5.2892, 5.4563, 5.2684, 5.0594, 4.7945, 4.4644, 4.0181, 3.6055, 3.2174, 2.7981, 2.0123, 1.4128, 0.9976, 0.6443, 0.0657, -0.4624, -0.9809, -1.6017, -2.2937]
T2 = [-1.350, -1.350, -1.400, -1.450, -1.450, -1.450, -1.400, -1.350, -1.280, -1.180, -1.080, -0.990, -0.910, -0.850, -0.770, -0.710, -0.670, -0.640, -0.580, -0.540, -0.500, -0.460, -0.400]

T3 = 0.1
T4 = 0.9
T5 = 0.0

T6 = [-0.0012, -0.0012, -0.0012, -0.0012, -0.0012, -0.0014, -0.0018, -0.0023, -0.0027, -0.0035, -0.0044, -0.0050, -0.0058, -0.0062, -0.0064, -0.0064, -0.0064, -0.0064, -0.0064, -0.0064, -0.0064, -0.0064, -0.0064]
T7 = [1.0988, 1.0988, 1.2536, 1.4175, 1.3997, 1.3582, 1.1648, 0.9940, 0.8821, 0.7046, 0.5799, 0.5021, 0.3687, 0.1746, -0.0820, -0.2821, -0.4108, -0.4466, -0.4344, -0.4368, -0.4586, -0.4433, -0.4828]
T8 = [-1.42, -1.42, -1.65, -1.80, -1.80, -1.69, -1.49, -1.30, -1.18, -0.98, -0.82, -0.70, -0.54, -0.34, -0.05, 0.12, 0.25, 0.30, 0.30, 0.30, 0.30, 0.30, 0.30]

T9 = 0.4

T10 = [3.12, 3.12, 3.37, 3.37, 3.33, 3.25, 3.03, 2.80, 2.59, 2.20, 1.92, 1.70, 1.42, 1.10, 0.70, 0.70, 0.70, 0.70, 0.70, 0.70, 0.70, 0.70, 0.70]
T11 = [0.0130, 0.0130, 0.0130, 0.0130, 0.0130, 0.0130, 0.0129, 0.0129, 0.0128, 0.0127, 0.0125, 0.0124, 0.0120, 0.0114, 0.0100, 0.0085, 0.0069, 0.0054, 0.0027, 0.0005, -0.0013, -0.0033, -0.0060]
T12 = [0.980, 0.980, 1.288, 1.483, 1.613, 1.882, 2.076, 2.248, 2.348, 2.427, 2.399, 2.273, 1.993, 1.470, 0.408, -0.401, -0.723, -0.673, -0.627, -0.596, -0.566, -0.528, -0.504]
T13 = [-0.0135, -0.0135, -0.0138, -0.0142, -0.0145, -0.0153, -0.0162, -0.0172, -0.0183, -0.0206, -0.0231, -0.0256, -0.0296, -0.0363, -0.0493, -0.0610, -0.0711, -0.0798, -0.0935, -0.0980, -0.0980, -0.0980, -0.0980]
T14 = [-0.40, -0.40, -0.40, -0.40, -0.40, -0.40, -0.35, -0.31, -0.28, -0.23, -0.19, -0.16, -0.12, -0.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
T15 = [0.9969, 0.9969, 1.1030, 1.2732, 1.3042, 1.2600, 1.2230, 1.1600, 1.0500, 0.8000, 0.6620, 0.5800, 0.4800, 0.3300, 0.3100, 0.3000, 0.3000, 0.3000, 0.3000, 0.3000, 0.3000, 0.3000, 0.3000]
T16 = [-1.00, -1.00, -1.18, -1.36, -1.36, -1.30, -1.25, -1.17, -1.06, -0.78, -0.62, -0.50, -0.34, -0.14, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]

N = 1.18
C = 1.88

C1 = 7.8
C4 = 10
DC11 = -0.3
DC10 = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.143682921, 0.1, 0.073696559, 0.04150375, 0, -0.05849625, -0.1, -0.155033971, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2]
PHI = 0.60
TAU = 0.43
SIGMA = 0.74
SIGMASS = 0.60

def fmag(mag, pIdx, DC1):
    if mag <= C1 + DC1:
        return T4*(mag - (C1 + DC1)) + T13[pIdx]*((10 - mag)**2)
    else:
        return T5*(mag - (C1 + DC1)) + T13[pIdx]*((10 - mag)**2)

def fdepth(z, evt, pIdx):
    return T11[pIdx] * (z - 60) * evt

def ffaba(faba, evt, r, pIdx):
    if evt == 1:
        return (T7[pIdx] + T8[pIdx] * np.log(max(r, 85)/40.0))*faba
    else:
        return (T15[pIdx] + T16[pIdx] * np.log(max(r, 100)/40.0))*faba

def fsite(mag, r, evt, z, faba, vs30, pIdx):
    vlin = VLIN[pIdx]
    vs = 1000 if vs30 >= 1000 else vs30

    if vs30 >= vlin:
        return (T12[pIdx]  + B[pIdx] * N) * np.log(vs/vlin)
    else:
        pga1000 = np.exp(computeSpectra(mag, r, evt, z, faba, 1000.0, 0.01))
        return T12[pIdx] * np.log(vs/vlin) - B[pIdx] * np.log(pga1000 + C) + B[pIdx] * np.log(pga1000 + C*((vs/vlin)**N))

def computeSpectra(mag, r, evt, z, faba, vs30, per):

    pIdx = np.nonzero(PERIODS == per)[0][0]
    DC1 = 0
    if evt == 1:
        DC1 = DC11
    else:
        DC1 = DC10[pIdx]
    # Find out the value of dC1
    lnSa = T1[pIdx] + T4*DC1 + (T2[pIdx] + T14[pIdx]*evt + T3*(mag - 7.8)) * np.log(r + C4 * np.exp((mag - 6)*T9)) + T6[pIdx]*r + T10[pIdx]*evt + fmag(mag, pIdx, DC1) + fdepth(z, evt, pIdx) + ffaba(faba, evt, r, pIdx) + fsite(mag, r, evt, z, faba, vs30, pIdx)
    #print 'For Vs30 = ', vs30
    #print 'fmag', fmag(mag, pIdx)
    #print 'fdepth', fdepth(z, evt, pIdx)
    #print 'ffaba', ffaba(faba, evt, r, pIdx)
    #print 'fsite', fsite(mag, r, evt, z, faba, vs30, pIdx), 'Vs', vs30
    #print 'otherTerms', T1[pIdx] + T4*DC1 + (T2[pIdx] + T14[pIdx]*evt + T3*(mag - 7.8)) * np.log(r + C4 * np.exp((mag - 6)*T9)) + T6[pIdx]*r + T10[pIdx]*evt
    #print 'T1', T1[pIdx] + T4*DC1
    #print 'T2', (T2[pIdx] + T14[pIdx]*evt + T3*(mag - 7.8)) * np.log(r + C4 * np.exp((mag - 6)*T9))
    #print 'T3', T6[pIdx]*r
    #print 'lnSa', lnSa

    return lnSa

def interEventSigma():
    return TAU

def intraEventSigma():
    return PHI
