/*
 * File: imc.c
 *
 * MATLAB Coder version            : 2.8
 * C/C++ source code generated on  : 03-Aug-2016 03:44:06
 */

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

/* Type Definitions */
#ifndef struct_emxArray__common
#define struct_emxArray__common

struct emxArray__common
{
  void *data;
  int *size;
  int allocatedSize;
  int numDimensions;
  boolean_T canFreeData;
};

#endif                                 /*struct_emxArray__common*/

#ifndef typedef_emxArray__common
#define typedef_emxArray__common

typedef struct emxArray__common emxArray__common;

#endif                                 /*typedef_emxArray__common*/

#ifndef struct_emxArray_boolean_T
#define struct_emxArray_boolean_T

struct emxArray_boolean_T
{
  boolean_T *data;
  int *size;
  int allocatedSize;
  int numDimensions;
  boolean_T canFreeData;
};

#endif                                 /*struct_emxArray_boolean_T*/

#ifndef typedef_emxArray_boolean_T
#define typedef_emxArray_boolean_T

typedef struct emxArray_boolean_T emxArray_boolean_T;

#endif                                 /*typedef_emxArray_boolean_T*/

/* Function Declarations */
static void b_abs(const emxArray_real_T *x, emxArray_real_T *y);
static void b_emxInit_real_T(emxArray_real_T **pEmxArray, int numDimensions);
static void b_exp(emxArray_real_T *x);
static void b_log(emxArray_real_T *x);
static void c_log(emxArray_real_T *x);
static void diff(const emxArray_real_T *x, emxArray_real_T *y);
static void emxEnsureCapacity(emxArray__common *emxArray, int oldNumel, int
  elementSize);
static void emxFree_boolean_T(emxArray_boolean_T **pEmxArray);
static void emxFree_real_T(emxArray_real_T **pEmxArray);
static void emxInit_boolean_T(emxArray_boolean_T **pEmxArray, int numDimensions);
static void emxInit_real_T(emxArray_real_T **pEmxArray, int numDimensions);
static void interp1(const emxArray_real_T *varargin_1, const emxArray_real_T
                    *varargin_2, const emxArray_real_T *varargin_3,
                    emxArray_real_T *Vq);
static double logncdf_(double x, double mu, double sigma);
static void rdivide(const emxArray_real_T *x, const emxArray_real_T *y,
                    emxArray_real_T *z);
static double rt_powd_snf(double u0, double u1);
static double sum(const emxArray_real_T *x);

/* Function Definitions */

/*
 * Arguments    : const emxArray_real_T *x
 *                emxArray_real_T *y
 * Return Type  : void
 */
static void b_abs(const emxArray_real_T *x, emxArray_real_T *y)
{
  unsigned int uv0[2];
  int k;
  for (k = 0; k < 2; k++) {
    uv0[k] = (unsigned int)x->size[k];
  }

  k = y->size[0] * y->size[1];
  y->size[0] = 1;
  y->size[1] = (int)uv0[1];
  emxEnsureCapacity((emxArray__common *)y, k, (int)sizeof(double));
  for (k = 0; k < x->size[1]; k++) {
    y->data[k] = fabs(x->data[k]);
  }
}

/*
 * Arguments    : emxArray_real_T **pEmxArray
 *                int numDimensions
 * Return Type  : void
 */
static void b_emxInit_real_T(emxArray_real_T **pEmxArray, int numDimensions)
{
  emxArray_real_T *emxArray;
  int i;
  *pEmxArray = (emxArray_real_T *)malloc(sizeof(emxArray_real_T));
  emxArray = *pEmxArray;
  emxArray->data = (double *)NULL;
  emxArray->numDimensions = numDimensions;
  emxArray->size = (int *)malloc((unsigned int)(sizeof(int) * numDimensions));
  emxArray->allocatedSize = 0;
  emxArray->canFreeData = true;
  for (i = 0; i < numDimensions; i++) {
    emxArray->size[i] = 0;
  }
}

/*
 * Arguments    : emxArray_real_T *x
 * Return Type  : void
 */
static void b_exp(emxArray_real_T *x)
{
  int i4;
  int k;
  i4 = x->size[1];
  for (k = 0; k < i4; k++) {
    x->data[k] = exp(x->data[k]);
  }
}

/*
 * Arguments    : emxArray_real_T *x
 * Return Type  : void
 */
static void b_log(emxArray_real_T *x)
{
  int i2;
  int k;
  i2 = x->size[0];
  for (k = 0; k < i2; k++) {
    x->data[k] = log(x->data[k]);
  }
}

/*
 * Arguments    : emxArray_real_T *x
 * Return Type  : void
 */
static void c_log(emxArray_real_T *x)
{
  int i3;
  int k;
  i3 = x->size[1];
  for (k = 0; k < i3; k++) {
    x->data[k] = log(x->data[k]);
  }
}

/*
 * Arguments    : const emxArray_real_T *x
 *                emxArray_real_T *y
 * Return Type  : void
 */
static void diff(const emxArray_real_T *x, emxArray_real_T *y)
{
  int orderForDim;
  int iyLead;
  double work_data_idx_0;
  int m;
  double tmp1;
  double tmp2;
  if (x->size[1] == 0) {
    orderForDim = y->size[0] * y->size[1];
    y->size[0] = 1;
    y->size[1] = 0;
    emxEnsureCapacity((emxArray__common *)y, orderForDim, (int)sizeof(double));
  } else {
    if (x->size[1] - 1 <= 1) {
      orderForDim = x->size[1] - 1;
    } else {
      orderForDim = 1;
    }

    if (orderForDim < 1) {
      orderForDim = y->size[0] * y->size[1];
      y->size[0] = 1;
      y->size[1] = 0;
      emxEnsureCapacity((emxArray__common *)y, orderForDim, (int)sizeof(double));
    } else {
      orderForDim = y->size[0] * y->size[1];
      y->size[0] = 1;
      y->size[1] = x->size[1] - 1;
      emxEnsureCapacity((emxArray__common *)y, orderForDim, (int)sizeof(double));
      orderForDim = x->size[1] - 1;
      if (!(orderForDim == 0)) {
        orderForDim = 1;
        iyLead = 0;
        work_data_idx_0 = x->data[0];
        for (m = 2; m <= x->size[1]; m++) {
          tmp1 = x->data[orderForDim];
          tmp2 = work_data_idx_0;
          work_data_idx_0 = tmp1;
          tmp1 -= tmp2;
          orderForDim++;
          y->data[iyLead] = tmp1;
          iyLead++;
        }
      }
    }
  }
}

/*
 * Arguments    : emxArray__common *emxArray
 *                int oldNumel
 *                int elementSize
 * Return Type  : void
 */
static void emxEnsureCapacity(emxArray__common *emxArray, int oldNumel, int
  elementSize)
{
  int newNumel;
  int i;
  void *newData;
  newNumel = 1;
  for (i = 0; i < emxArray->numDimensions; i++) {
    newNumel *= emxArray->size[i];
  }

  if (newNumel > emxArray->allocatedSize) {
    i = emxArray->allocatedSize;
    if (i < 16) {
      i = 16;
    }

    while (i < newNumel) {
      i <<= 1;
    }

    newData = calloc((unsigned int)i, (unsigned int)elementSize);
    if (emxArray->data != NULL) {
      memcpy(newData, emxArray->data, (unsigned int)(elementSize * oldNumel));
      if (emxArray->canFreeData) {
        free(emxArray->data);
      }
    }

    emxArray->data = newData;
    emxArray->allocatedSize = i;
    emxArray->canFreeData = true;
  }
}

/*
 * Arguments    : emxArray_boolean_T **pEmxArray
 * Return Type  : void
 */
