complex<typename boost::units::subtract_typeof_helper<X,Y>::type>
operator-(const complex<X>& x,const complex<Y>& y)
{
    typedef typename boost::units::subtract_typeof_helper<X,Y>::type    type;
    
    return complex<type>(x.real()-y.real(),x.imag()-y.imag());
}

void EMLTracker::integrate(complex<double> in)
{
    localReplica.tick();

    // First bring the signal level of the input to the same as the
    // local replicas
    in *= baseGain;

    complex<double> carrier = localReplica.getCarrier();
    complex<double> code = localReplica.getCode() ? plusOne : minusOne;

    // mix in the carrier local replica
    complex<double> m0 = in * conj(carrier);

    // and sum them up.. (yea, the conj of the codes should be a NoOp)
    code = conj(code);
    early.process(m0, code);
    prompt.process(m0, code);
    late.process(m0, code);

    // Update our sums for normalizing things...
    complex<double> lr = conj(carrier) * code;
    inSumSq += in.real()*in.real() + in.imag()*in.imag();
    lrSumSq += lr.real()*lr.real() + lr.imag()*lr.imag();
}

complex complex::operator *(complex x)
{
    complex y;
    y.Re=Re*x.take_Re()-Im*x.take_Im();
    y.Im=Re*x.take_Im()+Im*x.take_Re();
    return y;
}

/** Implement Exponential Integral function, E1(z), based on formulaes in
 *Abramowitz and Stegun (A&S)
 *  In fact this implementation returns exp(z)*E1(z) where z is a complex number
 *
 *  @param z :: input
 *  @return exp(z)*E1(z)
 */
complex<double> exponentialIntegral(const complex<double> &z) {
  double z_abs = abs(z);

  if (z_abs == 0.0) {
    // Basically function is infinite in along the real axis for this case
    return complex<double>(std::numeric_limits<double>::max(), 0.0);
  } else if (z_abs < 10.0) // 10.0 is a guess based on formula 5.1.55 (no idea
                           // how good it really is)
  {
    // use formula 5.1.11 in A&S. rewrite last term in 5.1.11 as
    // x*sum_n=0 (-1)^n x^n / (n+1)*(n+1)! and then calculate the
    // terms in the sum recursively
    complex<double> z1(1.0, 0.0);
    complex<double> z2(1.0, 0.0);
    for (int i = 1; i <= 100; i++) // is max of 100 here best?
    {
      z2 = -static_cast<double>(i) * (z2 * z) / ((i + 1.0) * (i + 1.0));
      z1 += z2;
      if (abs(z2) < 1.0E-10 * abs(z1))
        break; // i.e. if break loop if little change to term added
    }
    return exp(z) * (-log(z) - M_EULER + z * z1);
  } else {
    // use formula 5.1.22 in A&S. See discussion page 231
    complex<double> z1(0.0, 0.0);
    for (int i = 20; i >= 1; i--) // not sure if 20 is always sensible?
      z1 = static_cast<double>(i) / (1.0 + static_cast<double>(i) / (z + z1));
    complex<double> retVal = 1.0 / (z + z1);
    if (z.real() <= 0.0 && z.imag() == 0.0)
      retVal -= exp(z) * complex<double>(0.0, 1.0) * M_PI; // see formula 5.1.7
    return retVal;
  }
}

complex<double> cos(complex<double> x) { 
  double cosh(double x);
  double cos(double x);
  double r = cosh(x.imag())*cos(x.real());
  double i = sinh(x.imag())*sin(x.real());
  return(complex<double>(r,i));
}

