R/zea_2016_gmpe.R

In wltcwpf/GMPE: A pacakge of commonly used Ground Motion Prediction Equations (GMPE)

#' The GMPE of Zea 2016 for shallow crustal or upper mantle event
#'
#' This function calculates the ground motion median values and standard deviations
#' @param mag Moment magnitude, a numeric value
#' @param T Period (sec) (range should be between 0.005 s and 5 s).
#' Use 1000 (by default) for output the array of Sa with original Zhao et al per-defined periods
#' @param fDepth Depth(km) to the top of ruptured plane (km)
#' @param dist The closest distance (km) to rupture plane when is available,
#' otherwise hypocentral distance.
#' @param DistV The distance of volcanic path (km)
#' @param eqkType The indicator of earthquake type: 1-shallow crustal; 2-upper mantel;
#' only two values are valid
#' @param FaultMech The indicator of fault mechanism: 1-Reverse; 2-Strike Slip; 3-Normal
#' @param SiteClass The site class: 1 for SC I, 2 for SC II, 3 for SC III 4 for SC IV sites. The relationship
#' between site class and Vs30 is given by Zhao et al., 2016 (https://doi.org/10.1785/0120150056).
#' SC I: Vs30 > 600 m/s; SC II: 300 < Vs30 =< 600 m/s; SC III: 200 < Vs30 =< 300 m/s; SC IV: Vs30 <= 200 m/s
#' @return A list of six elements is returned: Sa - median spectral acceleration prediction (in g);
#' RockSa - median spectral acceleration prediction on Rock reference site (in g);
#' SiteAmp - Site amplication response spectral (subtraction in log(g)); sigma - totla standard deviation (log);
#' phiSS - single station standard deviation (log); period - the corresponding oscillator periods
#' @examples zea_2017_SC_UM(mag = 7.5, T = 1000, fDepth = 0.0, dist = 40.0, DistV = 0.0,
#' eqkType = 2, FaultMech = 2, SiteClass = 4)
#'
#' zea_2017_SC_UM(mag = 7.5, T = c(0.025, 0.5), fDepth = 0.0, dist = 40.0,
#' DistV = 0.0, eqkType = 1, FaultMech = 2, SiteClass = 4)
#' @references Zhao, et al. (2017). Ground‐Motion Prediction Equations for the Vertical Component
#' from Shallow Crustal and Upper‐Mantle Earthquakes in Japan Using Site Classes as the Site Term.
#' Bulletin of the Seismological Society of America. 107(5), 2310-2327.
#' @export
#' @importFrom stats approx
zea_2017_SC_UM <- function(mag, T = 1000, fDepth, dist, DistV,
                                 eqkType, FaultMech, SiteClass){

  period <- c(0.005,0.01,0.02,0.03,0.04,0.05,0.06,
              0.07,0.08,0.09,0.1,0.12,0.14,0.15,0.16,0.18,0.2,0.25,0.3,
              0.35,0.4,0.45,0.5,0.6,0.7,0.8,0.9,1.0,1.25,1.5,2.0,2.5,
              3.0,3.5,4.0,4.5,5.0)

  g <- 980.7

  Numperiod <- length(period)

  SpeCalLog <- rep(0, Numperiod)

  RockSpeCalLog <- rep(0, Numperiod)

  SiteAmpCalLog <- rep(0, Numperiod)

  SigmaT <- rep(0, Numperiod)

  sigmaSS <- rep(0, Numperiod)

  # model coefficients
  alfa <- c(-3.224,-3.357,-3.552,-3.640,-3.758,-3.826,
            -3.890,-3.965,-4.055,-4.153,-4.255,-4.466,-4.677,-4.781,
            -4.883,-5.085,-5.233,-5.229,-5.226,-5.223,-5.221,-5.218,
            -5.216,-5.213,-5.210,-5.208,-5.206,-5.204,-5.200,-5.196,
            -5.191,-5.187,-5.183,-5.181,-5.178,-5.176,-5.174)
  beta <- c(0.900,0.909,0.927,0.937,0.944,0.948,0.956,
            0.967,0.980,0.995,1.009,1.040,1.070,1.085,1.100,1.129,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151)
  ccr <- c(1.07312,1.0785,1.07254,1.05736,1.03568,
           1.00578,0.98413,0.98059,0.98634,0.99121,1.00033,1.03440,
           1.08388,1.10602,1.12674,1.16455,1.19837,1.27000,1.32852,
           1.37801,1.42087,1.45868,1.49250,1.55102,1.60051,1.64337,
           1.68118,1.71500,1.78663,1.84515,1.93750,2.00913,2.06765,
           2.12357,2.13734,2.13734,2.13734)
  dcr <- c(0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,
           0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.19,0.178,0.162,
           0.148,0.136,0.125,0.101,0.083,0.053,0.02967,0.01078,0.0,0.0,
           0.0,0.0)
  FN_Cr <- c(0.31282,0.31568,0.31854,0.32022,0.32141,
             0.32185,0.32471,0.32860,0.32984,0.33637,0.34096,0.34562,
             0.34596,0.34500,0.34347,0.33911,0.33355,0.31694,0.29912,
             0.28167,0.26521,0.24997,0.23598,0.21157,0.19133,0.17456,
             0.16062,0.14900,0.12761,0.11382,0.09918,0.09329,0.09106,
             0.09003,0.08897,0.08722,0.08444)
  FRV_um <- c(-0.20236,-0.21425,-0.22125,-0.22306,
              -0.22334,-0.22297,-0.22229,-0.22145,-0.22053,-0.21956,-0.21858,
              -0.21661,-0.21469,-0.21374,-0.21282,-0.21102,-0.20929,-0.20526,
              -0.20159,-0.19822,-0.19511,-0.19221,-0.18950,-0.18454,-0.18008,
              -0.17602,-0.17229,-0.16883,-0.16113,-0.15447,-0.14326,-0.13399,
              -0.12604,-0.11905,-0.11280,-0.10713,-0.10194)
  FNS_um <- c(0.25194,0.25837,0.26070,0.26015,0.25891,
              0.25746,0.25596,0.25446,0.25301,0.25160,0.25023,0.24765,
              0.24524,0.24410,0.24299,0.24088,0.23890,0.23437,0.23036,
              0.22674,0.22345,0.22041,0.21760,0.21250,0.20797,0.20388,
              0.20015,0.19671,0.18911,0.18259,0.17171,0.16279,0.15518,
              0.14853,0.14260,0.13724,0.13235)
  bcr <- c(0.00907,0.00907,0.00907,0.00907,0.00907,
           0.00907,0.00907,0.00907,0.00907,0.00907,0.00907,0.00958,
           0.01055,0.01122,0.01170,0.01233,0.01346,0.01617,0.01831,
           0.01980,0.02078,0.02138,0.02168,0.02161,0.02094,0.01987,
           0.01854,0.01704,0.01292,0.00863,0.00042,-0.00687,-0.01315,
           -0.01849,-0.02298,-0.02672,-0.02980)
  gcr <- c(-1.26034,-1.26949,-1.28775,-1.29619,-1.25147,
           -1.14724,-1.09126,-1.04676,-1.01364,-0.99232,-0.98315,-0.97311,
           -0.98324,-0.99261,-1.00423,-1.03299,-1.06493,-1.15136,-1.23723,
           -1.31736,-1.38989,-1.45461,-1.51219,-1.60834,-1.68374,-1.74306,
           -1.78967,-1.82603,-1.88672,-1.91705,-1.93178,-1.91787,-1.89632,
           -1.87695,-1.86168,-1.85415,-1.85292)
  gum <- c(-1.09985,-1.10720,-1.11766,-1.11392,-1.07971,
           -0.98606,-0.93901,-0.90036,-0.87064,-0.85080,-0.84175,-0.82925,
           -0.83484,-0.84153,-0.85033,-0.87328,-0.89945,-0.97276,-1.04809,
           -1.12034,-1.18735,-1.24859,-1.30439,-1.40097,-1.48074,-1.54711,
           -1.60263,-1.64914,-1.73810,-1.79802,-1.86912,-1.90393,-1.92241,
           -1.93474,-1.94360,-1.95343,-1.96352)
  gcrN <- c(-0.49919,-0.48684,-0.44645,-0.42175,
            -0.37622,-0.43576,-0.46789,-0.50260,-0.53802,-0.57301,-0.60839,
            -0.67012,-0.72524,-0.74977,-0.77240,-0.81285,-0.84625,-0.90575,
            -0.93933,-0.95403,-0.95484,-0.94513,-0.92771,-0.87661,-0.81164,
            -0.73894,-0.66214,-0.58347,-0.40237,-0.24053,-0.03590,0.08743,
            0.17408,0.24501,0.31173,0.37949,0.45204)
  gcrL <- c(1.26564,1.24149,1.19894,1.18654,1.14213,
            1.14137,1.15749,1.18679,1.22338,1.26404,1.30526,1.39284,1.47770,
            1.51881,1.55884,1.63476,1.70605,1.86247,1.99162,2.09823,2.18630,
            2.25934,2.31973,2.41072,2.47179,2.51100,2.53418,2.54538,2.54707,
            2.52581,2.49751,2.47889,2.46474,2.45267,2.44271,2.43391,2.42673)
  ecr <- c(-0.00794,-0.00772,-0.00756,-0.00788,-0.00863,
           -0.00988,-0.01075,-0.01130,-0.01164,-0.01182,-0.01184,
           -0.01174,-0.01142,-0.01123,-0.01101,-0.01054,-0.01007,
           -0.00890,-0.00784,-0.00691,-0.00610,-0.00541,-0.00482,
           -0.00387,-0.00315,-0.00261,-0.00220,-0.00189,-0.00139,
           -0.00116,-0.00109,-0.00125,-0.00144,-0.00159,-0.00171,
           -0.00176,-0.00177)
  eum <- c(-0.010830,-0.010580,-0.010505,-0.011005,
           -0.011488,-0.012292,-0.012739,-0.013080,-0.013313,-0.013455,
           -0.013479,-0.013488,-0.013335,-0.013231,-0.013115,-0.012840,
           -0.012558,-0.011817,-0.011100,-0.010436,-0.009831,-0.009285,
           -0.008788,-0.007922,-0.007193,-0.006568,-0.006025,-0.005551,
           -0.004580,-0.003842,-0.002838,-0.002226,-0.001844,-0.001604,
           -0.001488,-0.001452,-0.001499)
  ecrv <- c(-0.00628,-0.00629,-0.00634,-0.00647,
            -0.00681,-0.00710,-0.00724,-0.00734,-0.00741,-0.00746,
            -0.00749,-0.00751,-0.00748,-0.00746,-0.00743,-0.00735,
            -0.00725,-0.00696,-0.00661,-0.00624,-0.00586,-0.00548,
            -0.00511,-0.00439,-0.00374,-0.00315,-0.00261,-0.00214,
            -0.00119,-0.00054,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
  gamma1 <- c(-9.08724,-9.05206,-8.85775,-8.62794,
              -8.39822,-8.15307,-8.00327,-8.04783,-8.19691,-8.35079,-8.5243,
              -9.05478,-9.67640,-9.96595,-10.24160,-10.75200,-11.22492,
              -12.26275,-13.14181,-13.90254,-14.57355,-15.17384,-15.71567,
              -16.66219,-17.46728,-18.16453,-18.77798,-19.32562,-20.46705,
              -21.38822,-22.81170,-23.90512,-24.79448,-25.58872,-26.04635,
              -26.37008,-26.67797)
  SCTerm <- c(0.28877,0.29991,0.29778,0.22535,0.15865,
              0.08257,0.05542,0.04177,0.06040,0.09304,0.15210,0.26161,0.35467,
              0.39410,0.42523,0.47642,0.51169,0.55177,0.55325,0.53555,0.50855,
              0.47913,0.45094,0.39849,0.35254,0.31707,0.28896,0.26685,0.22889,
              0.20778,0.18675,0.17749,0.16807,0.15988,0.15433,0.15116,0.15275,
              0.12210,0.10321,0.11109,0.08972,0.04037,  # start the second column
              -0.05934,-0.10956,-0.13262,-0.13537,-0.12513,-0.07090,0.00261,
              0.11202,0.15670,0.19678,0.27414,0.33785,0.47868,0.57391,
              0.63824,0.67945,0.70517,0.72065,0.72820,0.71582,0.69660,
              0.67366,0.64985,0.59394,0.55120,0.50247,0.48879,0.48245,
              0.47612,0.46978,0.46344,0.45711, # start the third column
              0.20813,0.21925,0.21713,0.15948,0.07029,
              -0.03529,-0.09017,-0.09583,-0.07321,-0.03629,0.01110,0.10298,
              0.21434,0.25574,0.29811,0.36360,0.42568,0.55643,0.65844,
              0.73536,0.79334,0.83853,0.87505,0.92621,0.95627,0.97619,
              0.98819,0.99485,0.99639,0.98814,0.96004,0.92903,0.89351,
              0.85603,0.81676,0.77379,0.72802)
  SCTerm <- matrix(SCTerm, nrow = Numperiod)
  Sigma <- c(0.556,0.556,0.555,0.553,0.558,0.564,0.577,
             0.599,0.616,0.629,0.641,0.657,0.663,0.666,0.671,0.680,0.692,
             0.694,0.688,0.675,0.667,0.665,0.664,0.669,0.670,0.674,0.673,
             0.670,0.660,0.656,0.631,0.605,0.591,0.578,0.556,0.542,0.538)
  tauEvent <- c(0.391,0.390,0.396,0.408,0.438,0.460,
                0.481,0.488,0.481,0.477,0.460,0.437,0.414,0.402,0.384,0.380,
                0.359,0.340,0.344,0.353,0.363,0.359,0.361,0.362,0.370,0.378,
                0.386,0.390,0.399,0.397,0.381,0.376,0.363,0.363,0.376,0.377,
                0.395)
  Imav <- c(3.73,3.07,2.76,3.02)
  phaiSC <- c(3.50,3.00,2.50,3.00)
  gammaSC <- c(0.8,1.00,0.90,0.60)
  LnSCIAm <- c(0.323,0.205,0.083,0.041,0.034,0.046,0.069,
               0.098,0.132,0.169,0.208,0.288,0.370,0.412,0.453,0.535,0.606,
               0.670,0.710,0.71894,0.70561,0.69289,0.68075,0.65804,0.63715,
               0.61781,0.59979,0.58290,0.54475,0.51118,0.45382,0.40562,
               0.36383,0.32682,0.29350,0.26316,0.23525)
  LnAmax1D <- c(0.65022,0.65181,0.65362,0.65467,
                0.65285,0.67264,0.69966,0.71713,0.71603,0.72561,0.74200,
                0.76236,0.75215,0.73819,0.71911,0.65408,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395, # start the second column
                0.70973,0.70679,0.69465,0.68755,
                0.69892,0.70137,0.72445,0.74343,0.78598,0.79721,0.81668,
                0.84523,0.78296,0.79480,0.80861,0.84331,0.87770,0.93767,
                0.95000,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                1.0,1.0,1.0,1.0,1.0,1.0, # start the third column
                0.64434,0.64624,0.63865,0.63421,
                0.60604,0.61716,0.63797,0.65437,0.68019,0.70889,0.71881,
                0.72581,0.74525,0.76103,0.76813,0.75690,0.71785,0.65470,
                0.69619,0.77907,0.82776,0.87645,0.92514,0.97383,1.02252,
                1.07122,1.11991,1.16860,1.21729,1.26598,1.31467,1.36336,
                1.41205,1.46075,1.50944,1.55813,1.60682, # start the forth column
                0.40428,0.40428,0.38789,0.37830,
                0.31737,0.30934,0.32530,0.35412,0.39282,0.42184,0.43736,
                0.47208,0.51278,0.53432,0.55022,0.57279,0.59674,0.61136,
                0.62638,0.63012,0.64773,0.64152,0.65582,0.68668,0.70560,
                0.71429,0.70388,0.67813,0.61119,0.54736,0.45944,0.40846,
                0.36421,0.32984,0.30912,0.29251,0.54736)
  LnAmax1D <- matrix(LnAmax1D, nrow = Numperiod)
  Src1D <- c(8.429,8.090,6.992,6.350,4.883,5.043,
             6.271,7.667,9.034,11.251,14.817,14.817,14.817,14.817,
             14.817,14.817,14.817,14.817,14.817,14.817,14.817,14.817,
             14.817,14.817,14.817,14.817,14.817,14.817,14.817,14.817,
             14.817,14.817,14.817,14.817,14.817,14.817,14.817, # start the second column, these large numbers of 14.817 are not used
             1.91368,1.88256,1.77861,1.71781,
             2.05234,2.38713,2.83399,3.29447,3.99091,4.46576,5.04561,
             5.89960,5.05353,5.20490,5.38694,5.87165,6.57391,8.50000,
             10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,
             10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,
             10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,
             10.67030, # start the third column, these large numbers of 10.67030 are not used
             1.11714,1.11444,1.12437,1.13017,
             1.15080,1.23971,1.34819,1.45181,1.58315,1.73292,1.84134,
             2.03029,2.28133,2.44413,2.58017,2.74161,2.82587,2.71893,
             2.41759,2.30375,2.23625,2.21678,2.24338,2.80535,6.65839,
             30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,
             30.0, # start the forth column, these large numbers of 30.0 are not used
             0.83644,0.83644,0.83000,0.82624,
             0.76758,0.78632,0.83775,0.92616,1.02228,1.11802,1.16578,
             1.28551,1.39808,1.44327,1.47177,1.54694,1.64401,1.79013,
             1.82345,1.79037,1.76844,1.67539,1.62539,1.52453,1.39724,
             1.32029,1.26637,1.22680,1.22065,1.31805,2.12485,14.38181,
             14.38181,14.38181,14.38181,14.38181,14.38181) # these large numbers of 14.38181 are not used
  Src1D <- matrix(Src1D, nrow = Numperiod)
  fsrSC <- c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
             1.0,1.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.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.0, # start the second column, When fsr is zero only linear site terms will be used
             1.0,1.2,1.3,1.253,1.064,1.120,
             1.207,1.238,1.360,1.355,1.341,1.195,0.835,0.781,0.738,0.684,
             0.654,0.683,0.691,0.800,0.876,0.966,1.034,1.206,1.314,1.357,
             1.357,1.357,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, # start the third column
             1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
             1.0,1.080,1.093,0.948,0.908,0.862,0.745,0.623,0.436,0.400,
             0.433,0.460,0.496,0.529,0.578,0.578,0.578,0.578,0.000,0.0,
             0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, # start the forth column
             1.0,1.0,0.949,0.550,0.477,0.492,
             0.531,0.613,0.693,0.780,0.816,0.998,0.954,0.942,0.927,0.927,
             0.949,0.970,0.961,0.948,0.965,0.958,0.984,1.055,1.114,1.175,
             1.216,1.230,1.192,0.942,0.0,0.0,0.0,0.0,0.0,0.0,0.0) # When fsr is zero only linear site terms will be used
  fsrSC <- matrix(fsrSC, nrow = Numperiod)
  fsrUM <- c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
             1.0,1.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.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.0, # start the second column
             1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
             1.0,1.0,1.0,0.767,0.708,0.657,0.568,0.509,0.433,0.337,0.337,
             0.337,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,
             0.0,0.0,0.0, # start the third column
             1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.83,
             0.706,0.8,0.759,0.715,0.686,0.644,0.549,0.434,0.292,0.275,
             0.295,0.293,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,0.0,0.0,0.0, # start the forth column
             1.0,1.0,1.15,0.8,0.613,0.542,0.534,
             0.583,0.643,0.674,0.694,0.763,0.684,0.645,0.61,0.532,0.468,
             0.291,0.275,0.295,0.293,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,0.0,0.0,0.0) # When fsr is zero only linear site terms will be used
  fsrUM <- matrix(fsrUM, nrow = Numperiod)
  SigmaSC <- c(0.4153,0.4142,0.4123,0.4089,0.4072,
               0.4059,0.4066,0.4171,0.4255,0.4315,0.4371,0.4515,0.4559,
               0.4594,0.4647,0.4721,0.4722,0.4780,0.4906,0.4918,0.4951,
               0.4927,0.4890,0.4792,0.4665,0.4600,0.4541,0.4475,0.4404,
               0.4284,0.4148,0.4003,0.4033,0.4026,0.3804,0.3736,0.3831, # start the second column
               0.4363,0.4356,0.4344,0.4307,0.4275,
               0.4294,0.4304,0.4329,0.4361,0.4417,0.4478,0.4476,0.4553,
               0.4547,0.4607,0.4697,0.4800,0.4882,0.4991,0.5043,0.5074,
               0.5074,0.5149,0.5177,0.5044,0.4954,0.4834,0.4720,0.4595,
               0.4499,0.4374,0.4343,0.4328,0.4214,0.4264,0.4290,0.4378, # start the third column
               0.4235,0.4230,0.4226,0.4207,0.4208,
               0.4168,0.4156,0.4215,0.4220,0.4161,0.4241,0.4332,0.4347,
               0.4380,0.4419,0.4514,0.4493,0.4626,0.4833,0.4783,0.4720,
               0.4944,0.4949,0.4848,0.4816,0.4751,0.4554,0.4415,0.4380,
               0.4219,0.4172,0.3955,0.3995,0.4131,0.4060,0.3921,0.3779, # start the forth column
               0.4426,0.4419,0.4418,0.4415,0.4435,
               0.4446,0.4415,0.4407,0.4351,0.4374,0.4402,0.4594,0.4697,
               0.4707,0.4744,0.4784,0.4825,0.4843,0.4930,0.4914,0.4783,
               0.4787,0.4764,0.4854,0.4707,0.4837,0.4729,0.4709,0.4627,
               0.4564,0.4323,0.4128,0.4170,0.4265,0.4259,0.4256,0.4195)
  SigmaSC <- matrix(SigmaSC, nrow = Numperiod)
  tauSC <- c(0.3581,0.3581,0.3581,0.3625,0.3832,
             0.4052,0.4341,0.4517,0.4745,0.4884,0.4932,0.4927,0.4906,
             0.4815,0.4725,0.4640,0.4708,0.4650,0.4595,0.4432,0.4296,
             0.4310,0.4323,0.4402,0.4390,0.4491,0.4579,0.4663,0.4644,
             0.4546,0.4190,0.4084,0.4052,0.3861,0.3756,0.3595,0.3594, # start the second column
             0.3620,0.3627,0.3646,0.3580,0.3666,
             0.3760,0.3911,0.4057,0.4247,0.4445,0.4485,0.4641,0.4758,
             0.4817,0.4820,0.4830,0.4932,0.5436,0.5246,0.4872,0.4722,
             0.4651,0.4493,0.4472,0.4515,0.4504,0.4559,0.4593,0.4606,
             0.4688,0.4530,0.4310,0.4152,0.3953,0.3681,0.3589,0.3298, # start the third column
             0.3486,0.3482,0.3483,0.3544,0.3544,
             0.3508,0.3453,0.3605,0.3821,0.3958,0.4229,0.4768,0.4592,
             0.4519,0.4617,0.4660,0.4452,0.4192,0.4356,0.4689,0.4821,
             0.4638,0.4327,0.4577,0.4712,0.4622,0.4781,0.4743,0.4495,
             0.4633,0.4718,0.4530,0.4030,0.3803,0.3228,0.3069,0.3076, # start the forth column
             0.2778,0.2804,0.2823,0.2857,0.2881,
             0.2945,0.3082,0.3378,0.3557,0.3784,0.4182,0.4216,0.3869,
             0.3867,0.3788,0.3903,0.4010,0.3942,0.3725,0.3421,0.3466,
             0.3594,0.3998,0.4454,0.4982,0.5183,0.5358,0.5240,0.4932,
             0.4942,0.5195,0.5059,0.4886,0.4928,0.4510,0.4386,0.3774)
  tauSC <- matrix(tauSC, nrow = Numperiod)

  # calculation
  Alfa1D <- 2.0
  Beta1D <- 0.6
  xcr0 <- 2.0
  magc <- 7.1
  Imin <- 1.0
  Imax <- 12.0
  DepthC <- 25.0
  if((fDepth > 25 && eqkType == 1) || eqkType > 2 || eqkType == 0){
    print('The depth does not match eqk type or eqk type is not correct!')
    stop()
  }
  for(i in seq(1, Numperiod)){
    # Bilinear magnitude scalling term and magnitude dependent dist cons
    ############ source effect ##############
    if(mag <= magc){
      r0 <- xcr0 + dist + exp(alfa[i] + beta[i]*mag)
      SpeCalLog[i]=ccr[i]*mag
    }else{
      r0 <- xcr0 + dist + exp(alfa[i]+beta[i]*magc)
      SpeCalLog[i] <- ccr[i]*magc+dcr[i]*(mag-magc)
    }
    ## also add on some path effects "gcr/gum", "ecr/eum"
    if(eqkType == 1){  # shallow crustal
      SpeCalLog[i]=SpeCalLog[i]+bcr[i]*fDepth+gcr[i]*log(r0)+ecr[i]*dist # Reverse and Strike Slip
      if(FaultMech == 3) SpeCalLog[i] <- SpeCalLog[i]+FN_Cr[i] # Normal Fault
    }else{ # upper mantle
      if(FaultMech == 1) SpeCalLog[i] <- SpeCalLog[i]+FRV_um[i] # Reverse Fault
      if (FaultMech >= 2) SpeCalLog[i] <- SpeCalLog[i]+FNS_um[i] # Normal and Strike Slip
      SpeCalLog[i] <- SpeCalLog[i] + gum[i]*log(r0) + eum[i]*dist
    }
    ######################################
    ########### add other path effects #############
    # "grcL", "gamma1-gamma_cr", "gcrN-gN", "ecrv"
    SpeCalLog[i] <- SpeCalLog[i]+gcrL[i]*log(dist+200.0)+gamma1[i]
    if(dist <= 30){
      SpeCalLog[i] <- SpeCalLog[i]+gcrN[i]*log(xcr0+dist+exp(alfa[i]+6.5*beta[i]))
    }else{
      SpeCalLog[i]=SpeCalLog[i]+gcrN[i]*log(xcr0+30.0+exp(alfa[i]+6.5*beta[i]))
    }
    if(DistV > 0) SpeCalLog[i]=SpeCalLog[i]+ecrv[i]*DistV
    ######################################
    ########### add Site class terms / Site amplification ###############
    LnANmax=LnSCIAm[i] # linear site amplification
    if(SiteClass >= 2) LnANmax <- LnANmax+SCTerm[i,SiteClass-1]
    RockSpeCalLog[i] <- SpeCalLog[i]-LnSCIAm[i] # reference rock spectral response
    nSReff <- exp(RockSpeCalLog[i])/phaiSC[SiteClass]
    nSReffc <- Src1D[i,SiteClass]/phaiSC[SiteClass]
    if(Imav[SiteClass] > Imax){
      nSReff <- nSReff*(1.0+gammaSC[SiteClass]*(Imax-Imin))
      nSReffc <- nSReffc*(1.0+gammaSC[SiteClass]*(Imax-Imin))
    }else if(Imav[SiteClass] > Imin){
      nSReff <- nSReff*(1.0+gammaSC[SiteClass]*(Imav[SiteClass]-Imin))
      nSReffc <- nSReffc*(1.0+gammaSC[SiteClass]*(Imav[SiteClass]-Imin))
    }
    if(SiteClass == 1){
      LnSF <- LnSCIAm[i]-LnAmax1D[i,SiteClass]
    }else{
      LnSF <- LnSCIAm[i]+SCTerm[i,SiteClass-1]-LnAmax1D[i,SiteClass]
    }
    if(exp(LnANmax) < 1.25){
      ca <- LnAmax1D[i,SiteClass]/(log(Beta1D)-log(nSReffc^Alfa1D+Beta1D))
      cb <- -ca*log(nSReffc^Alfa1D+Beta1D)
      SNC=exp((ca*(Alfa1D-1.0)*log(Beta1D)*log(10.0*Beta1D)-
                 log(10.0)*(cb+LnSF))/(ca*(Alfa1D*log(10.0*Beta1D)-log(Beta1D))))
    }else{
      SNC=(exp((LnANmax*log(nSReffc^Alfa1D+Beta1D)-
                  LnSF*log(Beta1D))/LnAmax1D[i,SiteClass])-Beta1D)^(1.0/Alfa1D)
    }
    if(eqkType == 1) SMR <- nSReff*SNC/nSReffc*fsrSC[i,SiteClass]
    if(eqkType == 2) SMR <- nSReff*SNC/nSReffc*fsrUM[i,SiteClass]
    if(SMR != 0){
      SpeCalLog[i]=RockSpeCalLog[i]+LnANmax-LnAmax1D[i,SiteClass]*
        (log(SMR^Alfa1D+Beta1D)-log(Beta1D))/(log(nSReffc^Alfa1D+Beta1D)-log(Beta1D))
      SiteAmpCalLog[i] <- LnANmax-LnAmax1D[i,SiteClass]*
        (log(SMR^Alfa1D+Beta1D)-log(Beta1D))/(log(nSReffc^Alfa1D+Beta1D)-log(Beta1D))
    }else{
      SpeCalLog[i] <- RockSpeCalLog[i]+LnANmax
      SiteAmpCalLog[i] <- LnANmax
    }
    SigmaT[i] <- sqrt(Sigma[i]^2+tauEvent[i]^2)
    sigmaSS[i] <- sqrt(SigmaSC[i,SiteClass]^2+tauEvent[i]^2)
  }


  # Compute Sa and sigma with pre-defined period
  if (length(T) == 1 & T[1] == 1000) {

    T <- period

  } else {

    SpeCalLog <- approx(x = period, y = SpeCalLog, xout = T, rule = 2)$y

    RockSpeCalLog <- approx(x = period, y = RockSpeCalLog, xout = T, rule = 2)$y

    SiteAmpCalLog <- approx(x = period, y = SiteAmpCalLog, xout = T, rule = 2)$y

    SigmaT <- approx(x = period, y = SigmaT, xout = T, rule = 2)$y

    sigmaSS <- approx(x = period, y = sigmaSS, xout = T, rule = 2)$y
  }

  # output
  res <- list()

  res$Sa <- exp(SpeCalLog)

  res$RockSa <- exp(RockSpeCalLog)

  res$SiteAmp <- SiteAmpCalLog

  res$sigma <- SigmaT

  res$phiSS <- sigmaSS

  res$period <- T

  return(res)
}