static void emxFree_boolean_T(emxArray_boolean_T **pEmxArray)
{
  if (*pEmxArray != (emxArray_boolean_T *)NULL) {
    if (((*pEmxArray)->data != (boolean_T *)NULL) && (*pEmxArray)->canFreeData)
    {
      free((void *)(*pEmxArray)->data);
    }

    free((void *)(*pEmxArray)->size);
    free((void *)*pEmxArray);
    *pEmxArray = (emxArray_boolean_T *)NULL;
  }
}

/*
 * Arguments    : emxArray_real_T **pEmxArray
 * Return Type  : void
 */
static void emxFree_real_T(emxArray_real_T **pEmxArray)
{
  if (*pEmxArray != (emxArray_real_T *)NULL) {
    if (((*pEmxArray)->data != (double *)NULL) && (*pEmxArray)->canFreeData) {
      free((void *)(*pEmxArray)->data);
    }

    free((void *)(*pEmxArray)->size);
    free((void *)*pEmxArray);
    *pEmxArray = (emxArray_real_T *)NULL;
  }
}

/*
 * Arguments    : emxArray_boolean_T **pEmxArray
 *                int numDimensions
 * Return Type  : void
 */
static void emxInit_boolean_T(emxArray_boolean_T **pEmxArray, int numDimensions)
{
  emxArray_boolean_T *emxArray;
  int i;
  *pEmxArray = (emxArray_boolean_T *)malloc(sizeof(emxArray_boolean_T));
  emxArray = *pEmxArray;
  emxArray->data = (boolean_T *)NULL;
  emxArray->numDimensions = numDimensions;
  emxArray->size = (int *)malloc((unsigned int)(sizeof(int) * numDimensions));
  emxArray->allocatedSize = 0;
  emxArray->canFreeData = true;
  for (i = 0; i < numDimensions; i++) {
    emxArray->size[i] = 0;
  }
}

/*
 * Arguments    : emxArray_real_T **pEmxArray
 *                int numDimensions
 * Return Type  : void
 */
static void emxInit_real_T(emxArray_real_T **pEmxArray, int numDimensions)
{
  emxArray_real_T *emxArray;
  int i;
  *pEmxArray = (emxArray_real_T *)malloc(sizeof(emxArray_real_T));
  emxArray = *pEmxArray;
  emxArray->data = (double *)NULL;
  emxArray->numDimensions = numDimensions;
  emxArray->size = (int *)malloc((unsigned int)(sizeof(int) * numDimensions));
  emxArray->allocatedSize = 0;
  emxArray->canFreeData = true;
  for (i = 0; i < numDimensions; i++) {
    emxArray->size[i] = 0;
  }
}

/*
 * Arguments    : const emxArray_real_T *varargin_1
 *                const emxArray_real_T *varargin_2
 *                const emxArray_real_T *varargin_3
 *                emxArray_real_T *Vq
 * Return Type  : void
 */
