C++ GSL_IMAG函数代码示例

gsl_complex gsl_complex_sqrt(gsl_complex a)
{				/* z=sqrt(a) */
    gsl_complex z;

    if (GSL_REAL(a) == 0.0 && GSL_IMAG(a) == 0.0) {
	GSL_SET_COMPLEX(&z, 0, 0);
    } else {
	double x = fabs(GSL_REAL(a));
	double y = fabs(GSL_IMAG(a));
	double w;

	if (x >= y) {
	    double t = y / x;
	    w = sqrt(x) * sqrt(0.5 * (1.0 + sqrt(1.0 + t * t)));
	} else {
	    double t = x / y;
	    w = sqrt(y) * sqrt(0.5 * (t + sqrt(1.0 + t * t)));
	}

	if (GSL_REAL(a) >= 0.0) {
	    double ai = GSL_IMAG(a);
	    GSL_SET_COMPLEX(&z, w, ai / (2.0 * w));
	} else {
	    double ai = GSL_IMAG(a);
	    double vi = (ai >= 0) ? w : -w;
	    GSL_SET_COMPLEX(&z, ai / (2.0 * vi), vi);
	}
    }

    return z;
}

int
gsl_sf_complex_psi_e(
  const double x,
  const double y,
  gsl_sf_result * result_re,
  gsl_sf_result * result_im
  )
{
  if(x >= 0.0)
  {
    gsl_complex z = gsl_complex_rect(x, y);
    return psi_complex_rhp(z, result_re, result_im);
  }
  else
  {
    /* reflection formula [Abramowitz+Stegun, 6.3.7] */
    gsl_complex z = gsl_complex_rect(x, y);
    gsl_complex omz = gsl_complex_rect(1.0 - x, -y);
    gsl_complex zpi = gsl_complex_mul_real(z, M_PI);
    gsl_complex cotzpi = gsl_complex_cot(zpi);
    int ret_val = psi_complex_rhp(omz, result_re, result_im);

    if(GSL_IS_REAL(GSL_REAL(cotzpi)) && GSL_IS_REAL(GSL_IMAG(cotzpi)))
    {
      result_re->val -= M_PI * GSL_REAL(cotzpi);
      result_im->val -= M_PI * GSL_IMAG(cotzpi);
      return ret_val;
    }
    else
    {
      GSL_ERROR("singularity", GSL_EDOM);
    }
  }
}

//
// we update shadowing and pathloss coefficients for all the links tx -> rx
//
void  MCPMPCoeffs::SpatialChannelUpdate() {

	for (int i=0; i<_M; i++) { // user i (tx)
		for (int ii=0;ii<=i;ii++) { // user ii (rx)

			// we are considering the spatial channel tx-rx, so we get the sostheta(j) vector
			// and we compute the shadowing coefficient for the position of rx (rx_lon,rx_lat)

			gsl_complex txPos =  gsl_vector_complex_get(geoPositions,i); // lat,lon
			gsl_complex rxPos =  gsl_vector_complex_get(geoPositions,ii+_M); // lat,lon

			double txPosLat = GSL_REAL(txPos);
			double txPosLon = GSL_IMAG(txPos);
			double rxPosLat = GSL_REAL(rxPos);
			double rxPosLon = GSL_IMAG(rxPos);

			double deltaLat = rxPosLat-txPosLat; // ok
			double deltaLon = rxPosLon-txPosLon; // be careful around 180E/W
			double meanLat = (rxPosLat+txPosLat)/2.0; // be careful around poles

			if (deltaLon > 180) // es Lon1 =-179 Lon2=179 D=2
				deltaLon -= 360;
			if (deltaLon < -180)
				deltaLon += 360;

			double x= deltaLon * 111111 * gsl_sf_cos(meanLat*M_PI_OVER_180); // cartesian x position around tx (in m)
			double y= deltaLat * 111111; // cartesian y position arond tx (in m)

			// shadowing i->ii
			double shadowdb=0;
			for (int j=0; j<SOSN; j++){
				double th = gsl_matrix_get(sostheta,i*_M+ii,j);
				double fxn = gsl_vector_get(sosfxn,j);
				double fyn = gsl_vector_get(sosfyn,j);
				shadowdb += sosc * gsl_sf_cos(2.0 * M_PI * ( fxn * x + fyn * y + th));
			} // for j

			//
			// PATHLOSS MODEL
			//
			double dist = sqrt(x*x+y*y); // metres

			double plossdb = mudisp::OkumuraHataCitydB(dist,1500);
			//double plossdb = mudisp::OkumuraHataMetrodB(dist,1500);


			//cout << "rx[" << i << "] x=" << x << " y=" << y << " ploss="
			//	<< -plossdb << " shadow=" << shadowdb << " dist=" << dist << endl;

			//
		    double cx = mudisp::dbtolin(-0.5 * plossdb + SOSsigma() * shadowdb );

		     gsl_matrix_set(pathLoss,i,ii,cx);
		     gsl_matrix_set(pathLoss,ii,i,cx); // it's symmetric !

		} // for ii
	} // for i
}