complex<Float> Comms::GlobalSum(complex<Float> z)
{

  if (!initialized) { Initialize(); }

#ifdef USE_MPI
  Float sendbuf[2];
  Float recvbuf[2];

  sendbuf[0] = z.real();
  sendbuf[1] = z.imag();

  MPI_Barrier(MPI_COMM_WORLD);

#ifdef USE_SINGLE
  MPI_Reduce( sendbuf, recvbuf, 2, MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD );
#endif
#ifdef USE_DOUBLE
  MPI_Reduce( sendbuf, recvbuf, 2, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
#endif
#ifdef USE_LONG_DOUBLE
  MPI_Reduce( sendbuf, recvbuf, 2, MPI_LONG_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
#endif

  return complex<Float>(recvbuf[0], recvbuf[1]);
#else
  return z;
#endif

}

void write_mandel(std::string name, int x, int y)
{
  GIF_image gif(x, y);

  for(int i=0; i<256;i++)
  {
	  gif.setColor(i, i, i, (char)255);
  }

  for (int i = 0; i < y; i++) 
  { 
    for (int j = 0; j < x; j++) 
    {
      complex<double> c(c1.real()+j*(c2-c1).real()/X, c1.imag()+i*(c2-c1).imag()/Y);
      complex<double> z = 0;
      int iter = 0;
      
      while (abs(z) < 2 && iter < MAX_ITER) 
      {
		z = z*z+c;
		iter++;
      }
      gif.setPixel(j, i, ((MAX_ITER-iter)*255)/40);
    }
  }
  
  gif.writeToFile(name.c_str());
}

complex operator- (const complex& c1, const complex& c2)
{
    double res_re = c1.get_re() - c2.get_re();
    double res_im = c1.get_im() - c2.get_im();
    complex res{res_re, res_im};

    return res;
}

complex<Y> 
pow(const complex<Y>& x,const Y& y)
{
    std::complex<Y> tmp(x.real(),x.imag());
    
    tmp = std::pow(tmp,y);
    
    return complex<Y>(tmp.real(),tmp.imag());
}

 static type value(const complex<quantity<Unit,Y> >& x)  
 { 
     const complex<value_type>   tmp =
         root<static_rational<N,D> >(complex<Y>(x.real().value(),
                                                x.imag().value()));
     
     return type(quantity_type::from_value(tmp.real()),
                quantity_type::from_value(tmp.imag()));
 }

complex<double> complexAdd(complex<double> x1, complex<double> x2)
{
    int x1_real = int(x1.real());
    int x2_real = int(x2.real());
    int x1_imag = int(x1.imag());
    int x2_imag = int(x2.imag());
    int r1 = int(bitset<16>(x1_real+x2_real).to_ulong());
    int r2 = int(bitset<16>(x1_imag+x2_imag).to_ulong());
    return complex<double>(double(r1), double(r2));
}

QTextStream& operator<<(QTextStream& ts, KCValue value)
{
    ts << value.type();
    switch (value.type()) {
    case KCValue::Empty:   break;

    case KCValue::Boolean:
        ts << ": ";
        if (value.asBoolean()) ts << "TRUE";
        else ts << "FALSE"; break;

    case KCValue::Integer:
        ts << ": " << value.asInteger(); break;

    case KCValue::Float:
        ts << ": " << (double) numToDouble(value.asFloat()); break;

    case KCValue::Complex: {
        const complex<KCNumber> complex(value.asComplex());
        ts << ": " << (double) numToDouble(complex.real());
        if (complex.imag() >= 0.0)
            ts << '+';
        ts << (double) numToDouble(complex.imag()) << 'i';
        break;
    }

    case KCValue::String:
        ts << ": " << value.asString(); break;

    case KCValue::Array: {
        ts << ": {" << value.asString();
        const int cols = value.columns();
        const int rows = value.rows();
        for (int row = 0; row < rows; ++row) {
            for (int col = 0; col < cols; ++col) {
                ts << value.element(col, row);
                if (col < cols - 1)
                    ts << ';';
            }
            if (row < rows - 1)
                ts << '|';
        }
        ts << '}';
        break;
    }

    case KCValue::Error:
        ts << '(' << value.errorMessage() << ')'; break;

    default: break;
    }
    return ts;
}

std::string fmtvalue(std::true_type, const complex<T>& x)
{
    std::string restr = as_string(fmt<'g', number_width, number_precision>(x.real()));
    if (restr.size() > number_width)
        restr = as_string(fmt<'g', number_width, number_precision_short>(x.real()));

    std::string imstr = as_string(fmt<'g', -1, number_precision>(std::abs(x.imag())));
    if (imstr.size() > number_width)
        imstr = as_string(fmt<'g', -1, number_precision_short>(std::abs(x.imag())));

    return restr + (x.imag() < T(0) ? "-" : "+") + padleft(number_width, imstr + "j");
}

/* Complex root forward deflation.
 */
static int complexdeflation( const register int n, double a[], const complex<double> z )
   {
   register int i;
   double r, u;

   r = -2.0 * z.real();
   u = z.real() * z.real() + z.imag() * z.imag();
   a[ 1 ] -= r * a[ 0 ];
   for( i = 2; i < n - 1; i++ )
      a[ i ] = a[ i ] - r * a[ i - 1 ] - u * a[ i - 2 ];

   return n - 2;
}

complex<long double> sum_to_all(complex<long double> in) {
  complex<long double> out = in;
#ifdef HAVE_MPI
  if (MPI_LONG_DOUBLE == MPI_DATATYPE_NULL) {
    complex<double> dout;
    dout = sum_to_all(complex<double>(double(in.real()), double(in.imag())));
    out = complex<long double>(dout.real(), dout.imag());
  }
  else
    MPI_Allreduce(&in,&out,2,MPI_LONG_DOUBLE,MPI_SUM,mycomm);
#endif
  return out;
}

void transmissioncom::sendToTransmissionSolver(complex value)
{
  uint8_t data[sizeof(double)*2];
  double realPart=value.Re();
  double imagineryPart=value.Im();
  memcpy(&data,&realPart,sizeof(double));
  memcpy(&data[sizeof(double)],&imagineryPart,sizeof(double));
  
  OBJECT *hdr = OBJECTHDR(this);
  Message *msg=new Message(hdr->name,this->connectedBus,gl_globalclock,(uint8_t*)data,sizeof(double)*2);
  msg->setDelayThroughComm(false);
  myinterface->send(msg);
}

complex<double> complexMultiply(complex<double> x1, complex<double> x2)
{
    int x1_real = int(x1.real());
    int x2_real = int(x2.real());
    int x1_imag = int(x1.imag());
    int x2_imag = int(x2.imag());
    int r1 = int(bitset<16>(x1_real*x2_real).to_ulong());
    int r2 = int(bitset<16>(x1_imag*x2_imag).to_ulong());
    int r3 = int(bitset<16>(x1_real*x2_imag).to_ulong());
    int r4 = int(bitset<16>(x1_imag*x2_real).to_ulong());
    int m_real = int(bitset<16>(r1-r2).to_ulong());
    int m_imag = int(bitset<16>(r3+r4).to_ulong());
    return complex<double>(double(m_real), double(m_imag));
}

void Biquad::setZeroPolePairs(const complex<double> &zero, const complex<double> &pole)
{
    double b0 = 1;
    double b1 = -2 * zero.real();

    double zeroMag = abs(zero);
    double b2 = zeroMag * zeroMag;

    double a1 = -2 * pole.real();

    double poleMag = abs(pole);
    double a2 = poleMag * poleMag;
    setNormalizedCoefficients(b0, b1, b2, 1, a1, a2);
}

  void put(index_type i, complex<T> val)
  {
    assert(this->format_ != no_user_format);

    if (this->format_ == array_format)
      this->u_.data_[i] = val;
    else if (this->format_ == interleaved_format)
    {
      this->u_.split_.real_[2*i+0] = val.real();
      this->u_.split_.real_[2*i+1] = val.imag();
    }
    else // if (format_ == split_format)
    {
      this->u_.split_.real_[i] = val.real();
      this->u_.split_.imag_[i] = val.imag();
    }
  }

double microturbine::determine_frequency(complex power_out){
	//double Kx = 3 * Rinternal / 3.14159;
	//used linear approximation which does better job than textbook equation...
	//assumes power in W
	double f = power_out.Mag() * 1.5 + 55000;\
	f = f/60;
	return f;
}