自学教程：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;}

intgsl_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; }

voidtest_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);}

voidFUNCTION (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);}

intmain (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_complexgsl_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 intmatlab_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_complexgsl_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_complexgsl_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;}

gsl_complexgsl_complex_add_real (gsl_complex a, double x){				/* z=a+x */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_REAL (a) + x, GSL_IMAG (a));  return z;}

double gsl_complex_abs2(gsl_complex z){				/* return |z|^2 */    double x = GSL_REAL(z);    double y = GSL_IMAG(z);    return (x * x + y * y);}

gsl_complexgsl_complex_div_imag (gsl_complex a, double y){				/* z=a/(iy) */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_IMAG (a) / y,  - GSL_REAL (a) / y);  return z;}

gsl_complexgsl_complex_div_real (gsl_complex a, double x){				/* z=a/x */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_REAL (a) / x, GSL_IMAG (a) / x);  return z;}

gsl_complexgsl_complex_mul_imag (gsl_complex a, double y){				/* z=a*iy */  gsl_complex z;  GSL_SET_COMPLEX (&z, -y * GSL_IMAG (a), y * GSL_REAL (a));  return z;}

voidgsl_complex_cosh (complex_t const *a, complex_t *res){                               /* z = cosh(a) */        gnm_float R = GSL_REAL (a), I = GSL_IMAG (a);	complex_init (res, cosh (R) * gnm_cos (I), gnm_sinh (R) * gnm_sin (I));}

gsl_complexgsl_complex_sub_imag (gsl_complex a, double y){				/* z=a-iy */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_REAL (a), GSL_IMAG (a) - y);  return z;}

gsl_complexgsl_complex_sub_real (gsl_complex a, double x){				/* z=a-x */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_REAL (a) - x, GSL_IMAG (a));  return z;}

gsl_complexgsl_complex_add_imag (gsl_complex a, double y){				/* z=a+iy */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_REAL (a), GSL_IMAG (a) + y);  return z;}

gsl_complexgsl_complex_conjugate (gsl_complex a){				/* z=conj(a) */  gsl_complex z;  GSL_SET_COMPLEX (&z, GSL_REAL (a), -GSL_IMAG (a));  return z;}

static VALUE rb_gsl_sf_lngamma_complex_e(int argc, VALUE *argv, VALUE obj){  gsl_sf_result *lnr, *arg;  gsl_complex *z;  double re, im;  VALUE vlnr, varg;  int status;  switch (argc) {  case 1:    CHECK_COMPLEX(argv[0]);    Data_Get_Struct(argv[0], gsl_complex, z);    re = GSL_REAL(*z);    im = GSL_IMAG(*z);    break;  case 2:    Need_Float(argv[0]); Need_Float(argv[1]);    re = NUM2DBL(argv[0]);    im = NUM2DBL(argv[1]);  default:    rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);  }  vlnr = Data_Make_Struct(cgsl_sf_result, gsl_sf_result, 0, free, lnr);  varg = Data_Make_Struct(cgsl_sf_result, gsl_sf_result, 0, free, arg);  status = gsl_sf_lngamma_complex_e(re, im, lnr, arg);  return rb_ary_new3(3, vlnr, varg, INT2FIX(status));}

gsl_complexgsl_complex_negative (gsl_complex a){				/* z=-a */  gsl_complex z;  GSL_SET_COMPLEX (&z, -GSL_REAL (a), -GSL_IMAG (a));  return z;}

gsl_complexgsl_complex_mul_real (gsl_complex a, double x){				/* z=a*x */  gsl_complex z;  GSL_SET_COMPLEX (&z, x * GSL_REAL (a), x * GSL_IMAG (a));  return z;}

voidgsl_complex_inverse (complex_t const *a, complex_t *res){                               /* z=1/a */        gnm_float s = 1.0 / complex_mod (a);	complex_init (res, (GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) * s);}