#' The GMPE of Zea 2016 for interface event
#'
#' This function calculates the ground motion median values and standard deviations
#' @param mag Moment magnitude, a numeric value
#' @param T Period (sec) (range should be between 0.005 s and 5 s).
#' Use 1000 (by default) for output the array of Sa with original Zhao et al per-defined periods
#' @param fDepth Depth(km) to the top of ruptured plane (km)
#' @param dist The closest distance (km) to rupture plane when is available,
#' otherwise hypocentral distance.
#' @param DistV The distance of volcanic path (km)
#' @param SiteClass The site class: 1 for SC I, 2 for SC II, 3 for SC III 4 for SC IV sites. The relationship
#' between site class and Vs30 is given by Zhao et al., 2016 (https://doi.org/10.1785/0120150056).
#' SC I: Vs30 > 600 m/s; SC II: 300 < Vs30 =< 600 m/s; SC III: 200 < Vs30 =< 300 m/s; SC IV: Vs30 <= 200 m/s
#' @return A list of six elements is returned: Sa - median spectral acceleration prediction (in g);
#' RockSa - median spectral acceleration prediction on Rock reference site (in g);
#' SiteAmp - Site amplication response spectral (subtraction in log(g)); sigma - totla standard deviation;
#' phiSS - single station standard deviation; period - the corresponding oscillator periods
#' @examples zea_2016_Sub_Interf(mag = 7.5, T = 1000, fDepth = 0.0, dist = 40.0,
#' DistV = 0.0, SiteClass = 4)
#'
#' zea_2016_Sub_Interf(mag = 7.5, T = c(0.025, 0.5), fDepth = 0.0, dist = 40.0,
#' DistV = 0.0, SiteClass = 4)
#' @references Zhao, et al. (2016). Ground‐Motion Prediction Equations for Subduction Interface Earthquakes
#' in Japan Using Site Class and Simple Geometric Attenuation Functions
#' Bulletin of the Seismological Society of America. 106(4): 1518-1534
#' @export
#' @importFrom stats approx
zea_2016_Sub_Interf <- function(mag, T = 1000, fDepth, dist, DistV, SiteClass){

  period <- c(0.005,0.01,0.02,0.03,0.04,0.05,0.06,
              0.07,0.08,0.09,0.1,0.12,0.14,0.15,0.16,0.18,0.2,0.25,0.3,
              0.35,0.4,0.45,0.5,0.6,0.7,0.8,0.9,1.0,1.25,1.5,2.0,2.5,
              3.0,3.5,4.0,4.5,5.0)

  g <- 980.7

  Numperiod <- length(period)

  SpeCalLog <- rep(0, Numperiod)

  RockSpeCalLog <- rep(0, Numperiod)

  SiteAmpCalLog <- rep(0, Numperiod)

  SigmaT <- rep(0, Numperiod)

  sigmaSS <- rep(0, Numperiod)

  # model coefficients
  alfa <- c(-5.30119,-5.28844,-5.27568,-5.26822,
            -5.26293,-5.25882,-5.25547,-5.25263,-5.25017,-5.24801,
            -5.24607,-5.24271,-5.23988,-5.23861,-5.23742,-5.23525,
            -5.23331,-5.22921,-5.22585,-5.22302,-5.22056,-5.21839,
            -5.21645,-5.21310,-5.21026,-5.20781,-5.20564,-5.20370,
            -5.19959,-5.19624,-5.19095,-5.18684,-5.18349,-5.18065,
            -5.17819,-5.17602,-5.17409)
  beta <- c(1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151)
  cInt <- c(1.09973,1.09848,1.09227,1.10690,1.11578,
            1.10234,1.08611,1.07291,1.06383,1.05863,1.05673,1.06051,
            1.07135,1.07856,1.08664,1.10470,1.12443,1.17689,1.22970,
            1.28058,1.32873,1.37394,1.41630,1.49310,1.56069,1.62055,
            1.67390,1.72171,1.82188,1.90081,2.01482,2.08892,2.13574,
            2.16254,2.17393,2.17301,2.16199)
  cIntS <- c(1.31479,1.31739,1.31919,1.34096,1.38051,
             1.43246,1.46238,1.47120,1.46432,1.44695,1.42323,1.36833,
             1.31558,1.29277,1.27319,1.24828,1.23715,1.22387,1.22846,
             1.24219,1.26077,1.28190,1.30432,1.35019,1.39523,1.43824,
             1.47881,1.51685,1.60148,1.67280,1.78372,1.86241,1.91711,
             1.95322,1.97450,1.98361,1.98255)
  dInt1 <- c(0.553,0.553,0.553,0.553,0.553,0.553,0.553,
             0.553,0.553,0.553,0.553,0.553,0.553,0.553,0.553,0.553,
             0.553,0.553,0.553,0.553,0.553,0.553,0.553,0.553,0.560,
             0.580,0.602,0.622,0.667,0.705,0.768,0.820,0.863,0.902,
             0.935,0.966,0.994)
  gammaIntS <- c(-3.89528,-3.89528,-3.89528,-3.89528,
                 -3.89528,-3.89528,-3.89528,-3.89528,-3.89461,-3.90179,
                 -3.90765,-3.91644,-3.92273,-3.92526,-3.92753,-3.93130,
                 -3.93446,-3.94068,-3.94547,-3.94943,-3.95273,-3.95558,
                 -3.95795,-3.96181,-3.96483,-3.96729,-3.96962,-3.97202,
                 -3.97949,-3.99048,-4.02652,-4.08299,-4.15939,-4.25418,
                 -4.36584,-4.49267,-4.63313)
  bInt <- c(0.01999,0.01999,0.01999,0.02066,0.02308,
            0.02709,0.02972,0.03207,0.03196,0.02972,0.02789,0.02470,
            0.02117,0.01951,0.01793,0.01505,0.01255,0.00769,0.00438,
            0.00215,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,0.0,0.0,0.0,0.0)
  gint <- c(-2.05587,-2.06565,-2.10231,-2.19232,
            -2.24640,-2.29341,-2.31172,-2.31101,-2.28778,-2.24680,-2.20410,
            -2.12011,-2.04337,-2.01088,-1.98301,-1.94607,-1.92702,-1.89876,
            -1.89141,-1.89299,-1.89531,-1.90578,-1.91467,-1.92743,-1.93452,
            -1.93739,-1.93725,-1.93505,-1.92471,-1.91189,-1.88859,-1.87252,
            -1.86345,-1.85969,-1.85953,-1.86151,-1.86451)
  gintLD <- c(0.54541,0.54975,0.56172,0.57887,
              0.49333,0.49100,0.50848,0.52750,0.54595,0.56307,0.57615,
              0.59258,0.60983,0.61961,0.63085,0.66196,0.69978,0.78447,
              0.85939,0.92338,0.98011,1.02220,1.05874,1.11803,1.16302,
              1.19735,1.22362,1.24367,1.27245,1.28541,1.28830,1.27727,
              1.26049,1.24107,1.22026,1.19859,1.17629)
  gintLS <- c(1.13364,1.13365,1.13365,1.13364,0.98811,
              0.90444,0.88766,0.90491,0.94206,0.98652,1.03553,1.13529,
              1.23424,1.28134,1.32659,1.41134,1.48854,1.65209,1.78127,
              1.88438,1.96763,2.03554,2.09142,2.17640,2.23600,2.27831,
              2.30847,2.32985,2.35853,2.36653,2.35536,2.33106,2.30409,
              2.27794,2.25371,2.23161,2.21149)
  eintv <- c(-0.01123,-0.01125,-0.01127,-0.01158,
             -0.01203,-0.01256,-0.01312,-0.01359,-0.01382,-0.01393,-0.01395,
             -0.01381,-0.01351,-0.01333,-0.01312,-0.01269,-0.01223,-0.01108,
             -0.00998,-0.00898,-0.00808,-0.00727,-0.00656,-0.00534,-0.00437,
             -0.00359,-0.00296,-0.00244,-0.00153,-0.00097,-0.00043,-0.00023,
             -0.00016,0.0,0.0,0.0,0.0)
  eintS <- c(-0.00628,-0.00625,-0.00616,-0.00572,
             -0.00532,-0.00503,-0.00528,-0.00569,-0.00619,-0.00673,-0.00718,
             -0.00793,-0.00853,-0.00879,-0.00902,-0.00927,-0.00942,-0.00959,
             -0.00952,-0.00933,-0.00911,-0.00888,-0.00866,-0.00824,-0.00787,
             -0.00755,-0.00726,-0.00700,-0.00644,-0.00597,-0.00518,-0.00451,
             -0.00393,-0.00344,-0.00302,-0.00267,-0.00240)
  Imav <- c(3.73,3.07,2.76,3.02)
  phaiSC <- c(3.50,3.00,2.50,3.00)
  gammaSC <- c(0.8,1.00,0.90,0.60)
  gammaInt <- c(-4.49858,-4.45894,-4.25807,-3.91800,
                -3.11423,-2.76035,-2.64088,-2.65622,-2.75266,-2.89920,-3.07698,
                -3.48283,-3.91612,-4.13481,-4.35243,-4.78030,-5.19439,-6.15803,
                -7.02003,-7.79153,-8.49551,-9.11351,-9.68515,-10.68946,
                -11.54596,-12.28722,-12.93630,-13.51000,-14.69034,-15.60298,
                -16.90010,-17.73658,-18.27135,-18.59264,-18.75474,-18.79351,
                -18.73393)
  SCTerm <- c(0.31288,0.30850,0.29297,0.22869,0.16316,
              0.12129,0.12345,0.13970,0.16385,0.20501,0.24449,0.32284,0.40120,
              0.43622,0.46736,0.51201,0.53926,0.58602,0.60465,0.60640,0.60277,
              0.58043,0.55692,0.50969,0.46501,0.42440,0.38841,0.35703,0.29673,
              0.25785,0.22262,0.21840,0.21595,0.21595,0.21595,0.21595,0.21595, # start the second column
              -0.00431,-0.01115,-0.02166,-0.11291,
              -0.18870,-0.22831,-0.21920,-0.19009,-0.15499,-0.11250,-0.07508,
              0.01500,0.09699,0.14588,0.18793,0.25154,0.30298,0.42692,0.51620,
              0.56951,0.62366,0.65814,0.68671,0.71220,0.71239,0.69938,0.67996,
              0.65828,0.61928,0.58288,0.52620,0.48720,0.45698,0.42809,0.39524,
              0.35488,0.30468, # start the third column
              0.22838,0.22305,0.20892,0.13314,0.06959,
              0.02845,0.00070,-0.00949,-0.00485,0.03047,0.06080,0.14231,
              0.22696,0.25758,0.30748,0.35974,0.40313,0.50765,0.57779,0.63821,
              0.70323,0.75082,0.79378,0.84953,0.87979,0.89539,0.90256,0.90528,
              0.91793,0.92130,0.91709,0.90551,0.88668,0.85880,0.82025,0.76990,
              0.70705, # start the forth column
              0.31288,0.30850,0.29297,0.22869,0.16316,
              0.12129,0.12345,0.13970,0.16385,0.20501,0.24449,0.32284,0.40120,
              0.43622,0.46736,0.51201,0.53926,0.58602,0.60465,0.60640,0.60277,
              0.58043,0.55692,0.50969,0.46501,0.42440,0.38841,0.35703,0.29673,
              0.25785,0.22262,0.21840,0.21595,0.21595,0.21595,0.21595,0.21595, # start the fifth column
              -0.00431,-0.01115,-0.02166,-0.08291,
              -0.18870,-0.22831,-0.21920,-0.19009,-0.15499,-0.11250,-0.07508,
              0.01500,0.09699,0.13588,0.17293,0.24154,0.30298,0.42692,0.51620,
              0.57951,0.63366,0.65814,0.68671,0.71220,0.71239,0.69938,0.67996,
              0.65828,0.61928,0.58288,0.52620,0.48720,0.45698,0.42809,0.39524,
              0.35488,0.30468, # start the sixth column
              0.22838,0.22305,0.20892,0.16314,0.05959,
              0.00845,0.00070,0.00551,0.01415,0.04047,0.07080,0.14031,0.20196,
              0.23258,0.26248,0.31974,0.37313,0.48765,0.57779,0.64821,0.71323,
              0.75082,0.79378,0.84953,0.87979,0.89539,0.90256,0.90528,0.91793,
              0.92130,0.91709,0.90551,0.88668,0.85880,0.82025,0.76990,0.70705)
  SCTerm <- matrix(SCTerm, nrow = Numperiod)
  Sigma <- c(0.553,0.554,0.553,0.555,0.565,0.570,0.583,
             0.602,0.614,0.625,0.637,0.646,0.654,0.659,0.663,0.672,0.678,
             0.659,0.640,0.634,0.627,0.620,0.612,0.613,0.625,0.628,0.628,
             0.633,0.636,0.644,0.635,0.619,0.599,0.581,0.568,0.551,0.563)
  tauEvent <- c(0.377,0.377,0.384,0.397,0.428,0.463,
                0.488,0.501,0.501,0.495,0.478,0.453,0.412,0.404,0.398,0.387,
                0.382,0.365,0.348,0.36,0.354,0.363,0.364,0.379,0.393,0.396,
                0.397,0.403,0.404,0.392,0.382,0.393,0.385,0.376,0.377,0.376,
                0.374)
  LnSCIAm <- c(0.30591,0.22112,0.13918,0.09281,0.06270,
               0.03601,0.03723,0.04916,0.09833,0.17538,0.24472,0.33670,0.42199,
               0.45620,0.48578,0.53379,0.57004,0.62479,0.65053,0.66174,0.66495,
               0.66543,0.66353,0.65646,0.64752,0.63782,0.62753,0.61669,0.58711,
               0.55414,0.48232,0.41085,0.34770,0.29811,0.26524,0.25089,0.23654)
  LnAmax1D <- c(0.65022,0.65181,0.65362,0.65467,
                0.65285,0.67264,0.69966,0.71713,0.71603,0.72561,0.74200,0.76236,
                0.75215,0.73819,0.71911,0.65408,0.58395,0.58395,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,
                0.58395, # start the second column
                0.70973,0.70679,0.69465,0.68755,
                0.69892,0.70137,0.72445,0.74343,0.78598,0.79721,0.81668,0.84523,
                0.78296,0.79480,0.80861,0.84331,0.87770,0.93767,0.95000,1.0,1.0,
                1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, # start the third column
                0.64434,0.64624,0.63865,0.63421,
                0.60604,0.61716,0.63797,0.65437,0.68019,0.70889,0.71881,0.72581,
                0.74525,0.76103,0.76813,0.75690,0.71785,0.65470,0.69619,0.77907,
                0.82776,0.87645,0.92514,0.97383,1.02252,1.07122,1.11991,1.16860,
                1.21729,1.26598,1.31467,1.36336,1.41205,1.46075,1.50944,1.55813,
                1.60682, # start the forth column
                0.40428,0.40428,0.38789,0.37830,
                0.31737,0.30934,0.32530,0.35412,0.39282,0.42184,0.43736,0.47208,
                0.51278,0.53432,0.55022,0.57279,0.59674,0.61136,0.62638,0.63012,
                0.64773,0.64152,0.65582,0.68668,0.70560,0.71429,0.70388,0.67813,
                0.61119,0.54736,0.45944,0.40846,0.36421,0.32984,0.30912,0.29251,
                0.54736)
  LnAmax1D <- matrix(LnAmax1D, nrow = Numperiod)
  Src1D <- c(8.42900,8.09039,6.99211,6.34965,
             4.88339,5.04257,6.27132,7.66694,9.03436,11.25066,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739, # start the second column, nearly all of these large numbers of 14.81739 are not used
             1.91368,1.88256,1.77861,1.71781,
             2.05234,2.38713,2.83399,3.29447,3.99091,4.46576,5.04561,5.89960,
             5.05353,5.20490,5.38694,5.87165,6.57391,8.50000,10.67030,
             10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,
             10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,10.67030,
             10.67030,10.67030,10.67030,10.67030, # start the third column, nearly all of these large numbers of 10.67030, are not used
             1.11714,1.11444,1.12437,1.13017,
             1.15080,1.23971,1.34819,1.45181,1.58315,1.73292,1.84134,2.03029,
             2.28133,2.44413,2.58017,2.74161,2.82587,2.71893,2.41759,2.30375,
             2.23625,2.21678,2.24338,2.80535,6.65839,30.0,30.0,30.0,30.0,
             30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0, # start the forth column, These large numbers of 30.0, are not used
             0.83644,0.83644,0.83000,0.82624,
             0.76758,0.78632,0.83775,0.92616,1.02228,1.11802,1.16578,1.28551,
             1.39808,1.44327,1.47177,1.54694,1.64401,1.79013,1.82345,1.79037,
             1.76844,1.67539,1.62539,1.52453,1.39724,1.32029,1.26637,1.22680,
             1.22065,1.31805,2.12485,14.38181,14.38181,14.38181,14.38181,
             14.38181,14.38181) # nearly all of these large numbers of 14.38181, are not used
  Src1D <- matrix(Src1D, nrow = Numperiod)
  fsr <- c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
           1.0,1.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.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.0, # start the second column, When fsr is zero only linear site terms will be used
           1.0,1.0,1.0,1.0,0.843,0.663,0.841,
           1.029,1.235,1.144,1.092,0.945,0.624,0.577,0.545,0.527,0.546,
           0.596,0.623,0.701,0.708,0.737,0.748,0.728,0.634,0.0,0.0,0.0,
           0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, # start the third column
           1.165,0.944,1.012,1.1,0.959,0.889,
           0.946,1.006,1.065,1.093,1.077,1.036,0.895,0.860,0.822,0.743,
           0.663,0.487,0.447,0.473,0.487,0.511,0.536,0.540,0.477,0.0,
           0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, # start the foth column
           1.0,1.0,1.0,1.0,0.557,0.543,0.574,
           0.648,0.721,0.809,1.015,0.972,0.967,0.963,0.953,0.967,1.005,
           1.045,1.035,1.008,1.007,0.981,0.990,1.016,1.022,1.023,0.997,
           0.948,0.802,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0) # When fsr is zero only linear site terms will be used
  fsr <- matrix(fsr, nrow = Numperiod)
  SigmaSC <- c(0.398,0.397,0.395,0.389,0.387,0.387,
               0.397,0.403,0.413,0.422,0.429,0.440,0.441,0.450,0.452,0.454,
               0.462,0.474,0.472,0.468,0.457,0.450,0.445,0.448,0.440,0.441,
               0.435,0.427,0.413,0.416,0.409,0.398,0.390,0.386,0.377,0.360,
               0.361, # start the second column
               0.417,0.417,0.417,0.418,0.420,0.422,
               0.416,0.412,0.412,0.412,0.418,0.429,0.435,0.440,0.444,0.447,
               0.450,0.470,0.475,0.478,0.484,0.477,0.470,0.460,0.460,0.460,
               0.456,0.448,0.442,0.448,0.437,0.429,0.423,0.410,0.410,0.412,
               0.447, # start the third column
               0.409,0.409,0.408,0.409,0.413,0.409,
               0.401,0.394,0.389,0.391,0.394,0.428,0.433,0.420,0.424,0.446,
               0.441,0.459,0.437,0.444,0.453,0.477,0.472,0.459,0.461,0.457,
               0.449,0.442,0.428,0.418,0.406,0.387,0.369,0.376,0.362,0.368,
               0.381, # start the forth column
               0.415,0.415,0.416,0.417,0.420,0.422,
               0.423,0.420,0.423,0.427,0.426,0.445,0.442,0.440,0.437,0.436,
               0.432,0.432,0.432,0.432,0.416,0.408,0.409,0.407,0.403,0.407,
               0.408,0.408,0.412,0.418,0.413,0.414,0.412,0.402,0.394,0.385,
               0.380)
  SigmaSC <- matrix(SigmaSC, nrow = Numperiod)
  tauSC <- c(0.511,0.517,0.518,0.537,0.572,0.586,
             0.613,0.633,0.643,0.631,0.626,0.614,0.613,0.602,0.599,0.599,
             0.590,0.555,0.532,0.505,0.479,0.461,0.451,0.436,0.430,0.427,
             0.431,0.440,0.440,0.449,0.441,0.425,0.396,0.386,0.373,0.364,
             0.359, # start the second column
             0.449,0.450,0.449,0.449,0.456,0.479,
             0.507,0.532,0.555,0.568,0.566,0.570,0.597,0.599,0.597,0.601,
             0.596,0.601,0.557,0.518,0.488,0.479,0.471,0.467,0.464,0.456,
             0.466,0.466,0.476,0.473,0.462,0.427,0.407,0.401,0.386,0.373,
             0.322, # start the third column
             0.431,0.431,0.431,0.430,0.428,0.429,
             0.449,0.456,0.455,0.475,0.507,0.548,0.507,0.490,0.494,0.500,
             0.491,0.450,0.499,0.532,0.546,0.529,0.495,0.470,0.473,0.457,
             0.439,0.442,0.432,0.461,0.472,0.484,0.447,0.434,0.406,0.383,
             0.316, # start the forth column
             0.422,0.418,0.425,0.422,0.431,0.439,
             0.445,0.473,0.491,0.518,0.559,0.576,0.561,0.555,0.552,0.553,
             0.562,0.507,0.489,0.463,0.466,0.461,0.467,0.434,0.430,0.454,
             0.456,0.462,0.444,0.443,0.439,0.427,0.431,0.413,0.391,0.366,
             0.324)
  tauSC <- matrix(tauSC, nrow = Numperiod)

  Alfa1D=2.0
  Beta1D=0.6
  xcr0=10.0
  magc=7.1
  Imin=1.0
  Imax=12.0
  DepthC=25.0

  # calculation
  for(i in seq(1, Numperiod)){
    # Bilinear magnitude scalling term and magnitude dependent dist cons
    ############ source effect ##############
    if(mag <= magc){
      r0 <- xcr0+dist + exp(alfa[i]+beta[i]*mag)
      if(fDepth <= DepthC){
        SpeCalLog[i]=cIntS[i]*mag
      }else{
        SpeCalLog[i]=cInt[i]*mag
      }
    }else{
      r0 <- xcr0+dist + exp(alfa[i]+beta[i]*magc)
      if(fDepth <= DepthC){
        SpeCalLog[i]=cIntS[i]*magc+dInt1[i]*(mag-magc)
      }else{
        SpeCalLog[i]=cInt[i]*magc+dInt1[i]*(mag-magc)
      }
    }
    SpeCalLog[i]=SpeCalLog[i]+bInt[i]*fDepth
    #################################
    ############ add path effects ############
    SpeCalLog[i]=SpeCalLog[i]+gint[i]*log(r0)+eintv[i]*DistV+gammaInt[i]
    if(fDepth <= DepthC){
      SpeCalLog[i]=SpeCalLog[i]+gammaIntS[i]+gintLS[i]*log(dist+200.0)+eintS[i]*dist
    }else{
      SpeCalLog[i]=SpeCalLog[i]+gintLD[i]*log(dist+200.0)
    }
    ################################
    ############ add site amplification effects ############
    LnANmax=LnSCIAm[i]
    if(SiteClass >= 2){
      if(fDepth < DepthC){
        LnANmax=LnANmax+SCTerm[i,SiteClass-1]
      } else{
        LnANmax=LnANmax+SCTerm[i,SiteClass+2]
      }
    }
    RockSpeCalLog[i]=SpeCalLog[i]-LnSCIAm[i]
    nSReff=exp(RockSpeCalLog[i])/phaiSC[SiteClass]
    nSReffc=Src1D[i,SiteClass]/phaiSC[SiteClass]
    if(Imav[SiteClass] > Imax){
      nSReff=nSReff*(1.0+gammaSC[SiteClass]*(Imax-Imin))
      nSReffc=nSReffc*(1.0+gammaSC[SiteClass]*(Imax-Imin))
    }else if(Imav[SiteClass] > Imin){
      nSReff=nSReff*(1.0+gammaSC[SiteClass]*(Imav[SiteClass]-Imin))
      nSReffc=nSReffc*(1.0+gammaSC[SiteClass]*(Imav[SiteClass]-Imin))
    }
    if(SiteClass == 1){
      LnSF=LnSCIAm[i]-LnAmax1D[i,SiteClass]
    }else{
      if(SiteClass >= 2){
        if(fDepth < DepthC){
          LnSF=LnSCIAm[i]+SCTerm[i,SiteClass-1]-LnAmax1D[i,SiteClass]
        }else{
          LnSF=LnSCIAm[i]+SCTerm[i,SiteClass+2]-LnAmax1D[i,SiteClass]
        }
      }
    }
    if(exp(LnANmax) < 1.25){
      ca=LnAmax1D[i,SiteClass]/(log(Beta1D)-log(nSReffc^Alfa1D+Beta1D))
      cb=-ca*log(nSReffc^Alfa1D+Beta1D)
      SNC=exp((ca*(Alfa1D-1.0)*log(Beta1D)*log(10.0*Beta1D)-
                 log(10.0)*(cb+LnSF))/(ca*(Alfa1D*log(10.0*Beta1D)-log(Beta1D))))
    }else{
      SNC=(exp((LnANmax*log(nSReffc^Alfa1D+Beta1D)-
                  LnSF*log(Beta1D))/LnAmax1D[i,SiteClass])-Beta1D)^(1.0/Alfa1D)
    }
    SMR=nSReff*SNC/nSReffc*fsr[i,SiteClass]
    if(SMR != 0){
      SpeCalLog[i]=RockSpeCalLog[i]+LnANmax-LnAmax1D[i,SiteClass]*(log(SMR^Alfa1D+Beta1D)-log(Beta1D))/(log(nSReffc^Alfa1D+Beta1D)-log(Beta1D))
      SiteAmpCalLog[i] <- LnANmax-LnAmax1D[i,SiteClass]*
        (log(SMR^Alfa1D+Beta1D)-log(Beta1D))/(log(nSReffc^Alfa1D+Beta1D)-log(Beta1D))
    }else{
      SpeCalLog[i]=RockSpeCalLog[i]+LnANmax
      SiteAmpCalLog[i] <- LnANmax
    }
    SigmaT[i]=sqrt(Sigma[i]^2+tauEvent[i]^2)
    sigmaSS[i]=sqrt(SigmaSC[i,SiteClass]^2+tauEvent[i]^2)

  }

  # Compute Sa and sigma with pre-defined period
  if (length(T) == 1 & T[1] == 1000) {

    T <- period

  } else {

    SpeCalLog <- approx(x = period, y = SpeCalLog, xout = T, rule = 2)$y

    RockSpeCalLog <- approx(x = period, y = RockSpeCalLog, xout = T, rule = 2)$y

    SiteAmpCalLog <- approx(x = period, y = SiteAmpCalLog, xout = T, rule = 2)$y

    SigmaT <- approx(x = period, y = SigmaT, xout = T, rule = 2)$y

    sigmaSS <- approx(x = period, y = sigmaSS, xout = T, rule = 2)$y
  }

  # output
  res <- list()

  res$Sa <- exp(SpeCalLog)

  res$RockSa <- exp(RockSpeCalLog)

  res$SiteAmp <- SiteAmpCalLog

  res$sigma <- SigmaT

  res$phiSS <- sigmaSS

  res$period <- T

  return(res)
}


