/*
 * SP_1996.c
 *
 * Code generation for function 'SP_1996'
 *
 * C source code generated on: Wed Aug 26 14:59:33 2015
 *
 */

/* Include files */
#include "rt_nonfinite.h"
#include "Select_Ground_Motions.h"
#include "SP_1996.h"
#include "interp1.h"
#include "Select_Ground_Motions_data.h"

/* Variable Definitions */
static emlrtRSInfo ad_emlrtRSI = { 91, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

static emlrtRSInfo bd_emlrtRSI = { 87, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

static emlrtRSInfo cd_emlrtRSI = { 86, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

static emlrtRSInfo dd_emlrtRSI = { 85, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

static emlrtRSInfo ed_emlrtRSI = { 84, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

static emlrtRSInfo fd_emlrtRSI = { 82, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

static emlrtRSInfo gd_emlrtRSI = { 81, "SP_1996",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/SP_1996.m"
};

/* Function Definitions */
void SP_1996(real_T T, real_T M, real_T R, real_T *Sa, real_T *sigma)
{
  real_T f;
  int32_T i3;
  static const real_T x[14] = { 2.0844908830972888, 2.0844908830972888,
    2.0844908830972888, 2.0653801558105296, 2.0323570109362219,
    2.0090928126087282, 1.9498445997580451, 1.9054607179632472,
    1.8197008586099834, 1.7139573075084253, 1.6595869074375607,
    1.6143585568264862, 1.5848931924611136, 1.5667510701081491 };

  real_T sigmaSP[14];
  static const real_T aSP[14] = { -2.5, -2.25, -1.9, -1.647, -1.28, -1.0, -0.595,
    -0.281, 0.1, 0.296, 0.222, -0.019, -0.312, -0.817 };

  static const real_T fSP[14] = { 0.25, 0.33, 0.5, 0.67, 1.0, 1.33, 2.0, 2.5,
    3.33, 5.0, 6.67, 10.0, 15.0, 25.0 };

  real_T a;
  static const real_T bSP[14] = { 0.725, 0.715, 0.687, 0.66, 0.612, 0.57, 0.5,
    0.445, 0.377, 0.323, 0.31, 0.304, 0.304, 0.33 };

  real_T b;
  static const real_T e1SP[14] = { 0.0, 0.0, 0.0, 0.01, 0.05, 0.12, 0.23, 0.222,
    0.185, 0.161, 0.161, 0.161, 0.161, 0.161 };

  real_T e1;
  static const real_T e2SP[14] = { 0.1, 0.108, 0.15, 0.175, 0.208, 0.19, 0.124,
    0.078, 0.02, 0.0, 0.0, 0.0, 0.0, 0.0 };

  real_T e2;
  static const real_T hSP[14] = { 2.6, 3.0, 3.6, 4.0, 4.4, 4.7, 5.0, 5.2, 5.4,
    5.7, 5.9, 6.2, 6.3, 4.7 };

  real_T h;

  /*  Sabetta & Pugliese ground motion prediction model, 1996 */
  /*  Marko Brozovi� */
  /*  2.8.2011 */
  /*  */
  /*  Useful for 4.6<=M<=6.8 and R<=100 km */
  /*  Only for larger horizontal component */
  /*  */
  /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
  /*  Input variables */
  /*  T = period (sec) */
  /*  M = magnitude */
  /*  R = fault distance (km) */
  /*  Site Geology */
  /*  S1: 1-shallow alluvium, 0-otherwise; 1-plitva naplavina H<20m, 0-skala(400  < stri�na hitrost < 800m/s), 0-skala(stri�na hitrost > 800m/s) */
  /*  S2: 1-deep alluvium, 0-otherwise 1-globoka naplavina H>20m, 0-skala(400  < stri�na hitrost < 800m/s), 0-skala(stri�na hitrost > 800m/s) */
  /*  coef: 1-fault distance, raw; 2-fault distance, smooth; 3-epicentral */
  /*        distance, raw; 4-epicentral distance, smooth; */
  /*  */
  /*  Output variables */
  /*  Sa: median spectral acceleration prediction */
  /*  sigma: logarithmic standard deviation of spectral acceleration */
  /*           prediction */
  /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
  /*  Period and frequency */
  if (T < 0.04) {
    T = 0.04;
  }

  /*  (s) */
  f = 1.0 / T;

  /*  (Hz) */
  /*  Regression coefficients */
  /*  epicentral distance, smooth coefficients (Table 2) */
  /*  Frequency (Hz) */
  /*  Constant term a */
  /*  Magnitude coefficient b */
  /*  Distance coefficient c */
  /*  Site coefficient e1 */
  /*  Site coefficient e2 */
  /*  h */
  for (i3 = 0; i3 < 14; i3++) {
    sigmaSP[i3] = muDoubleScalarLog(x[i3]);
  }

  /*  sigma */
  /*  Interpolation of coefficients */
  emlrtPushRtStackR2012b(&gd_emlrtRSI, emlrtRootTLSGlobal);
  a = interp1(fSP, aSP, f);
  emlrtPopRtStackR2012b(&gd_emlrtRSI, emlrtRootTLSGlobal);
  emlrtPushRtStackR2012b(&fd_emlrtRSI, emlrtRootTLSGlobal);
  b = interp1(fSP, bSP, f);
  emlrtPopRtStackR2012b(&fd_emlrtRSI, emlrtRootTLSGlobal);
  emlrtPushRtStackR2012b(&ed_emlrtRSI, emlrtRootTLSGlobal);
  e1 = interp1(fSP, e1SP, f);
  emlrtPopRtStackR2012b(&ed_emlrtRSI, emlrtRootTLSGlobal);
  emlrtPushRtStackR2012b(&dd_emlrtRSI, emlrtRootTLSGlobal);
  e2 = interp1(fSP, e2SP, f);
  emlrtPopRtStackR2012b(&dd_emlrtRSI, emlrtRootTLSGlobal);
  emlrtPushRtStackR2012b(&cd_emlrtRSI, emlrtRootTLSGlobal);
  h = interp1(fSP, hSP, f);
  emlrtPopRtStackR2012b(&cd_emlrtRSI, emlrtRootTLSGlobal);
  emlrtPushRtStackR2012b(&bd_emlrtRSI, emlrtRootTLSGlobal);
  *sigma = interp1(fSP, sigmaSP, f);
  emlrtPopRtStackR2012b(&bd_emlrtRSI, emlrtRootTLSGlobal);

  /*  Calculation of spectral acceleration */
  /*  gravitational acceleration (m/s2) */
  emlrtPushRtStackR2012b(&ad_emlrtRSI, emlrtRootTLSGlobal);
  f = R * R + h * h;
  if (f < 0.0) {
    emlrtPushRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal);
    e_eml_error();
    emlrtPopRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal);
  }

  *Sa = 6.2831853071795862 / T / 100.0 / 9.81 * muDoubleScalarPower(10.0, (((a +
    b * M) + -muDoubleScalarLog10(muDoubleScalarSqrt(f))) + e1) + e2 * 0.0);
  emlrtPopRtStackR2012b(&ad_emlrtRSI, emlrtRootTLSGlobal);

  /*  (g) */
}

void e_eml_error(void)
{
  static char_T cv19[4][1] = { { 's' }, { 'q' }, { 'r' }, { 't' } };

  emlrtPushRtStackR2012b(&jd_emlrtRSI, emlrtRootTLSGlobal);
  emlrtErrorWithMessageIdR2012b(emlrtRootTLSGlobal, &ee_emlrtRTEI,
    "Coder:toolbox:ElFunDomainError", 3, 4, 4, cv19);
  emlrtPopRtStackR2012b(&jd_emlrtRSI, emlrtRootTLSGlobal);
}

/* End of code generation (SP_1996.c) */