#include <cmath>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#define CREAL_T
typedef float real32_T;
typedef double real64_T;
typedef struct {
  real32_T re;
  real32_T im;
} creal32_T;

typedef struct {
  real64_T re;
  real64_T im;
} creal_T; 
// Function Declarations
static double rt_powd_snf(double u0, double u1);

void CalculateAppxIndex(double r, double epsilon_s, double epsilon_h, double//Noncompliant
  lambda, double f, creal_T *eps_eff, creal_T *mu_eff)
{
  double n_h;
  double m;
  double a;
  double alpha_e_re;
  double alpha_e_im;
  double alpha_m_re;
  double alpha_m_im;
  int l;
  double br;

  //  sphere radius is 1 micron;
  //  the refractive index of inclusion
  n_h = std::sqrt(epsilon_h);

  //  the refractive index of host
  m = std::sqrt(epsilon_s) / n_h;
  n_h = 6.2831853071795862 * n_h * r / lambda;

  //  size parameter;
  a = (n_h + m) + 2.0;
  m = n_h - m;

  //  polarizability
  n_h = 6.2831853071795862 * rt_powd_snf(r, 3.0) / rt_powd_snf(n_h, 3.0);
  alpha_e_re = 0.0;
  alpha_e_im = 0.0;
  alpha_m_re = 0.0;
  alpha_m_im = 0.0;
  for (l = 0; l < 5; l++) {
    alpha_e_re += 0.0 * ((2.0 * (1.0 + (double)l) + 1.0) * a);
    alpha_e_im += (2.0 * (1.0 + (double)l) + 1.0) * a;
    alpha_m_re += 0.0 * ((2.0 * (1.0 + (double)l) + 1.0) * (m + 2.0));
    alpha_m_im += (2.0 * (1.0 + (double)l) + 1.0) * (m + 2.0);

    //  alpha = alpha + 1i* [(2*l+1)*(an(l) + bn(l))];
  }

  alpha_e_re *= n_h;
  alpha_e_im *= n_h;
  alpha_m_re *= n_h;
  alpha_m_im *= n_h;

  // alpha = aa*alpha;
  n_h = f / (4.1887902047863905 * rt_powd_snf(r, 3.0));
  alpha_e_re *= n_h;
  alpha_e_im *= n_h;
  alpha_m_re *= n_h;
  alpha_m_im *= n_h;
  if (alpha_e_im == 0.0) {
    m = alpha_e_re / 3.0;
    n_h = 0.0;
  } else if (alpha_e_re == 0.0) {
    m = 0.0;
    n_h = alpha_e_im / 3.0;
  } else {
    m = alpha_e_re / 3.0;
    n_h = alpha_e_im / 3.0;
  }

  br = 1.0 - m;
  m = 0.0 - n_h;
  if (m == 0.0) {
    if (alpha_e_im == 0.0) {
      m = alpha_e_re / br;
      alpha_e_im = 0.0;
    } else if (alpha_e_re == 0.0) {
      m = 0.0;
      alpha_e_im /= br;
    } else {
      m = alpha_e_re / br;
      alpha_e_im /= br;
    }
  } else {
    n_h = std::abs(m);
    if (br > n_h) {
      a = m / br;
      n_h = br + a * m;
      m = (alpha_e_re + a * alpha_e_im) / n_h;
      alpha_e_im = (alpha_e_im - a * alpha_e_re) / n_h;
    } else if (n_h == br) {
      if (br > 0.0) {
        a = 0.5;
      } else {
        a = -0.5;
      }

      if (m > 0.0) {
        n_h = 0.5;
      } else {
        n_h = -0.5;
      }

      m = alpha_e_re * a + alpha_e_im * n_h;
      alpha_e_im = alpha_e_im * a - alpha_e_re * n_h;
    } else {
      a = br / m;
      n_h = m + a * br;
      m = (a * alpha_e_re + alpha_e_im) / n_h;
      alpha_e_im = (a * alpha_e_im - alpha_e_re) / n_h;
    }
  }

  eps_eff->re = epsilon_h * (1.0 + m);
  eps_eff->im = epsilon_h * alpha_e_im;
  if (alpha_m_im == 0.0) {
    m = alpha_m_re / 3.0;
    n_h = 0.0;
  } else if (alpha_m_re == 0.0) {
    m = 0.0;
    n_h = alpha_m_im / 3.0;
  } else {
    m = alpha_m_re / 3.0;
    n_h = alpha_m_im / 3.0;
  }

  br = 1.0 - m;
  m = 0.0 - n_h;
  if (m == 0.0) {
    if (alpha_m_im == 0.0) {
      m = alpha_m_re / br;
      alpha_m_im = 0.0;
    } else if (alpha_m_re == 0.0) {
      m = 0.0;
      alpha_m_im /= br;
    } else {
      m = alpha_m_re / br;
      alpha_m_im /= br;
    }
  } else {
    n_h = std::abs(m);
    if (br > n_h) {
      a = m / br;
      n_h = br + a * m;
      m = (alpha_m_re + a * alpha_m_im) / n_h;
      alpha_m_im = (alpha_m_im - a * alpha_m_re) / n_h;
    } else if (n_h == br) {
      if (br > 0.0) {
        a = 0.5;
      } else {
        a = -0.5;
      }

      if (m > 0.0) {
        n_h = 0.5;
      } else {
        n_h = -0.5;
      }

      m = alpha_m_re * a + alpha_m_im * n_h;
      alpha_m_im = alpha_m_im * a - alpha_m_re * n_h;
    } else {
      a = br / m;
      n_h = m + a * br;
      m = (a * alpha_m_re + alpha_m_im) / n_h;
      alpha_m_im = (a * alpha_m_im - alpha_m_re) / n_h;
    }
  }

  mu_eff->re = 1.0 + m;
  mu_eff->im = alpha_m_im;
}