#' The GMPE of Zea 2016 for slab event
#'
#' This function calculates the ground motion median values and standard deviations
#' @param mag Moment magnitude, a numeric value
#' @param T Period (sec) (range should be between 0.005 s and 5 s).
#' Use 1000 (by default) for output the array of Sa with original Zhao et al per-defined periods
#' @param fDepth Depth(km) to the top of ruptured plane (km)
#' @param dist The closest distance (km) to rupture plane when is available,
#' otherwise hypocentral distance.
#' @param DistV The distance of volcanic path (km)
#' @param SiteClass The site class: 1 for SC I, 2 for SC II, 3 for SC III 4 for SC IV sites. The relationship
#' between site class and Vs30 is given by Zhao et al., 2016 (https://doi.org/10.1785/0120150056).
#' SC I: Vs30 > 600 m/s; SC II: 300 < Vs30 =< 600 m/s; SC III: 200 < Vs30 =< 300 m/s; SC IV: Vs30 <= 200 m/s
#' @return A list of six elements is returned: Sa - median spectral acceleration prediction (in g);
#' RockSa - median spectral acceleration prediction on Rock reference site (in g);
#' SiteAmp - Site amplication response spectral (subtraction in log(g)); sigma - totla standard deviation;
#' phiSS - single station standard deviation; period - the corresponding oscillator periods
#' @examples zea_2016_Sub_Slab(mag = 7.5, T = 1000, fDepth = 0.0, dist = 40.0,
#' DistV = 0.0, SiteClass = 4)
#'
#' zea_2016_Sub_Slab(mag = 7.5, T = c(0.025, 0.5), fDepth = 0.0, dist = 40.0,
#' DistV = 0.0, SiteClass = 4)
#' @references Zhao, et al. (2016). Ground‐Motion Prediction Equations for Subduction Slab Earthquakes in
#' Japan Using Site Class and Simple Geometric Attenuation Functions.
#' Bulletin of the Seismological Society of America. 106(4): 1535-1551.
#' @export
#' @importFrom stats approx
zea_2016_Sub_Slab <- function(mag, T = 1000, fDepth, dist, DistV, SiteClass){

  period <- c(0.005,0.01,0.02,0.03,0.04,0.05,0.06,
              0.07,0.08,0.09,0.1,0.12,0.14,0.15,0.16,0.18,0.2,0.25,0.3,
              0.35,0.4,0.45,0.5,0.6,0.7,0.8,0.9,1.0,1.25,1.5,2.0,2.5,
              3.0,3.5,4.0,4.5,5.0)

  g <- 980.7

  Numperiod <- length(period)

  SpeCalLog <- rep(0, Numperiod)

  RockSpeCalLog <- rep(0, Numperiod)

  SiteAmpCalLog <- rep(0, Numperiod)

  SigmaT <- rep(0, Numperiod)

  sigmaSS <- rep(0, Numperiod)

  # model coefficients
  alfa <- c(-5.30119,-5.28844,-5.27568,-5.26822,
            -5.26293,-5.25882,-5.25547,-5.25263,-5.25017,-5.24801,-5.24607,
            -5.24271,-5.23988,-5.23861,-5.23742,-5.23525,-5.23331,-5.22921,
            -5.22585,-5.22302,-5.22056,-5.21839,-5.21645,-5.21310,-5.21026,
            -5.20781,-5.20564,-5.20370,-5.19959,-5.19624,-5.19095,-5.18684,
            -5.18349,-5.18065,-5.17819,-5.17602,-5.17409)
  beta <- c(1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,
            1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151,1.151)
  cSL <- c(1.44758,1.45400,1.46625,1.49246,1.50129,
           1.51051,1.51380,1.51111,1.50406,1.49423,1.48300,1.45559,
           1.44277,1.43314,1.43253,1.43710,1.44781,1.48260,1.51881,
           1.55291,1.58443,1.61360,1.64075,1.69020,1.73450,1.77474,
           1.81162,1.84561,1.92015,1.98274,2.08214,2.15841,2.22046,
           2.27406,2.32307,2.37009,2.37009)
  cSL2 <- c(0.37625,0.38099,0.39101,0.41976,0.45746,
            0.48601,0.50311,0.50704,0.50004,0.48071,0.45759,0.41355,
            0.37828,0.36308,0.34919,0.32464,0.30358,0.26174,0.23036,
            0.20580,0.18597,0.16960,0.15585,0.13405,0.11757,0.10476,
            0.09458,0.08636,0.07173,0.06258,0.05327,0.05036,0.04536,
            0.04536,0.04536,0.04536,0.04536)
  dSL <- c(0.42646,0.42075,0.40055,0.36433,0.32072,
           0.30000,0.31147,0.32673,0.34289,0.35921,0.37000,0.40606,
           0.43450,0.45000,0.46055,0.48439,0.509,0.555,0.593,0.625,
           0.652,0.675,0.695,0.729,0.756,0.778,0.796,0.812,0.841,
           0.861,0.884,0.9,0.9,0.9,0.9,0.9,0.9)
  bSLH <- c(0.01826,0.01826,0.01826,0.01826,0.01826,
            0.01826,0.01826,0.01826,0.01826,0.01826,0.01826,0.01826,
            0.01826,0.01826,0.01826,0.01826,0.01826,0.01826,0.01826,
            0.01826,0.01826,0.01826,0.01826,0.01826,0.01826,0.01826,
            0.01826,0.01826,0.01808,0.01786,0.01718,0.01628,0.01549,
            0.01489,0.01458,0.01459,0.01459)
  gSL <- c(-1.98471,-1.96360,-1.91839,-1.89271,
           -1.87260,-1.85351,-1.83395,-1.81345,-1.79189,-1.76931,-1.74581,
           -1.73746,-1.74463,-1.74972,-1.76259,-1.78989,-1.82110,-1.90412,
           -1.98439,-2.05756,-2.12282,-2.18047,-2.23118,-2.31475,-2.37885,
           -2.42769,-2.46450,-2.49170,-2.52758,-2.53359,-2.49565,-2.42623,
           -2.34726,-2.27002,-2.19947,-2.12528,-2.02646)
  gSLL <- c(1.12071,1.03278,0.94715,0.93420,0.97168,
            1.01492,1.06854,1.13401,1.20364,1.25808,1.30112,1.39137,
            1.47084,1.50784,1.54326,1.60985,1.67146,1.80738,1.92242,
            2.02102,2.10642,2.18097,2.24651,2.35602,2.44331,2.51391,
            2.57166,2.61931,2.70638,2.76244,2.82205,2.84475,2.84988,
            2.84667,2.83992,2.82802,2.82521)
  eSLV <- c(-0.01499,-0.01503,-0.01517,-0.01567,
            -0.01616,-0.01676,-0.01722,-0.01752,-0.01768,-0.01772,-0.01768,
            -0.01742,-0.01700,-0.01676,-0.01649,-0.01594,-0.01537,-0.01395,
            -0.01261,-0.01139,-0.01029,-0.00931,-0.00843,-0.00694,-0.00574,
            -0.00477,-0.00398,-0.00333,-0.00215,-0.00142,-0.00067,-0.00039,
            -0.00030,-0.00026,-0.00021,-0.00021,-0.00021)
  eSL <- c(-0.00340,-0.00331,-0.00345,-0.00391,
           -0.00454,-0.00510,-0.00552,-0.00588,-0.00615,-0.00635,-0.00652,
           -0.00660,-0.00652,-0.00647,-0.00636,-0.00614,-0.00590,-0.00526,
           -0.00468,-0.00415,-0.00369,-0.00327,-0.00290,-0.00227,-0.00178,
           -0.00139,-0.00109,-0.00086,-0.00052,-0.00043,-0.00070,-0.00127,
           -0.00198,-0.00271,-0.00341,-0.00421,-0.00500)
  eSLH <- c(-0.00050,-0.00050,-0.00050,-0.00050,
            -0.00050,-0.00050,-0.00050,-0.00049,-0.00048,-0.00048,-0.00048,
            -0.00049,-0.00051,-0.00052,-0.00053,-0.00056,-0.00059,-0.00067,
            -0.00075,-0.00083,-0.00091,-0.00099,-0.00107,-0.00124,-0.00139,
            -0.00154,-0.00166,-0.00178,-0.00199,-0.00213,-0.00225,-0.00219,
            -0.00207,-0.00193,-0.00180,-0.00170,-0.00158)
  gammaSL <- c(-9.87956,-9.51269,-9.26626,-9.33150,
               -9.50798,-9.72858,-9.96628,-10.22583,-10.55111,-10.80721,
               -11.02190,-11.36530,-11.73039,-11.88013,-12.05637,-12.42044,
               -12.78542,-13.63537,-14.38086,-15.03511,-15.61599,-16.13830,
               -16.61324,-17.45298,-18.18095,-18.82499,-19.40313,-19.92766,
               -21.05818,-21.99633,-23.48839,-24.64741,-25.59713,-26.40997,
               -27.13181,-27.79299,-28.31346)
  Imav <- c(3.73,3.07,2.76,3.02)
  phaiSC <- c(3.50,3.00,2.50,3.00)
  gammaSC <- c(0.8,1.00,0.90,0.60)
  SCTerm <- c(0.23200,0.22886,0.21825,0.18737,
              0.12332,0.07207,0.02701,-0.00621,0.01565,0.05089,0.09560,
              0.20037,0.30372,0.34284,0.37400,0.42700,0.46297,0.50856,
              0.50776,0.49710,0.48065,0.46159,0.44224,0.40537,0.37342,
              0.34623,0.32364,0.30479,0.27026,0.24831,0.22529,0.21540,
              0.21154,0.20976,0.20875,0.20774,0.20672, # start the second column
              0.14371,0.13978,0.12600,0.06164,
              -0.01705,-0.06331,-0.10103,-0.14684,-0.14479,-0.12666,-0.09319,
              -0.00878,0.08926,0.13602,0.17751,0.25309,0.32005,0.45304,
              0.54875,0.61713,0.66634,0.70111,0.72558,0.75294,0.76245,
              0.76119,0.75384,0.74279,0.70833,0.67256,0.61067,0.56403,
              0.52612,0.49766,0.47685,0.46223,0.45267, # start the third column
              0.14704,0.13284,0.14431,0.06600,
              -0.01714,-0.07312,-0.11955,-0.16006,-0.12434,-0.07293,-0.01458,
              0.08252,0.17151,0.20932,0.24117,0.29896,0.34591,0.44231,
              0.51782,0.57596,0.62239,0.65976,0.69066,0.73796,0.77226,
              0.79736,0.81616,0.83009,0.85036,0.85732,0.84991,0.82757,
              0.79911,0.76782,0.73594,0.70407,0.67220)
  SCTerm <- matrix(SCTerm, nrow = Numperiod)
  Sigma <- c(0.587,0.587,0.587,0.588,0.600,0.607,0.623,
             0.638,0.651,0.663,0.674,0.690,0.692,0.696,0.697,0.704,0.713,
             0.711,0.683,0.665,0.657,0.647,0.640,0.633,0.633,0.636,0.636,
             0.637,0.635,0.645,0.633,0.608,0.582,0.562,0.540,0.525,0.522)
  tauEvent <- c(0.457,0.459,0.465,0.480,0.520,
                0.555,0.584,0.598,0.598,0.585,0.567,0.534,0.504,0.486,0.465,
                0.430,0.406,0.385,0.365,0.373,0.383,0.391,0.403,0.412,0.432,
                0.436,0.437,0.436,0.444,0.448,0.424,0.413,0.407,0.394,0.381,
                0.365,0.378)
  LnSCIAm <- c(0.323,0.205,0.083,0.041,0.034,0.046,
               0.069,0.098,0.132,0.169,0.208,0.288,0.37,0.412,0.453,0.535,
               0.606,0.670,0.710,0.71893,0.70561,0.69289,0.68075,0.65804,
               0.63715,0.61781,0.59979,0.5829,0.54475,0.51118,0.45382,
               0.405624,0.36383,0.32682,0.29350,0.263,0.235)
  LnAmax1D <- c(0.65022,0.65181,0.65362,0.65467,
                0.65285,0.67264,0.69966,0.71713,0.71603,0.72561,0.74200,
                0.76236,0.75215,0.73819,0.71911,0.65408,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,0.58395,
                0.58395,0.58395,0.58395,0.58395,0.58395, # start the second column
                0.70973,0.70679,0.69465,0.68755,
                0.69892,0.70137,0.72445,0.74343,0.78598,0.79721,0.81668,
                0.84523,0.78296,0.79480,0.80861,0.84331,0.87770,0.93767,
                0.95,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                1.0,1.0,1.0,1.0,1.0, # start the third column
                0.64434,0.64624,0.63865,0.63421,
                0.60604,0.61716,0.63797,0.65437,0.68019,0.70889,0.71881,
                0.72581,0.74525,0.76103,0.76813,0.75690,0.71785,0.65470,
                0.69619,0.77907,0.82776,0.87645,0.92514,0.97383,1.02252,
                1.07122,1.11991,1.16860,1.21729,1.26598,1.31467,1.36336,
                1.41205,1.46075,1.50944,1.55813,1.60682, # start the forth column
                0.40428,0.40428,0.38789,0.37830,
                0.31737,0.30934,0.32530,0.35412,0.39282,0.42184,0.43736,
                0.47208,0.51278,0.53432,0.55022,0.57279,0.59674,0.61136,
                0.62638,0.63012,0.64773,0.64152,0.65582,0.68668,0.70560,
                0.71429,0.70388,0.67813,0.61119,0.54736,0.45944,0.40846,
                0.36421,0.32984,0.30912,0.29251,0.54736)
  LnAmax1D <- matrix(LnAmax1D, nrow = Numperiod)
  Src1D <- c(8.42900,8.09039,6.99211,6.34965,
             4.88339,5.04257,6.27132,7.66694,9.03436,11.25066,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739,14.81739,14.81739,14.81739,14.81739,
             14.81739,14.81739, # start the second column, nearly all of these large numbers of 14.81739, are not used
             1.91368,1.88256,1.77861,1.71781,
             2.05234,2.38713,2.83399,3.29447,3.99091,4.46576,5.04561,
             5.89960,5.05353,5.20490,5.38694,5.87165,6.57391,8.50000,
             10.6703,10.6703,10.6703,10.6703,10.6703,10.6703,10.6703,
             10.6703,10.6703,10.6703,10.6703,10.6703,10.6703,10.6703,
             10.6703,10.6703,10.6703,10.6703,10.6703, # start the third column, nearly all of these large numbers of 10.6703, are not used
             1.11714,1.11444,1.12437,1.13017,
             1.15080,1.23971,1.34819,1.45181,1.58315,1.73292,1.84134,
             2.03029,2.28133,2.44413,2.58017,2.74161,2.82587,2.71893,
             2.41759,2.30375,2.23625,2.21678,2.24338,2.80535,6.65839,
             30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0,30.0, # start the forth column, nearly all of these large numbers of 30.0, are not used
             0.83644,0.83644,0.83000,0.82624,
             0.76758,0.78632,0.83775,0.92616,1.02228,1.11802,1.16578,
             1.28551,1.39808,1.44327,1.47177,1.54694,1.64401,1.79013,
             1.82345,1.79037,1.76844,1.67539,1.62539,1.52453,1.39724,
             1.32029,1.26637,1.22680,1.22065,1.31805,2.12485,14.38181,
             14.38181,14.38181,14.38181,14.38181,14.38181) # nearly all of these large numbers of 14.38181 are not used
  Src1D <- matrix(Src1D, nrow = Numperiod)
  fsr <- c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
           1.0,1.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.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.0, # start the second column, When fsr is zero only linear site terms will be used
           1.0,1.0,1.0,1.0,1.006,0.851,0.803,
           0.918,1.062,1.106,1.071,0.9515,0.672,0.631,0.6,0.571,0.565,
           0.601,0.579,0.679,0.655,0.615,0.55,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,0.0, # start the third column
           1.0,1.0,1.0,1.0,1.0,1.0,1.044,0.975,
           0.964,0.98,0.97,1.022,0.889,0.861,0.831,0.748,0.65,0.479,0.449,
           0.482,0.499,0.515,0.53,0.53,0.499,0.369,0.3,0.2,0.0,0.0,0.0,
           0.0,0.0,0.0,0.0,0.0,0.0, # start the forth column
           1.0,1.0,1.05,0.58,0.482,0.472,0.506,
           0.587,0.683,0.782,0.823,1.029,0.991,0.983,0.973,0.979,1.006,
           1.027,1.021,1.003,1.01,0.985,0.99,1.006,1.0,1.0,0.96,0.904,
           0.738,0.535,0.358,0.0,0.0,0.0,0.0,0.0,0.0) # When fsr is zero only linear site terms will be used
  fsr <- matrix(fsr, nrow = Numperiod)
  SigmaSC <- c(0.398,0.397,0.395,0.389,0.387,0.387,
               0.397,0.403,0.413,0.422,0.429,0.440,0.441,0.450,0.452,0.454,
               0.462,0.474,0.472,0.468,0.457,0.450,0.445,0.448,0.440,0.441,
               0.435,0.427,0.413,0.416,0.409,0.398,0.390,0.386,0.377,0.360,
               0.361, # start the second column
               0.417,0.417,0.417,0.418,0.420,0.422,
               0.416,0.412,0.412,0.412,0.418,0.429,0.435,0.440,0.444,0.447,
               0.450,0.470,0.475,0.478,0.484,0.477,0.470,0.460,0.460,0.460,
               0.456,0.448,0.442,0.448,0.437,0.429,0.423,0.410,0.410,0.412,
               0.447, # start the third column
               0.409,0.409,0.408,0.409,0.413,0.409,
               0.401,0.394,0.389,0.391,0.394,0.428,0.433,0.420,0.424,0.446,
               0.441,0.459,0.437,0.444,0.453,0.477,0.472,0.459,0.461,0.457,
               0.449,0.442,0.428,0.418,0.406,0.387,0.369,0.376,0.362,0.368,
               0.381, # start the forth column
               0.415,0.415,0.416,0.417,0.420,0.422,
               0.423,0.420,0.423,0.427,0.426,0.445,0.442,0.440,0.437,0.436,
               0.432,0.432,0.432,0.432,0.416,0.408,0.409,0.407,0.403,0.407,
               0.408,0.408,0.412,0.418,0.413,0.414,0.412,0.402,0.394,0.385,
               0.380)
  SigmaSC <- matrix(SigmaSC, nrow = Numperiod)
  tauSC <- c(0.511,0.517,0.518,0.537,0.572,0.586,
             0.613,0.633,0.643,0.631,0.626,0.614,0.613,0.602,0.599,0.599,
             0.590,0.555,0.532,0.505,0.479,0.461,0.451,0.436,0.430,0.427,
             0.431,0.440,0.440,0.449,0.441,0.425,0.396,0.386,0.373,0.364,
             0.359, # start the second column
             0.449,0.450,0.449,0.449,0.456,0.479,
             0.507,0.532,0.555,0.568,0.566,0.570,0.597,0.599,0.597,0.601,
             0.596,0.601,0.557,0.518,0.488,0.479,0.471,0.467,0.464,0.456,
             0.466,0.466,0.476,0.473,0.462,0.427,0.407,0.401,0.386,0.373,
             0.322, # start the third column
             0.431,0.431,0.431,0.430,0.428,0.429,
             0.449,0.456,0.455,0.475,0.507,0.548,0.507,0.490,0.494,0.500,
             0.491,0.450,0.499,0.532,0.546,0.529,0.495,0.470,0.473,0.457,
             0.439,0.442,0.432,0.461,0.472,0.484,0.447,0.434,0.406,0.383,
             0.316, # start the forth column
             0.422,0.418,0.425,0.422,0.431,0.439,
             0.445,0.473,0.491,0.518,0.559,0.576,0.561,0.555,0.552,0.553,
             0.562,0.507,0.489,0.463,0.466,0.461,0.467,0.434,0.430,0.454,
             0.456,0.462,0.444,0.443,0.439,0.427,0.431,0.413,0.391,0.366,
             0.324)
  tauSC <- matrix(tauSC, nrow = Numperiod)
  Alfa1D=2.0
  Beta1D=0.6
  magc=7.1
  magcc=6.3
  Imin=1.0
  Imax=12.0
  DepthC=50.0
  Hcut=100.0

  # calculation
  for(i in seq(1, Numperiod)){
    # Bilinear magnitude scalling term and magnitude dependent dist cons
    ############ source effect ##############
    if(mag <= magc){
      r0 = dist+exp(alfa[i]+beta[i]*mag)
      SpeCalLog[i]=cSL[i]*mag+cSL2[i]*(mag-magcc)^2
    }else{
      r0 = dist + exp(alfa[i]+beta[i]*magc)
      SpeCalLog[i]=cSL[i]*magc+cSL2[i]*(magc-magcc)^2+dSL[i]*(mag-magc)
    }
    TmpD=SpeCalLog[i]+bSLH[i]*fDepth
    TmpF=gSL[i]*log(r0)
    TmpG=gSLL[i]*log(dist+200.0)
    if(fDepth >= Hcut){
      SpeCalLog[i]=SpeCalLog[i]+bSLH[i]*Hcut
    }else{
      SpeCalLog[i]=SpeCalLog[i]+bSLH[i]*fDepth
    }
    ########################################
    #################### add path effects ####################
    SpeCalLog[i]=SpeCalLog[i]+gSL[i]*log(r0)+eSL[i]*dist+
      eSLV[i]*DistV+gammaSL[i]+gSLL[i]*log(dist+200.0)
    if(fDepth >= DepthC) SpeCalLog[i]=SpeCalLog[i]+eSLH[i]*(0.02*fDepth-1.0)*dist
    ########################################
    ############ add site amplification effects ####################
    LnANmax=LnSCIAm[i]
    if(SiteClass >= 2) LnANmax=LnANmax+SCTerm[i,SiteClass-1]
    RockSpeCalLog[i] <- SpeCalLog[i]-LnSCIAm[i]
    nSReff=exp(RockSpeCalLog[i])/phaiSC[SiteClass]
    nSReffc=Src1D[i,SiteClass]/phaiSC[SiteClass]
    if(Imav[SiteClass] > Imax){
      nSReff=nSReff*(1.0+gammaSC[SiteClass]*(Imax-Imin))
      nSReffc=nSReffc*(1.0+gammaSC[SiteClass]*(Imax-Imin))
    }else if(Imav[SiteClass] > Imin){
      nSReff=nSReff*(1.0+gammaSC[SiteClass]*(Imav[SiteClass]-Imin))
      nSReffc=nSReffc*(1.0+gammaSC[SiteClass]*(Imav[SiteClass]-Imin))
    }
    if(SiteClass == 1){
      LnSF=LnSCIAm[i]-LnAmax1D[i,SiteClass]
    }else{
      LnSF=LnSCIAm[i]+SCTerm[i,SiteClass-1]-LnAmax1D[i,SiteClass]
    }
    if(exp(LnANmax) < 1.25){
      ca=LnAmax1D[i,SiteClass]/(log(Beta1D)-log(nSReffc^Alfa1D+Beta1D))
      cb=-ca*log(nSReffc^Alfa1D+Beta1D)
      SNC=exp((ca*(Alfa1D-1.0)*log(Beta1D)*log(10.0*Beta1D)-
                 log(10.0)*(cb+LnSF))/(ca*(Alfa1D*log(10.0*Beta1D)-log(Beta1D))))
    }else{
      SNC=(exp((LnANmax*log(nSReffc^Alfa1D+Beta1D)-LnSF*
                  log(Beta1D))/LnAmax1D[i,SiteClass])-Beta1D)^(1.0/Alfa1D)
    }
    SMR=nSReff*SNC/nSReffc*fsr[i,SiteClass]
    if(SMR != 0){
      SpeCalLog[i]=RockSpeCalLog[i]+LnANmax-LnAmax1D[i,SiteClass]*
        (log(SMR^Alfa1D+Beta1D)-log(Beta1D))/(log(nSReffc^Alfa1D+Beta1D)-log(Beta1D))
      SiteAmpCalLog[i] <- LnANmax-LnAmax1D[i,SiteClass]*
        (log(SMR^Alfa1D+Beta1D)-log(Beta1D))/(log(nSReffc^Alfa1D+Beta1D)-log(Beta1D))
    }else{
      SpeCalLog[i]=RockSpeCalLog[i]+LnANmax
      SiteAmpCalLog[i] <- LnANmax
    }
    SigmaT[i]=sqrt(Sigma[i]^2+tauEvent[i]^2)
    sigmaSS[i]=sqrt(SigmaSC[i,SiteClass]^2+tauEvent[i]^2)
  }

  # Compute Sa and sigma with pre-defined period
  if (length(T) == 1 & T[1] == 1000) {

    T <- period

  } else {

    SpeCalLog <- approx(x = period, y = SpeCalLog, xout = T, rule = 2)$y

    RockSpeCalLog <- approx(x = period, y = RockSpeCalLog, xout = T, rule = 2)$y

    SiteAmpCalLog <- approx(x = period, y = SiteAmpCalLog, xout = T, rule = 2)$y

    SigmaT <- approx(x = period, y = SigmaT, xout = T, rule = 2)$y

    sigmaSS <- approx(x = period, y = sigmaSS, xout = T, rule = 2)$y
  }

  # output
  res <- list()

  res$Sa <- exp(SpeCalLog)

  res$RockSa <- exp(RockSpeCalLog)

  res$SiteAmp <- SiteAmpCalLog

  res$sigma <- SigmaT

  res$phiSS <- sigmaSS

  res$period <- T

  return(res)
}