static void interp1(const emxArray_real_T *varargin_1, const emxArray_real_T
                    *varargin_2, const emxArray_real_T *varargin_3,
                    emxArray_real_T *Vq)
{
  emxArray_real_T *y;
  int nd2;
  int n;
  emxArray_real_T *x;
  int nyrows;
  int nx;
  unsigned int outsize[2];
  int k;
  int32_T exitg1;
  double minx;
  double maxx;
  int mid_i;
  double r;
  b_emxInit_real_T(&y, 1);
  nd2 = y->size[0];
  y->size[0] = varargin_2->size[0];
  emxEnsureCapacity((emxArray__common *)y, nd2, (int)sizeof(double));
  n = varargin_2->size[0];
  for (nd2 = 0; nd2 < n; nd2++) {
    y->data[nd2] = varargin_2->data[nd2];
  }

  b_emxInit_real_T(&x, 1);
  nd2 = x->size[0];
  x->size[0] = varargin_1->size[0];
  emxEnsureCapacity((emxArray__common *)x, nd2, (int)sizeof(double));
  n = varargin_1->size[0];
  for (nd2 = 0; nd2 < n; nd2++) {
    x->data[nd2] = varargin_1->data[nd2];
  }

  nyrows = varargin_2->size[0] - 1;
  nx = varargin_1->size[0];
  for (nd2 = 0; nd2 < 2; nd2++) {
    outsize[nd2] = (unsigned int)varargin_3->size[nd2];
  }

  nd2 = Vq->size[0] * Vq->size[1];
  Vq->size[0] = 1;
  emxEnsureCapacity((emxArray__common *)Vq, nd2, (int)sizeof(double));
  nd2 = Vq->size[0] * Vq->size[1];
  Vq->size[1] = (int)outsize[1];
  emxEnsureCapacity((emxArray__common *)Vq, nd2, (int)sizeof(double));
  n = (int)outsize[1];
  for (nd2 = 0; nd2 < n; nd2++) {
    Vq->data[nd2] = 0.0;
  }

  if (varargin_3->size[1] == 0) {
  } else {
    k = 1;
    do {
      exitg1 = 0;
      if (k <= nx) {
        if (rtIsNaN(varargin_1->data[k - 1])) {
          exitg1 = 1;
        } else {
          k++;
        }
      } else {
        if (varargin_1->data[1] < varargin_1->data[0]) {
          nd2 = nx >> 1;
          for (n = 1; n <= nd2; n++) {
            minx = x->data[n - 1];
            x->data[n - 1] = x->data[nx - n];
            x->data[nx - n] = minx;
          }

          if ((!(varargin_2->size[0] == 0)) && (varargin_2->size[0] > 1)) {
            n = varargin_2->size[0];
            nd2 = varargin_2->size[0] >> 1;
            for (k = 1; k <= nd2; k++) {
              minx = y->data[k - 1];
              y->data[k - 1] = y->data[n - k];
              y->data[n - k] = minx;
            }
          }
        }

        minx = x->data[0];
        maxx = x->data[x->size[0] - 1];
        for (k = 0; k + 1 <= varargin_3->size[1]; k++) {
          if (rtIsNaN(varargin_3->data[k])) {
            Vq->data[k] = rtNaN;
          } else if (varargin_3->data[k] > maxx) {
            if (nyrows > 0) {
              Vq->data[k] = y->data[nyrows] + (varargin_3->data[k] - maxx) /
                (maxx - x->data[x->size[0] - 2]) * (y->data[nyrows] - y->
                data[nyrows - 1]);
            }
          } else if (varargin_3->data[k] < minx) {
            Vq->data[k] = y->data[0] + (varargin_3->data[k] - minx) / (x->data[1]
              - minx) * (y->data[1] - y->data[0]);
          } else {
            n = 1;
            nd2 = 2;
            nx = x->size[0];
            while (nx > nd2) {
              mid_i = (n >> 1) + (nx >> 1);
              if (((n & 1) == 1) && ((nx & 1) == 1)) {
                mid_i++;
              }

              if (varargin_3->data[k] >= x->data[mid_i - 1]) {
                n = mid_i;
                nd2 = mid_i + 1;
              } else {
                nx = mid_i;
              }
            }

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

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

  emxFree_real_T(&x);
  emxFree_real_T(&y);
}

/*
 * LOGNCDF Lognormal cumulative distribution function (cdf).
 *    P = LOGNCDF(X,MU,SIGMA) returns values at X of the lognormal cdf with
 *    distribution parameters MU and SIGMA.  MU and SIGMA are the mean and
 *    standard deviation, respectively, of the associated normal distribution.
 *    The size of P is the common size of X, MU and SIGMA.  A scalar input
 *    functions as a constant matrix of the same size as the other inputs.
 *
 *    Default values for MU and SIGMA are 0 and 1, respectively.
 *
 *    [P,PLO,PUP] = LOGNCDF(X,MU,SIGMA,PCOV,ALPHA) returns confidence bounds
 *    for P when the input parameters MU and SIGMA are estimates.  PCOV is a
 *    2-by-2 matrix containing the covariance matrix of the estimated parameters.
 *    ALPHA has a default value of 0.05, and specifies 100*(1-ALPHA)% confidence
 *    bounds.  PLO and PUP are arrays of the same size as P containing the lower
 *    and upper confidence bounds.
 *
 *    See also ERF, ERFC, LOGNFIT, LOGNINV, LOGNLIKE, LOGNPDF, LOGNRND, LOGNSTAT.
 * Arguments    : double x
 *                double mu
 *                double sigma
 * Return Type  : double
 */
static double logncdf_(double x, double mu, double sigma)
{
  double p;
  double z;
  int i;
  double s;
  double b_z;
  double absx;
  double R;
  int eint;

  /*    References: */
  /*       [1] Abramowitz, M. and Stegun, I.A. (1964) Handbook of Mathematical */
  /*           Functions, Dover, New York, 1046pp., sections 7.1, 26.2. */
  /*       [2] Evans, M., Hastings, N., and Peacock, B. (1993) Statistical */
  /*           Distributions, 2nd ed., Wiley, 170pp. */
  /*    Copyright 1993-2011 The MathWorks, Inc. */
  /*  More checking if we need to compute confidence bounds. */
  /*  Return NaN for out of range parameters. */
  z = sigma;
  i = 0;
  while (i < 1) {
    if (sigma <= 0.0) {
      z = rtNaN;
    }

    i = 1;
  }

  /*  Negative data would create complex values, which erfc cannot handle. */
  s = x;
  i = 0;
  while (i < 1) {
    if (x < 0.0) {
      s = 0.0;
    }

    i = 1;
  }

  /* try */
  /* catch ME */
  /*     error(message('stats:logncdf:InputSizeMismatch')); */
  /* end */
  /*  Use the complementary error function, rather than .5*(1+erf(z/sqrt(2))), */
  /*  to produce accurate near-zero results for large negative x. */
  b_z = -((log(s) - mu) / z) / 1.4142135623730951;

  /* ========================== COPYRIGHT NOTICE ============================ */
  /*  The algorithms for calculating ERF(X) and ERFC(X) are derived           */
  /*  from FDLIBM, which has the following notice:                            */
  /*                                                                          */
  /*  Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.       */
  /*                                                                          */
  /*  Developed at SunSoft, a Sun Microsystems, Inc. business.                */
  /*  Permission to use, copy, modify, and distribute this                    */
  /*  software is freely granted, provided that this notice                   */
  /*  is preserved.                                                           */
  /* =============================    END    ================================ */
  absx = fabs(b_z);
  if (rtIsNaN(b_z)) {
    z = b_z;
  } else if (rtIsInf(b_z)) {
    if (b_z < 0.0) {
      z = 2.0;
    } else {
      z = 0.0;
    }
  } else if (absx < 0.84375) {
    if (absx < 1.3877787807814457E-17) {
      z = 1.0 - b_z;
    } else {
      z = b_z * b_z;
      z = (0.12837916709551256 + z * (-0.3250421072470015 + z *
            (-0.02848174957559851 + z * (-0.0057702702964894416 + z *
              -2.3763016656650163E-5)))) / (1.0 + z * (0.39791722395915535 + z *
        (0.0650222499887673 + z * (0.0050813062818757656 + z *
        (0.00013249473800432164 + z * -3.9602282787753681E-6)))));
      if (b_z < 0.25) {
        z = 1.0 - (b_z + b_z * z);
      } else {
        z = 0.5 - (b_z * z + (b_z - 0.5));
      }
    }
  } else if (absx < 1.25) {
    s = -0.0023621185607526594 + (absx - 1.0) * (0.41485611868374833 + (absx -
      1.0) * (-0.37220787603570132 + (absx - 1.0) * (0.31834661990116175 + (absx
      - 1.0) * (-0.11089469428239668 + (absx - 1.0) * (0.035478304325618236 +
      (absx - 1.0) * -0.0021663755948687908)))));
    z = 1.0 + (absx - 1.0) * (0.10642088040084423 + (absx - 1.0) *
      (0.540397917702171 + (absx - 1.0) * (0.071828654414196266 + (absx - 1.0) *
      (0.12617121980876164 + (absx - 1.0) * (0.013637083912029051 + (absx - 1.0)
      * 0.011984499846799107)))));
    if (b_z >= 0.0) {
      z = 0.15493708848953247 - s / z;
    } else {
      z = 1.0 + (0.84506291151046753 + s / z);
    }
  } else if (b_z < -6.0) {
    z = 2.0;
  } else if (b_z >= 28.0) {
    z = 0.0;
  } else {
    s = 1.0 / (absx * absx);
    if (absx < 2.8571414947509766) {
      R = -0.0098649440348471482 + s * (-0.69385857270718176 + s *
        (-10.558626225323291 + s * (-62.375332450326006 + s *
        (-162.39666946257347 + s * (-184.60509290671104 + s * (-81.2874355063066
        + s * -9.8143293441691455))))));
      s = 1.0 + s * (19.651271667439257 + s * (137.65775414351904 + s *
        (434.56587747522923 + s * (645.38727173326788 + s * (429.00814002756783
        + s * (108.63500554177944 + s * (6.5702497703192817 + s *
        -0.0604244152148581)))))));
    } else {
      R = -0.0098649429247001 + s * (-0.799283237680523 + s *
        (-17.757954917754752 + s * (-160.63638485582192 + s *
        (-637.56644336838963 + s * (-1025.0951316110772 + s * -483.5191916086514)))));
      s = 1.0 + s * (30.338060743482458 + s * (325.79251299657392 + s *
        (1536.729586084437 + s * (3199.8582195085955 + s * (2553.0504064331644 +
        s * (474.52854120695537 + s * -22.440952446585818))))));
    }

    if ((!rtIsInf(absx)) && (!rtIsNaN(absx))) {
      z = frexp(absx, &eint);
    } else {
      z = absx;
      eint = 0;
    }

    z = floor(z * 2.097152E+6) / 2.097152E+6 * rt_powd_snf(2.0, eint);
    z = exp(-z * z - 0.5625) * exp((z - absx) * (z + absx) + R / s) / absx;
    if (b_z < 0.0) {
      z = 2.0 - z;
    }
  }

  p = 0.5 * z;

  /*  Compute confidence bounds if requested. */
  /* if nargout>=2 */
  /*    zvar = (pcov(1,1) + 2*pcov(1,2)*z + pcov(2,2)*z.^2) ./ (sigma.^2); */
  /*    if all(all(zvar >= 0) == 1) */
  /*        normz = -norminv(alpha/2); */
  /*        halfwidth = normz * sqrt(zvar); */
  /*        zlo = z - halfwidth; */
  /*        zup = z + halfwidth; */
  /*  */
  /*        plo = 0.5 * erfc(-zlo./sqrt(2)); */
  /*        pup = 0.5 * erfc(-zup./sqrt(2)); */
  /*    else */
  /*       return; */
  /*       %error(message('stats:logncdf:BadCovarianceSymPos')); */
  /*    end */
  /*  */
  /* end */
  return p;
}

/*
 * Arguments    : const emxArray_real_T *x
 *                const emxArray_real_T *y
 *                emxArray_real_T *z
 * Return Type  : void
 */
static void rdivide(const emxArray_real_T *x, const emxArray_real_T *y,
                    emxArray_real_T *z)
{
  int i1;
  int loop_ub;
  i1 = z->size[0] * z->size[1];
  z->size[0] = 1;
  z->size[1] = x->size[1];
  emxEnsureCapacity((emxArray__common *)z, i1, (int)sizeof(double));
  loop_ub = x->size[0] * x->size[1];
  for (i1 = 0; i1 < loop_ub; i1++) {
    z->data[i1] = x->data[i1] / y->data[i1];
  }
}

/*
 * Arguments    : double u0
 *                double u1
 * Return Type  : double
 */
static double rt_powd_snf(double u0, double u1)
{
  double y;
  double d1;
  double d2;
  if (rtIsNaN(u0) || rtIsNaN(u1)) {
    y = rtNaN;
  } else {
    d1 = fabs(u0);
    d2 = fabs(u1);
    if (rtIsInf(u1)) {
      if (d1 == 1.0) {
        y = rtNaN;
      } else if (d1 > 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 (d2 == 0.0) {
      y = 1.0;
    } else if (d2 == 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;
}

/*
 * Arguments    : const emxArray_real_T *x
 * Return Type  : double
 */
static double sum(const emxArray_real_T *x)
{
  double y;
  int k;
  if (x->size[1] == 0) {
    y = 0.0;
  } else {
    y = x->data[0];
    for (k = 2; k <= x->size[1]; k++) {
      y += x->data[k - 1];
    }
  }

  return y;
}

/*
 * Arguments    : int numDimensions
 *                int *size
 * Return Type  : emxArray_real_T *
 */
emxArray_real_T *emxCreateND_real_T(int numDimensions, int *size)
{
  emxArray_real_T *emx;
  int numEl;
  int i;
  b_emxInit_real_T(&emx, numDimensions);
  numEl = 1;
  for (i = 0; i < numDimensions; i++) {
    numEl *= size[i];
    emx->size[i] = size[i];
  }

  emx->data = (double *)calloc((unsigned int)numEl, sizeof(double));
  emx->numDimensions = numDimensions;
  emx->allocatedSize = numEl;
  return emx;
}

/*
 * Arguments    : double *data
 *                int numDimensions
 *                int *size
 * Return Type  : emxArray_real_T *
 */
emxArray_real_T *emxCreateWrapperND_real_T(double *data, int numDimensions, int *
  size)
{
  emxArray_real_T *emx;
  int numEl;
  int i;
  b_emxInit_real_T(&emx, numDimensions);
  numEl = 1;
  for (i = 0; i < numDimensions; i++) {
    numEl *= size[i];
    emx->size[i] = size[i];
  }

  emx->data = data;
  emx->numDimensions = numDimensions;
  emx->allocatedSize = numEl;
  emx->canFreeData = false;
  return emx;
}

/*
 * Arguments    : double *data
 *                int rows
 *                int cols
 * Return Type  : emxArray_real_T *
 */
emxArray_real_T *emxCreateWrapper_real_T(double *data, int rows, int cols)
{
  emxArray_real_T *emx;
  int size[2];
  int numEl;
  int i;
  size[0] = rows;
  size[1] = cols;
  b_emxInit_real_T(&emx, 2);
  numEl = 1;
  for (i = 0; i < 2; i++) {
    numEl *= size[i];
    emx->size[i] = size[i];
  }

  emx->data = data;
  emx->numDimensions = 2;
  emx->allocatedSize = numEl;
  emx->canFreeData = false;
  return emx;
}

/*
 * Arguments    : int rows
 *                int cols
 * Return Type  : emxArray_real_T *
 */
emxArray_real_T *emxCreate_real_T(int rows, int cols)
{
  emxArray_real_T *emx;
  int size[2];
  int numEl;
  int i;
  size[0] = rows;
  size[1] = cols;
  b_emxInit_real_T(&emx, 2);
  numEl = 1;
  for (i = 0; i < 2; i++) {
    numEl *= size[i];
    emx->size[i] = size[i];
  }

  emx->data = (double *)calloc((unsigned int)numEl, sizeof(double));
  emx->numDimensions = 2;
  emx->allocatedSize = numEl;
  return emx;
}

/*
 * Arguments    : emxArray_real_T *emxArray
 * Return Type  : void
 */
void emxDestroyArray_real_T(emxArray_real_T *emxArray)
{
  emxFree_real_T(&emxArray);
}

/*
 * Arguments    : emxArray_real_T **pEmxArray
 *                int numDimensions
 * Return Type  : void
 */
void emxInitArray_real_T(emxArray_real_T **pEmxArray, int numDimensions)
{
  b_emxInit_real_T(pEmxArray, numDimensions);
}

/*
 * Dolo�itev ciljne intenzitete in pripadajo�ega potresnega scenarija
 * Arguments    : double lambdaCt
 *                const emxArray_real_T *kpn_x
 *                const emxArray_real_T *kpn_y
 *                double betaC
 *                double TRagR
 *                double *b_imc
 *                double *imr
 * Return Type  : void
 */
void imc(double lambdaCt, const emxArray_real_T *kpn_x, const emxArray_real_T
         *kpn_y, double betaC, double TRagR, double *b_imc, double *imr)
{
  int ixstart;
  int n;
  double mtmp;
  int ix;
  boolean_T exitg5;
  double anew;
  double apnd;
  double ndbl;
  double cdiff;
  emxArray_real_T *im1;
  int i0;
  int k;
  double kd;
  boolean_T exitg4;
  double absb;
  emxArray_real_T *im2;
  emxArray_real_T *PC;
  emxArray_real_T *r0;
  emxArray_real_T *m;
  emxArray_real_T *Hint;
  emxArray_real_T *imCtM;
  emxArray_real_T *lambdaCa;
  emxArray_real_T *err;
  unsigned int j;
  double imCheck;
  emxArray_real_T *x;
  emxArray_real_T *b_lambdaCa;
  emxArray_real_T *b_imCtM;
  emxArray_real_T *c_imCtM;
  emxArray_real_T *d_imCtM;
  emxArray_real_T *e_imCtM;
  emxArray_real_T *f_imCtM;
  emxArray_real_T *g_imCtM;
  boolean_T exitg3;
  double imMax;
  double PC1;
  double d0;
  boolean_T guard1 = false;
  double b_x;
  double h_imCtM;
  double i_imCtM;
  double j_imCtM;
  emxArray_boolean_T *c_x;
  double y;
  int ii_data[1];
  boolean_T exitg2;
  int indagR__data[1];
  int indagR;
  int indagR_[2];
  double varargin_1[2];
  int b_indagR_[2];
  double varargin_2[2];
  double imr_;
  int32_T exitg1;
  int b_indagR[2];
  double r;

  /* intenzitete */
  /* frekvence */
  /*  Izra�un mediane ciljne intenzitete z numeri�no integracijo ena�be tveganja */
  ixstart = 1;
  n = kpn_x->size[0];
  mtmp = kpn_x->data[0];
  if (kpn_x->size[0] > 1) {
    if (rtIsNaN(kpn_x->data[0])) {
      ix = 2;
      exitg5 = false;
      while ((!exitg5) && (ix <= n)) {
        ixstart = ix;
        if (!rtIsNaN(kpn_x->data[ix - 1])) {
          mtmp = kpn_x->data[ix - 1];
          exitg5 = true;
        } else {
          ix++;
        }
      }
    }

    if (ixstart < kpn_x->size[0]) {
      while (ixstart + 1 <= n) {
        if (kpn_x->data[ixstart] > mtmp) {
          mtmp = kpn_x->data[ixstart];
        }

        ixstart++;
      }
    }
  }

  if (rtIsNaN(mtmp)) {
    n = 0;
    anew = rtNaN;
    apnd = mtmp;
  } else if (mtmp < 0.001) {
    n = -1;
    anew = 0.001;
    apnd = mtmp;
  } else if (rtIsInf(mtmp)) {
    n = 0;
    anew = rtNaN;
    apnd = mtmp;
  } else {
    anew = 0.001;
    ndbl = floor((mtmp - 0.001) / 0.001 + 0.5);
    apnd = 0.001 + ndbl * 0.001;
    cdiff = apnd - mtmp;
    if (fabs(cdiff) < 4.4408920985006262E-16 * mtmp) {
      ndbl++;
      apnd = mtmp;
    } else if (cdiff > 0.0) {
      apnd = 0.001 + (ndbl - 1.0) * 0.001;
    } else {
      ndbl++;
    }

    if (ndbl >= 0.0) {
      n = (int)ndbl - 1;
    } else {
      n = -1;
    }
  }

  emxInit_real_T(&im1, 2);
  i0 = im1->size[0] * im1->size[1];
  im1->size[0] = 1;
  im1->size[1] = n + 1;
  emxEnsureCapacity((emxArray__common *)im1, i0, (int)sizeof(double));
  if (n + 1 > 0) {
    im1->data[0] = anew;
    if (n + 1 > 1) {
      im1->data[n] = apnd;
      ixstart = n / 2;
      for (k = 1; k < ixstart; k++) {
        kd = (double)k * 0.001;
        im1->data[k] = anew + kd;
        im1->data[n - k] = apnd - kd;
      }

      if (ixstart << 1 == n) {
        im1->data[ixstart] = (anew + apnd) / 2.0;
      } else {
        kd = (double)ixstart * 0.001;
        im1->data[ixstart] = anew + kd;
        im1->data[ixstart + 1] = apnd - kd;
      }
    }
  }

  ixstart = 1;
  n = kpn_x->size[0];
  mtmp = kpn_x->data[0];
  if (kpn_x->size[0] > 1) {
    if (rtIsNaN(kpn_x->data[0])) {
      ix = 2;
      exitg4 = false;
      while ((!exitg4) && (ix <= n)) {
        ixstart = ix;
        if (!rtIsNaN(kpn_x->data[ix - 1])) {
          mtmp = kpn_x->data[ix - 1];
          exitg4 = true;
        } else {
          ix++;
        }
      }
    }

    if (ixstart < kpn_x->size[0]) {
      while (ixstart + 1 <= n) {
        if (kpn_x->data[ixstart] > mtmp) {
          mtmp = kpn_x->data[ixstart];
        }

        ixstart++;
      }
    }
  }

  if (rtIsNaN(mtmp + 0.001)) {
    n = 0;
    anew = rtNaN;
    apnd = mtmp + 0.001;
  } else if (mtmp + 0.001 < 0.0005) {
    n = -1;
    anew = 0.0005;
    apnd = mtmp + 0.001;
  } else if (rtIsInf(mtmp + 0.001)) {
    n = 0;
    anew = rtNaN;
    apnd = mtmp + 0.001;
  } else {
    anew = 0.0005;
    ndbl = floor(((mtmp + 0.001) - 0.0005) / 0.001 + 0.5);
    apnd = 0.0005 + ndbl * 0.001;
    cdiff = apnd - (mtmp + 0.001);
    absb = fabs(mtmp + 0.001);
    if ((0.0005 >= absb) || rtIsNaN(absb)) {
      absb = 0.0005;
    }

    if (fabs(cdiff) < 4.4408920985006262E-16 * absb) {
      ndbl++;
      apnd = mtmp + 0.001;
    } else if (cdiff > 0.0) {
      apnd = 0.0005 + (ndbl - 1.0) * 0.001;
    } else {
      ndbl++;
    }

    if (ndbl >= 0.0) {
      n = (int)ndbl - 1;
    } else {
      n = -1;
    }
  }

  emxInit_real_T(&im2, 2);
  i0 = im2->size[0] * im2->size[1];
  im2->size[0] = 1;
  im2->size[1] = n + 1;
  emxEnsureCapacity((emxArray__common *)im2, i0, (int)sizeof(double));
  if (n + 1 > 0) {
    im2->data[0] = anew;
    if (n + 1 > 1) {
      im2->data[n] = apnd;
      ixstart = n / 2;
      for (k = 1; k < ixstart; k++) {
        kd = (double)k * 0.001;
        im2->data[k] = anew + kd;
        im2->data[n - k] = apnd - kd;
      }

      if (ixstart << 1 == n) {
        im2->data[ixstart] = (anew + apnd) / 2.0;
      } else {
        kd = (double)ixstart * 0.001;
        im2->data[ixstart] = anew + kd;
        im2->data[ixstart + 1] = apnd - kd;
      }
    }
  }

  b_emxInit_real_T(&PC, 1);

  /* Hint = exp(pchip_(log(x),log(y),log(im2))); */
  /* Hint = exp(pchip_(log(x),log(y),log(im2'))); */
  i0 = PC->size[0];
  PC->size[0] = kpn_x->size[0];
  emxEnsureCapacity((emxArray__common *)PC, i0, (int)sizeof(double));
  ixstart = kpn_x->size[0];
  for (i0 = 0; i0 < ixstart; i0++) {
    PC->data[i0] = kpn_x->data[i0];
  }

  b_emxInit_real_T(&r0, 1);
  b_log(PC);
  i0 = r0->size[0];
  r0->size[0] = kpn_y->size[0];
  emxEnsureCapacity((emxArray__common *)r0, i0, (int)sizeof(double));
  ixstart = kpn_y->size[0];
  for (i0 = 0; i0 < ixstart; i0++) {
    r0->data[i0] = kpn_y->data[i0];
  }

  emxInit_real_T(&m, 2);
  b_log(r0);
  i0 = m->size[0] * m->size[1];
  m->size[0] = 1;
  m->size[1] = im2->size[1];
  emxEnsureCapacity((emxArray__common *)m, i0, (int)sizeof(double));
  ixstart = im2->size[0] * im2->size[1];
  for (i0 = 0; i0 < ixstart; i0++) {
    m->data[i0] = im2->data[i0];
  }

  emxInit_real_T(&Hint, 2);
  emxInit_real_T(&imCtM, 2);
  emxInit_real_T(&lambdaCa, 2);
  emxInit_real_T(&err, 2);
  c_log(m);
  interp1(PC, r0, m, Hint);
  b_exp(Hint);
  diff(Hint, m);
  diff(im2, Hint);
  rdivide(m, Hint, im2);
  b_abs(im2, m);
  j = 1U;
  i0 = imCtM->size[0] * imCtM->size[1];
  imCtM->size[0] = 1;
  imCtM->size[1] = 1;
  emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
  imCtM->data[0] = 1.0;
  imCheck = 1.0;
  i0 = lambdaCa->size[0] * lambdaCa->size[1];
  lambdaCa->size[0] = 1;
  lambdaCa->size[1] = 0;
  emxEnsureCapacity((emxArray__common *)lambdaCa, i0, (int)sizeof(double));
  i0 = err->size[0] * err->size[1];
  err->size[0] = 1;
  err->size[1] = 0;
  emxEnsureCapacity((emxArray__common *)err, i0, (int)sizeof(double));

  /* % */
  emxInit_real_T(&x, 2);
  emxInit_real_T(&b_lambdaCa, 2);
  emxInit_real_T(&b_imCtM, 2);
  emxInit_real_T(&c_imCtM, 2);
  emxInit_real_T(&d_imCtM, 2);
  emxInit_real_T(&e_imCtM, 2);
  emxInit_real_T(&f_imCtM, 2);
  emxInit_real_T(&g_imCtM, 2);
  while (imCheck > 0.001) {
    ixstart = 1;
    n = kpn_x->size[0];
    mtmp = kpn_x->data[0];
    if (kpn_x->size[0] > 1) {
      if (rtIsNaN(kpn_x->data[0])) {
        ix = 2;
        exitg3 = false;
        while ((!exitg3) && (ix <= n)) {
          ixstart = ix;
          if (!rtIsNaN(kpn_x->data[ix - 1])) {
            mtmp = kpn_x->data[ix - 1];
            exitg3 = true;
          } else {
            ix++;
          }
        }
      }

      if (ixstart < kpn_x->size[0]) {
        while (ixstart + 1 <= n) {
          if (kpn_x->data[ixstart] > mtmp) {
            mtmp = kpn_x->data[ixstart];
          }

          ixstart++;
        }
      }
    }

    imMax = mtmp;
    PC1 = logncdf_(mtmp, log(imCtM->data[(int)j - 1]), betaC);
    while (PC1 < 0.99) {
      imMax++;
      PC1 = logncdf_(imMax, log(imCtM->data[(int)j - 1]), betaC);
      if (rtIsNaN(imMax)) {
        n = 0;
        anew = rtNaN;
        apnd = imMax;
      } else if (imMax < 0.001) {
        n = -1;
        anew = 0.001;
        apnd = imMax;
      } else if (rtIsInf(imMax)) {
        n = 0;
        anew = rtNaN;
        apnd = imMax;
      } else {
        anew = 0.001;
        ndbl = floor((imMax - 0.001) / 0.001 + 0.5);
        apnd = 0.001 + ndbl * 0.001;
        cdiff = apnd - imMax;
        if (fabs(cdiff) < 4.4408920985006262E-16 * imMax) {
          ndbl++;
          apnd = imMax;
        } else if (cdiff > 0.0) {
          apnd = 0.001 + (ndbl - 1.0) * 0.001;
        } else {
          ndbl++;
        }

        if (ndbl >= 0.0) {
          n = (int)ndbl - 1;
        } else {
          n = -1;
        }
      }

      i0 = im1->size[0] * im1->size[1];
      im1->size[0] = 1;
      im1->size[1] = n + 1;
      emxEnsureCapacity((emxArray__common *)im1, i0, (int)sizeof(double));
      if (n + 1 > 0) {
        im1->data[0] = anew;
        if (n + 1 > 1) {
          im1->data[n] = apnd;
          ixstart = n / 2;
          for (k = 1; k < ixstart; k++) {
            kd = (double)k * 0.001;
            im1->data[k] = anew + kd;
            im1->data[n - k] = apnd - kd;
          }

          if (ixstart << 1 == n) {
            im1->data[ixstart] = (anew + apnd) / 2.0;
          } else {
            kd = (double)ixstart * 0.001;
            im1->data[ixstart] = anew + kd;
            im1->data[ixstart + 1] = apnd - kd;
          }
        }
      }

      if (rtIsNaN(imMax + 0.001)) {
        n = 0;
        anew = rtNaN;
        apnd = imMax + 0.001;
      } else if (imMax + 0.001 < 0.0005) {
        n = -1;
        anew = 0.0005;
        apnd = imMax + 0.001;
      } else if (rtIsInf(imMax + 0.001)) {
        n = 0;
        anew = rtNaN;
        apnd = imMax + 0.001;
      } else {
        anew = 0.0005;
        ndbl = floor(((imMax + 0.001) - 0.0005) / 0.001 + 0.5);
        apnd = 0.0005 + ndbl * 0.001;
        cdiff = apnd - (imMax + 0.001);
        absb = fabs(imMax + 0.001);
        if ((0.0005 >= absb) || rtIsNaN(absb)) {
          absb = 0.0005;
        }

        if (fabs(cdiff) < 4.4408920985006262E-16 * absb) {
          ndbl++;
          apnd = imMax + 0.001;
        } else if (cdiff > 0.0) {
          apnd = 0.0005 + (ndbl - 1.0) * 0.001;
        } else {
          ndbl++;
        }

        if (ndbl >= 0.0) {
          n = (int)ndbl - 1;
        } else {
          n = -1;
        }
      }

      i0 = im2->size[0] * im2->size[1];
      im2->size[0] = 1;
      im2->size[1] = n + 1;
      emxEnsureCapacity((emxArray__common *)im2, i0, (int)sizeof(double));
      if (n + 1 > 0) {
        im2->data[0] = anew;
        if (n + 1 > 1) {
          im2->data[n] = apnd;
          ixstart = n / 2;
          for (k = 1; k < ixstart; k++) {
            kd = (double)k * 0.001;
            im2->data[k] = anew + kd;
            im2->data[n - k] = apnd - kd;
          }

          if (ixstart << 1 == n) {
            im2->data[ixstart] = (anew + apnd) / 2.0;
          } else {
            kd = (double)ixstart * 0.001;
            im2->data[ixstart] = anew + kd;
            im2->data[ixstart + 1] = apnd - kd;
          }
        }
      }

      /* Hint   = exp(pchip_(log(x),log(y),log(im2))); */
      /* m      = abs(diff(Hint)./diff(im2)); */
      /*          im2a=max(x)+dim/2:dim:imMax+dim; */
      i0 = PC->size[0];
      PC->size[0] = kpn_x->size[0];
      emxEnsureCapacity((emxArray__common *)PC, i0, (int)sizeof(double));
      ixstart = kpn_x->size[0];
      for (i0 = 0; i0 < ixstart; i0++) {
        PC->data[i0] = kpn_x->data[i0];
      }

      b_log(PC);
      i0 = r0->size[0];
      r0->size[0] = kpn_y->size[0];
      emxEnsureCapacity((emxArray__common *)r0, i0, (int)sizeof(double));
      ixstart = kpn_y->size[0];
      for (i0 = 0; i0 < ixstart; i0++) {
        r0->data[i0] = kpn_y->data[i0];
      }

      b_log(r0);
      i0 = m->size[0] * m->size[1];
      m->size[0] = 1;
      m->size[1] = im2->size[1];
      emxEnsureCapacity((emxArray__common *)m, i0, (int)sizeof(double));
      ixstart = im2->size[0] * im2->size[1];
      for (i0 = 0; i0 < ixstart; i0++) {
        m->data[i0] = im2->data[i0];
      }

      c_log(m);
      interp1(PC, r0, m, x);
      i0 = Hint->size[0] * Hint->size[1];
      Hint->size[0] = 1;
      Hint->size[1] = x->size[1];
      emxEnsureCapacity((emxArray__common *)Hint, i0, (int)sizeof(double));
      ixstart = x->size[0] * x->size[1];
      for (i0 = 0; i0 < ixstart; i0++) {
        Hint->data[i0] = x->data[i0];
      }

      for (k = 0; k < x->size[1]; k++) {
        Hint->data[k] = exp(Hint->data[k]);
      }

      diff(Hint, m);
      diff(im2, Hint);
      rdivide(m, Hint, im2);
      b_abs(im2, m);

      /*          m2=[m ma]; */
    }

    /*      if extrap== */
    /*          m2=m; */
    /*      end */
    i0 = Hint->size[0] * Hint->size[1];
    Hint->size[0] = 1;
    Hint->size[1] = 0;
    emxEnsureCapacity((emxArray__common *)Hint, i0, (int)sizeof(double));
    i0 = PC->size[0];
    PC->size[0] = im1->size[1];
    emxEnsureCapacity((emxArray__common *)PC, i0, (int)sizeof(double));
    ixstart = im1->size[1];
    for (i0 = 0; i0 < ixstart; i0++) {
      PC->data[i0] = 0.0;
    }

    for (ixstart = 0; ixstart < im1->size[1]; ixstart++) {
      PC->data[ixstart] = logncdf_(im1->data[ixstart], log(imCtM->data[(int)j -
        1]), betaC);
      ix = Hint->size[1];
      i0 = Hint->size[0] * Hint->size[1];
      Hint->size[1] = ix + 1;
      emxEnsureCapacity((emxArray__common *)Hint, i0, (int)sizeof(double));
      Hint->data[ix] = PC->data[ixstart] * m->data[ixstart] * 0.001;

      /* lambdaCa1(i)=PC(i)*m2(i)*dim; */
    }

    d0 = sum(Hint);
    i0 = b_lambdaCa->size[0] * b_lambdaCa->size[1];
    b_lambdaCa->size[0] = 1;
    b_lambdaCa->size[1] = lambdaCa->size[1] + 1;
    emxEnsureCapacity((emxArray__common *)b_lambdaCa, i0, (int)sizeof(double));
    ixstart = lambdaCa->size[1];
    for (i0 = 0; i0 < ixstart; i0++) {
      b_lambdaCa->data[b_lambdaCa->size[0] * i0] = lambdaCa->data[lambdaCa->
        size[0] * i0];
    }

    b_lambdaCa->data[b_lambdaCa->size[0] * lambdaCa->size[1]] = d0;
    i0 = lambdaCa->size[0] * lambdaCa->size[1];
    lambdaCa->size[0] = 1;
    lambdaCa->size[1] = b_lambdaCa->size[1];
    emxEnsureCapacity((emxArray__common *)lambdaCa, i0, (int)sizeof(double));
    ixstart = b_lambdaCa->size[1];
    for (i0 = 0; i0 < ixstart; i0++) {
      lambdaCa->data[lambdaCa->size[0] * i0] = b_lambdaCa->data[b_lambdaCa->
        size[0] * i0];
    }

    /* lambdaCa(j)=sum(lambdaCa1); */
    /*  lambdaCa=PC*m'*dim; */
    ix = err->size[1];
    i0 = err->size[0] * err->size[1];
    err->size[1] = ix + 1;
    emxEnsureCapacity((emxArray__common *)err, i0, (int)sizeof(double));
    err->data[ix] = lambdaCa->data[(int)j - 1] - lambdaCt;
    guard1 = false;
    if ((int)j == 1) {
      guard1 = true;
    } else {
      i0 = x->size[0] * x->size[1];
      x->size[0] = 1;
      x->size[1] = err->size[1];
      emxEnsureCapacity((emxArray__common *)x, i0, (int)sizeof(double));
      ixstart = err->size[0] * err->size[1];
      for (i0 = 0; i0 < ixstart; i0++) {
        x->data[i0] = err->data[i0];
      }

      for (k = 0; k < err->size[1]; k++) {
        if (x->data[k] < 0.0) {
          b_x = -1.0;
        } else if (x->data[k] > 0.0) {
          b_x = 1.0;
        } else if (x->data[k] == 0.0) {
          b_x = 0.0;
        } else {
          b_x = x->data[k];
        }

        x->data[k] = b_x;
      }

      b_x = sum(x);
      if (fabs(b_x) == j) {
        guard1 = true;
      } else if (err->data[(int)j - 1] > 0.0) {
        h_imCtM = imCtM->data[(int)j - 1];
        i_imCtM = imCtM->data[(int)j - 1];
        j_imCtM = imCtM->data[(int)j - 2];
        i0 = d_imCtM->size[0] * d_imCtM->size[1];
        d_imCtM->size[0] = 1;
        d_imCtM->size[1] = imCtM->size[1] + 1;
        emxEnsureCapacity((emxArray__common *)d_imCtM, i0, (int)sizeof(double));
        ixstart = imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          d_imCtM->data[d_imCtM->size[0] * i0] = imCtM->data[imCtM->size[0] * i0];
        }

        d_imCtM->data[d_imCtM->size[0] * imCtM->size[1]] = h_imCtM + fabs
          (i_imCtM - j_imCtM) / 2.0;
        i0 = imCtM->size[0] * imCtM->size[1];
        imCtM->size[0] = 1;
        imCtM->size[1] = d_imCtM->size[1];
        emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
        ixstart = d_imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          imCtM->data[imCtM->size[0] * i0] = d_imCtM->data[d_imCtM->size[0] * i0];
        }
      } else if (err->data[(int)j - 1] < 0.0) {
        h_imCtM = imCtM->data[(int)j - 1];
        i_imCtM = imCtM->data[(int)j - 1];
        j_imCtM = imCtM->data[(int)j - 2];
        i0 = c_imCtM->size[0] * c_imCtM->size[1];
        c_imCtM->size[0] = 1;
        c_imCtM->size[1] = imCtM->size[1] + 1;
        emxEnsureCapacity((emxArray__common *)c_imCtM, i0, (int)sizeof(double));
        ixstart = imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          c_imCtM->data[c_imCtM->size[0] * i0] = imCtM->data[imCtM->size[0] * i0];
        }

        c_imCtM->data[c_imCtM->size[0] * imCtM->size[1]] = h_imCtM - fabs
          (i_imCtM - j_imCtM) / 2.0;
        i0 = imCtM->size[0] * imCtM->size[1];
        imCtM->size[0] = 1;
        imCtM->size[1] = c_imCtM->size[1];
        emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
        ixstart = c_imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          imCtM->data[imCtM->size[0] * i0] = c_imCtM->data[c_imCtM->size[0] * i0];
        }
      } else {
        h_imCtM = imCtM->data[(int)j - 1];
        i0 = b_imCtM->size[0] * b_imCtM->size[1];
        b_imCtM->size[0] = 1;
        b_imCtM->size[1] = imCtM->size[1] + 1;
        emxEnsureCapacity((emxArray__common *)b_imCtM, i0, (int)sizeof(double));
        ixstart = imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          b_imCtM->data[b_imCtM->size[0] * i0] = imCtM->data[imCtM->size[0] * i0];
        }

        b_imCtM->data[b_imCtM->size[0] * imCtM->size[1]] = h_imCtM;
        i0 = imCtM->size[0] * imCtM->size[1];
        imCtM->size[0] = 1;
        imCtM->size[1] = b_imCtM->size[1];
        emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
        ixstart = b_imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          imCtM->data[imCtM->size[0] * i0] = b_imCtM->data[b_imCtM->size[0] * i0];
        }
      }
    }

    if (guard1) {
      if (err->data[(int)j - 1] > 0.0) {
        h_imCtM = imCtM->data[(int)j - 1];
        i0 = g_imCtM->size[0] * g_imCtM->size[1];
        g_imCtM->size[0] = 1;
        g_imCtM->size[1] = imCtM->size[1] + 1;
        emxEnsureCapacity((emxArray__common *)g_imCtM, i0, (int)sizeof(double));
        ixstart = imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          g_imCtM->data[g_imCtM->size[0] * i0] = imCtM->data[imCtM->size[0] * i0];
        }

        g_imCtM->data[g_imCtM->size[0] * imCtM->size[1]] = h_imCtM * 2.0;
        i0 = imCtM->size[0] * imCtM->size[1];
        imCtM->size[0] = 1;
        imCtM->size[1] = g_imCtM->size[1];
        emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
        ixstart = g_imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          imCtM->data[imCtM->size[0] * i0] = g_imCtM->data[g_imCtM->size[0] * i0];
        }
      } else if (err->data[(int)j - 1] < 0.0) {
        h_imCtM = imCtM->data[(int)j - 1];
        i0 = f_imCtM->size[0] * f_imCtM->size[1];
        f_imCtM->size[0] = 1;
        f_imCtM->size[1] = imCtM->size[1] + 1;
        emxEnsureCapacity((emxArray__common *)f_imCtM, i0, (int)sizeof(double));
        ixstart = imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          f_imCtM->data[f_imCtM->size[0] * i0] = imCtM->data[imCtM->size[0] * i0];
        }

        f_imCtM->data[f_imCtM->size[0] * imCtM->size[1]] = h_imCtM / 2.0;
        i0 = imCtM->size[0] * imCtM->size[1];
        imCtM->size[0] = 1;
        imCtM->size[1] = f_imCtM->size[1];
        emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
        ixstart = f_imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          imCtM->data[imCtM->size[0] * i0] = f_imCtM->data[f_imCtM->size[0] * i0];
        }
      } else {
        h_imCtM = imCtM->data[(int)j - 1];
        i0 = e_imCtM->size[0] * e_imCtM->size[1];
        e_imCtM->size[0] = 1;
        e_imCtM->size[1] = imCtM->size[1] + 1;
        emxEnsureCapacity((emxArray__common *)e_imCtM, i0, (int)sizeof(double));
        ixstart = imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          e_imCtM->data[e_imCtM->size[0] * i0] = imCtM->data[imCtM->size[0] * i0];
        }

        e_imCtM->data[e_imCtM->size[0] * imCtM->size[1]] = h_imCtM;
        i0 = imCtM->size[0] * imCtM->size[1];
        imCtM->size[0] = 1;
        imCtM->size[1] = e_imCtM->size[1];
        emxEnsureCapacity((emxArray__common *)imCtM, i0, (int)sizeof(double));
        ixstart = e_imCtM->size[1];
        for (i0 = 0; i0 < ixstart; i0++) {
          imCtM->data[imCtM->size[0] * i0] = e_imCtM->data[e_imCtM->size[0] * i0];
        }
      }
    }

    imCheck = fabs(imCtM->data[(int)((double)j + 1.0) - 1] - imCtM->data[(int)j
                   - 1]);
    j++;
  }

  emxFree_real_T(&g_imCtM);
  emxFree_real_T(&f_imCtM);
  emxFree_real_T(&e_imCtM);
  emxFree_real_T(&d_imCtM);
  emxFree_real_T(&c_imCtM);
  emxFree_real_T(&b_imCtM);
  emxFree_real_T(&b_lambdaCa);
  emxFree_real_T(&r0);
  emxFree_real_T(&x);
  emxFree_real_T(&PC);
  emxFree_real_T(&err);
  emxFree_real_T(&lambdaCa);
  emxFree_real_T(&m);
  emxFree_real_T(&Hint);
  emxFree_real_T(&im2);
  emxFree_real_T(&im1);
  emxInit_boolean_T(&c_x, 1);
  *b_imc = imCtM->data[imCtM->size[1] - 1];
  y = 1.0 / TRagR;
  i0 = c_x->size[0];
  c_x->size[0] = kpn_y->size[0];
  emxEnsureCapacity((emxArray__common *)c_x, i0, (int)sizeof(boolean_T));
  ixstart = kpn_y->size[0];
  emxFree_real_T(&imCtM);
  for (i0 = 0; i0 < ixstart; i0++) {
    c_x->data[i0] = (kpn_y->data[i0] < y);
  }

  k = (1 <= c_x->size[0]);
  ixstart = 0;
  ix = 1;
  exitg2 = false;
  while ((!exitg2) && (ix <= c_x->size[0])) {
    if (c_x->data[ix - 1]) {
      ixstart = 1;
      ii_data[0] = ix;
      exitg2 = true;
    } else {
      ix++;
    }
  }

  emxFree_boolean_T(&c_x);
  if (k == 1) {
    if (ixstart == 0) {
      k = 0;
    }
  } else {
    k = !(1 > ixstart);
  }

  i0 = 0;
  while (i0 <= k - 1) {
    indagR__data[0] = ii_data[0];
    i0 = 1;
  }

  indagR = indagR__data[0];
  if (indagR__data[0] > 1) {
    y = 1.0 / TRagR;
    indagR_[0] = (int)((double)indagR__data[0] - 1.0);
    indagR_[1] = indagR__data[0];
    for (i0 = 0; i0 < 2; i0++) {
      varargin_1[i0] = kpn_y->data[indagR_[i0] - 1];
    }

    b_indagR_[0] = (int)((double)indagR__data[0] - 1.0);
    b_indagR_[1] = indagR__data[0];
    for (i0 = 0; i0 < 2; i0++) {
      varargin_2[i0] = kpn_x->data[b_indagR_[i0] - 1];
    }

    imr_ = rtNaN;
    k = 1;
    do {
      exitg1 = 0;
      if (k < 3) {
        b_indagR[0] = (int)((double)indagR - 1.0);
        b_indagR[1] = indagR;
        if (rtIsNaN(kpn_y->data[b_indagR[k - 1] - 1])) {
          exitg1 = 1;
        } else {
          k++;
        }
      } else {
        if (kpn_y->data[indagR - 1] < kpn_y->data[(int)((double)indagR - 1.0) -
            1]) {
          varargin_1[0] = kpn_y->data[indagR - 1];
          varargin_1[1] = kpn_y->data[(int)((double)indagR__data[0] - 1.0) - 1];
          varargin_2[0] = kpn_x->data[indagR - 1];
          varargin_2[1] = kpn_x->data[(int)((double)indagR__data[0] - 1.0) - 1];
        }

        if (rtIsNaN(y) || (y > varargin_1[1]) || (y < varargin_1[0])) {
        } else {
          r = (y - varargin_1[0]) / (varargin_1[1] - varargin_1[0]);
          if (r == 0.0) {
            imr_ = varargin_2[0];
          } else if (r == 1.0) {
            imr_ = varargin_2[1];
          } else if (varargin_2[0] == varargin_2[1]) {
            imr_ = varargin_2[0];
          } else {
            imr_ = (1.0 - r) * varargin_2[0] + r * varargin_2[1];
          }
        }

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

    *imr = imr_;
  } else {
    *imr = kpn_x->data[0];
  }

  /*  */
  /*  % */
  /*  kvanAll = zeros(nAcc,1); */
  /*  for i=1:nAcc */
  /*     kvanAll(i) = (i-0.5)/nAcc; */
  /*  end */
  /*  [~,index] = min(abs(kvanAll-0.1587)); */
  /*   */
  /*  % Izracun ciljne intenzitete za izbran kvantil (0.1625 za skupino s 40 akcelerogrami) */
  /*  ss = size(imCtM,2); */
  /*  if ss ~= 1 */
  /*      imCt = imCtM(ss)*exp(norminv_(kvanAll(index),0,1)*betaC); */
  /*  end */
  /*   */
  /*  % Dolo�itev epsilona za pogojni spekter */
  /*  %[median_sa,sigma]=SP_1996(T1,M,R,1,0,4); */
  /*  %epsilon=(log(imCt16)-log(median_sa))/sigma; */
  /*   */
  /*  % Dolo�itev povratne dobe, ki pripada ciljni intenziteti */
  /*  %RP=1/exp(pchip_(log(x),log(y),log(imCt16))); */
  /*   */
  /*  %Himct */
  /*  Himct  = exp(pchip_(log(x),log(y),log(imCt))); */
  /*  PCimct = logncdf_(imCt,log(imCtM(ss)),betaC); */
  /*   */
  /*  saUser         = kpn_x'; */
  /*  hUser          = kpn_y'; */
  /*   */
  /*  isExtrap = extrap; */
  /*  if(isExtrap == 1) */
  /*  	% Shorten interpolated HF to max 200 points (or greater than user defined function) */
  /*  	while(size(Hint,2) > 200 && size(Hint,2) > size(hUser,2)) */
  /*  		Hint_tmp = Hint(2:2:end-1); */
  /*  		im2_tmp  = im2(2:2:end-1); */
  /*  		Hint     = [Hint(1,1) Hint_tmp Hint(1,end)]; */
  /*  		im2      = [im2(1,1) im2_tmp im2(1,end)]; */
  /*      end */
  /*   */
  /*  	% Modify interpolated Hazard function to include(Sact) */
  /*  	if ~any(Hint == Himct) */
  /*  	    hInterpolated  = [Hint(Hint > Himct) Himct Hint(Hint < Himct)]; */
  /*  	    saInterpolated = [im2(im2 < imCt) imCt im2(im2 > imCt)]; */
  /*  	else */
  /*  	    hInterpolated  = Hint; */
  /*  	    saInterpolated = im2; */
  /*      end */
  /*  	% Calc PC */
  /*  	PC_ = zeros(1,size(saInterpolated,2)); */
  /*  	for i = 1:size(saInterpolated,2) */
  /*  		PC_(1,i) = logncdf_(saInterpolated(1,i),log(imCtM(ss)),betaC); */
  /*  	end */
  /*  else */
  /*  	saInterpolated = []; */
  /*  	hInterpolated  = []; */
  /*  	if ~any(hUser == Himct) */
  /*  	    hUser  = [hUser(hUser > Himct) Himct hUser(hUser < Himct)]; */
  /*  	    saUser = [saUser(saUser < imCt) imCt saUser(saUser > imCt)]; */
  /*  	end */
  /*   */
  /*  	% Calc PC */
  /*  	PC_ = zeros(1,size(saUser,2)); */
  /*  	for i = 1:size(saUser,2) */
  /*  		PC_(1,i) = logncdf_(saUser(1,i),log(imCtM(ss)),betaC); */
  /*  	end */
  /*  end */
}

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

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

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