/*
 * File: eps_bar.c
 *
 * MATLAB Coder version            : 2.8
 * C/C++ source code generated on  : 03-Aug-2016 03:45:00
 */

/* Include Files */
#include "rt_nonfinite.h"
#include "eps_bar.h"

/* Function Declarations */
static double interp1(const double varargin_1[14], const double varargin_2[14],
                      double varargin_3);
static double rt_powd_snf(double u0, double u1);

/* Function Definitions */

/*
 * Arguments    : const double varargin_1[14]
 *                const double varargin_2[14]
 *                double varargin_3
 * Return Type  : double
 */
static double interp1(const double varargin_1[14], const double varargin_2[14],
                      double varargin_3)
{
  double Vq;
  double y[14];
  double x[14];
  int low_i;
  int32_T exitg1;
  double r;
  int low_ip1;
  int high_i;
  int mid_i;
  for (low_i = 0; low_i < 14; low_i++) {
    y[low_i] = varargin_2[low_i];
    x[low_i] = varargin_1[low_i];
  }

  Vq = rtNaN;
  low_i = 1;
  do {
    exitg1 = 0;
    if (low_i < 15) {
      if (rtIsNaN(varargin_1[low_i - 1])) {
        exitg1 = 1;
      } else {
        low_i++;
      }
    } else {
      if (varargin_1[1] < varargin_1[0]) {
        for (low_i = 0; low_i < 7; low_i++) {
          r = x[low_i];
          x[low_i] = x[13 - low_i];
          x[13 - low_i] = r;
        }

        for (low_i = 0; low_i < 7; low_i++) {
          r = y[low_i];
          y[low_i] = y[13 - low_i];
          y[13 - low_i] = r;
        }
      }

      if (rtIsNaN(varargin_3) || (varargin_3 > x[13]) || (varargin_3 < x[0])) {
      } else {
        low_i = 1;
        low_ip1 = 2;
        high_i = 14;
        while (high_i > low_ip1) {
          mid_i = (low_i + high_i) >> 1;
          if (varargin_3 >= x[mid_i - 1]) {
            low_i = mid_i;
            low_ip1 = mid_i + 1;
          } else {
            high_i = mid_i;
          }
        }

        r = (varargin_3 - x[low_i - 1]) / (x[low_i] - x[low_i - 1]);
        if (r == 0.0) {
          Vq = y[low_i - 1];
        } else if (r == 1.0) {
          Vq = y[low_i];
        } else if (y[low_i - 1] == y[low_i]) {
          Vq = y[low_i - 1];
        } else {
          Vq = (1.0 - r) * y[low_i - 1] + r * y[low_i];
        }
      }

      exitg1 = 1;
    }
  } while (exitg1 == 0);

  return Vq;
}

/*
 * Arguments    : double u0
 *                double u1
 * Return Type  : double
 */
static double rt_powd_snf(double u0, double u1)
{
  double y;
  double d0;
  double d1;
  if (rtIsNaN(u0) || rtIsNaN(u1)) {
    y = rtNaN;
  } else {
    d0 = fabs(u0);
    d1 = fabs(u1);
    if (rtIsInf(u1)) {
      if (d0 == 1.0) {
        y = rtNaN;
      } else if (d0 > 1.0) {
        if (u1 > 0.0) {
          y = rtInf;
        } else {
          y = 0.0;
        }
      } else if (u1 > 0.0) {
        y = 0.0;
      } else {
        y = rtInf;
      }
    } else if (d1 == 0.0) {
      y = 1.0;
    } else if (d1 == 1.0) {
      if (u1 > 0.0) {
        y = u0;
      } else {
        y = 1.0 / u0;
      }
    } else if (u1 == 2.0) {
      y = u0 * u0;
    } else if ((u1 == 0.5) && (u0 >= 0.0)) {
      y = sqrt(u0);
    } else if ((u0 < 0.0) && (u1 > floor(u1))) {
      y = rtNaN;
    } else {
      y = pow(u0, u1);
    }
  }

  return y;
}

/*
 * EPS_BAR Summary of this function goes here
 *    Detailed explanation goes here
 * Arguments    : double agC
 *                double Rls
 *                double T1
 *                double M_mean
 *                double R_mean
 *                double ecSoilCode
 * Return Type  : double
 */
double eps_bar(double agC, double Rls, double T1, double M_mean, double R_mean,
               double ecSoilCode)
{
  int S1;
  int S2;
  double T;
  double f;
  double aSP[14];
  double bSP[14];
  double e1SP[14];
  double e2SP[14];
  double hSP[14];
  double sigmaSP[14];
  int i0;
  static const double dv0[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 double dv1[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 };

  static const double dv2[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 };

  static const double dv3[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 };

  static const double dv4[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 };

  static const double 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 };

  static const double 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 };

  double h;
  if (ecSoilCode == 1.0) {
    S1 = 0;

    /*  is there shallow alluvium */
    S2 = 0;

    /*  is there deep alluvium */
    /*  1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef */
  } else if (ecSoilCode == 2.0) {
    S1 = 0;

    /*  is there shallow alluvium */
    S2 = 1;

    /*  is there deep alluvium */
    /*  1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef   */
  } else if (ecSoilCode == 3.0) {
    S1 = 1;

    /*  is there shallow alluvium */
    S2 = 0;

    /*  is there deep alluvium */
    /*  1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef   */
  } else {
    S1 = 0;

    /*  is there shallow alluvium */
    S2 = 1;

    /*  is there deep alluvium */
    /*  1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef   */
  }

  T = T1;

  /*  Period and frequency */
  /*  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 */
  /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
  if (T1 < 0.04) {
    T = 0.04;
  }

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

  /*  (Hz) */
  /*  Regression coefficients */
  /*  epicentral distance, smooth coefficients (Table 2) */
  /*  Frequency (Hz) */
  for (i0 = 0; i0 < 14; i0++) {
    aSP[i0] = dv0[i0];

    /*  Constant term a */
    bSP[i0] = dv1[i0];

    /*  Magnitude coefficient b */
    /*  Distance coefficient c */
    e1SP[i0] = dv2[i0];

    /*  Site coefficient e1 */
    e2SP[i0] = dv3[i0];

    /*  Site coefficient e2 */
    hSP[i0] = dv4[i0];

    /*  h */
    sigmaSP[i0] = log(x[i0]);
  }

  /*  sigma */
  /*  Interpolation of coefficients */
  h = interp1(fSP, hSP, f);

  /*  Calculation of spectral acceleration */
  /*  gravitational acceleration (m/s2) */
  /*  (g) */
  return (log(agC / Rls) - log(6.2831853071795862 / T / 100.0 / 9.81 *
           rt_powd_snf(10.0, (((interp1(fSP, aSP, f) + interp1(fSP, bSP, f) *
    M_mean) + -log10(sqrt(R_mean * R_mean + h * h))) + interp1(fSP, e1SP, f) *
             (double)S1) + interp1(fSP, e2SP, f) * (double)S2))) / interp1(fSP,
    sigmaSP, f);
}

/*
 * Arguments    : void
 * Return Type  : void
 */
void eps_bar_initialize(void)
{
  rt_InitInfAndNaN(8U);
}

/*
 * Arguments    : void
 * Return Type  : void
 */
void eps_bar_terminate(void)
{
  /* (no terminate code required) */
}

/*
 * File trailer for eps_bar.c
 *
 * [EOF]
 */