gsl_complex gsl_complex_add(gsl_complex a, gsl_complex b)
{				/* z=a+b */
    double ar = GSL_REAL(a), ai = GSL_IMAG(a);
    double br = GSL_REAL(b), bi = GSL_IMAG(b);

    gsl_complex z;
    GSL_SET_COMPLEX(&z, ar + br, ai + bi);
    return z;
}

gsl_complex gsl_complex_sub(gsl_complex a, gsl_complex b)
{				/* z=a-b */
    double ar = GSL_REAL(a), ai = GSL_IMAG(a);
    double br = GSL_REAL(b), bi = GSL_IMAG(b);

    gsl_complex z;
    GSL_SET_COMPLEX(&z, ar - br, ai - bi);
    return z;
}

gsl_complex gsl_complex_mul(gsl_complex a, gsl_complex b)
{				/* z=a*b */
    double ar = GSL_REAL(a), ai = GSL_IMAG(a);
    double br = GSL_REAL(b), bi = GSL_IMAG(b);

    gsl_complex z;
    GSL_SET_COMPLEX(&z, ar * br - ai * bi, ar * bi + ai * br);
    return z;
}

std::ostream& operator<<(std::ostream &os, const complex_spinor & c_spinor){

	int i;
	int nSitios = c_spinor._numSites;
	for(i=0;i<nSitios;i++){
	os << "[(" <<  GSL_REAL(c_spinor.complex_spinor_get(i,UP))<<" + i. "<<GSL_IMAG(c_spinor.complex_spinor_get(i,UP))<<")]"<<std::endl;
	os << "[(" <<  GSL_REAL(c_spinor.complex_spinor_get(i,DOWN))<<" + i. "<<GSL_IMAG(c_spinor.complex_spinor_get(i,DOWN))<<")]"<<std::endl;
	}
	return os; 
}

void
test_eigen_genherm_results (const gsl_matrix_complex * A, 
                            const gsl_matrix_complex * B,
                            const gsl_vector * eval, 
                            const gsl_matrix_complex * evec, 
                            size_t count,
                            const char * desc,
                            const char * desc2)
{
  const size_t N = A->size1;
  size_t i, j;

  gsl_vector_complex * x = gsl_vector_complex_alloc(N);
  gsl_vector_complex * y = gsl_vector_complex_alloc(N);

  /* check A v = lambda B v */
  for (i = 0; i < N; i++)
    {
      double ei = gsl_vector_get (eval, i);
      gsl_vector_complex_const_view vi =
        gsl_matrix_complex_const_column(evec, i);
      double norm = gsl_blas_dznrm2(&vi.vector);

      /* check that eigenvector is normalized */
      gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
                   "genherm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count,
                   desc, i, desc2);

      /* compute y = A z */
      gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, &vi.vector, GSL_COMPLEX_ZERO, y);

      /* compute x = B z */
      gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, B, &vi.vector, GSL_COMPLEX_ZERO, x);

      /* compute x = lambda B z */
      gsl_blas_zdscal(ei, x);

      /* now test if y = x */
      for (j = 0; j < N; j++)
        {
          gsl_complex xj = gsl_vector_complex_get (x, j);
          gsl_complex yj = gsl_vector_complex_get (y, j);

          gsl_test_rel(GSL_REAL(yj), GSL_REAL(xj), 1e9 * GSL_DBL_EPSILON, 
                       "genherm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
          gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e9 * GSL_DBL_EPSILON, 
                       "genherm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), imag, %s", N, count, desc, i, j, desc2);
        }
    }

  gsl_vector_complex_free(x);
  gsl_vector_complex_free(y);
}

void
FUNCTION (test, trap) (size_t stride, size_t N)
{
  TYPE (gsl_vector) * vc = FUNCTION (create, vector) (stride, N);

  BASE z = {{(ATOMIC)1.2, (ATOMIC)3.4}};
  BASE z1 = {{(ATOMIC)4.5, (ATOMIC)6.7}};

  size_t j = 0;

  status = 0;
  FUNCTION (gsl_vector, set) (vc, j - 1, z);
  gsl_test (!status,
            NAME (gsl_vector) "_set traps index below lower bound");

  status = 0;
  FUNCTION (gsl_vector, set) (vc, N + 1, z);
  gsl_test (!status,
            NAME (gsl_vector) "_set traps index above upper bound");

  status = 0;
  FUNCTION (gsl_vector, set) (vc, N, z);
  gsl_test (!status, NAME (gsl_vector) "_set traps index at upper bound");

  status = 0;
  z1 = FUNCTION (gsl_vector, get) (vc, j - 1);
  gsl_test (!status,
            NAME (gsl_vector) "_get traps index below lower bound");

  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_vector) "_get returns zero real below lower bound");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_vector) "_get returns zero imag below lower bound");

  status = 0;
  z1 = FUNCTION (gsl_vector, get) (vc, N + 1);
  gsl_test (!status,
            NAME (gsl_vector) "_get traps index above upper bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_vector) "_get returns zero real above upper bound");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_vector) "_get returns zero imag above upper bound");

  status = 0;
  z1 = FUNCTION (gsl_vector, get) (vc, N);
  gsl_test (!status, NAME (gsl_vector) "_get traps index at upper bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_vector) "_get returns zero real at upper bound");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_vector) "_get returns zero imag at upper bound");

  FUNCTION (gsl_vector, free) (vc);
}

int
main (void)
{
  double data[] = { -1.0, 1.0, -1.0, 1.0,
                    -8.0, 4.0, -2.0, 1.0,
                    27.0, 9.0, 3.0, 1.0,
                    64.0, 16.0, 4.0, 1.0 };

  gsl_matrix_view m 
    = gsl_matrix_view_array (data, 4, 4);

  gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
  gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);

  gsl_eigen_nonsymmv_workspace * w = 
    gsl_eigen_nonsymmv_alloc (4);
  
  gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);

  gsl_eigen_nonsymmv_free (w);

  gsl_eigen_nonsymmv_sort (eval, evec, 
                           GSL_EIGEN_SORT_ABS_DESC);
  
  {
    int i, j;

    for (i = 0; i < 4; i++)
      {
        gsl_complex eval_i 
           = gsl_vector_complex_get (eval, i);
        gsl_vector_complex_view evec_i 
           = gsl_matrix_complex_column (evec, i);

        printf ("eigenvalue = %g + %gi\n",
                GSL_REAL(eval_i), GSL_IMAG(eval_i));
        printf ("eigenvector = \n");
        for (j = 0; j < 4; ++j)
          {
            gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
            printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
          }
      }
  }

  gsl_vector_complex_free(eval);
  gsl_matrix_complex_free(evec);

  return 0;
}

gsl_complex
gsl_complex_pow_real (gsl_complex a, double b)
{                               /* z=a^b */
  gsl_complex z;

  if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0)
    {
      if (b == 0)
        {
          GSL_SET_COMPLEX (&z, 1, 0);
        }
      else
        {
          GSL_SET_COMPLEX (&z, 0, 0);
        }
    }
  else
    {
      double logr = gsl_complex_logabs (a);
      double theta = gsl_complex_arg (a);
      double rho = exp (logr * b);
      double beta = theta * b;
      GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));
    }

  return z;
}

bool VECTOR_ensure_not_complex(CVECTOR *_object)
{
	gsl_vector *v;
	int size = SIZE(THIS);
	int i;
	gsl_complex c;
	
	if (!COMPLEX(THIS))
		return FALSE;
	
	for (i = 0; i < size; i++)
	{
		c = gsl_vector_complex_get(CVEC(THIS), i);
		if (GSL_IMAG(c) != 0.0)
			return TRUE;
	}
	
	v = gsl_vector_alloc(size);
	
	for (i = 0; i < size; i++)
		gsl_vector_set(v, i, GSL_REAL(gsl_vector_complex_get(CVEC(THIS), i)));
	
	gsl_vector_complex_free(CVEC(THIS));
	THIS->vector = v;
	THIS->complex = FALSE;
	return FALSE;
}

static int
matlab_dlmwrite_complex(FILE * fp, const gsl_matrix_complex * A)
{
    const size_t M = A->size1;
    const size_t N = A->size2;
    size_t i, j;

    for (i = 0; i < M; ++i)
    {
        for (j = 0; j < N; ++j)
        {
            gsl_complex z = gsl_matrix_complex_get(A, i, j);
            double zr = GSL_REAL(z);
            double zi = GSL_IMAG(z);

            fprintf(fp, "%.12e%s%.12ei%s",
                    zr,
                    (zi < 0.0) ? "" : "+",
                    zi,
                    (j < N - 1) ? "," : "\n");
        }
    }

    return 0;
}

gsl_complex
gsl_complex_arccosh (gsl_complex a)
{				/* z = arccosh(a) */
  gsl_complex z = gsl_complex_arccos (a);
  z = gsl_complex_mul_imag (z, GSL_IMAG(z) > 0 ? -1.0 : 1.0);
  return z;
}

gsl_complex
gsl_complex_tan (gsl_complex a)
{				/* z = tan(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;

  if (fabs (I) < 1)
    {
      double D = pow (cos (R), 2.0) + pow (sinh (I), 2.0);

      GSL_SET_COMPLEX (&z, 0.5 * sin (2 * R) / D, 0.5 * sinh (2 * I) / D);
    }
  else
    {
      double u = exp (-I);
      double C = 2 * u / (1 - pow (u, 2.0));
      double D = 1 + pow (cos (R), 2.0) * pow (C, 2.0);

      double S = pow (C, 2.0);
      double T = 1.0 / tanh (I);

      GSL_SET_COMPLEX (&z, 0.5 * sin (2 * R) * S / D, T / D);
    }

  return z;
}