complexf.c

该程序是用vc开发的对动态数组进行管理的DLL

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/************************************************************************
COMPLEXF  - Subroutines to perform operations on fcomplex numbers.
		This set of functions complement the one in fcomplex.c
		of the CWP library

cipow		raise a fcomplex number to an integer power
crpow		raise a fcomplex number to a real power
rcpow		raise a real number to a fcomplex power
ccpow		raise a fcomplex number to a fcomplex power
ccos		compute the fcomplex cosine of a fcomplex angle
csin		compute the fcomplex sine of a fcomplex angle
ccosh		compute the fcomplex hyperbolic cosine of a fcomplex angle
csinh		compute the fcomplex hyperbolic sine of a fcomplex angle
cexp1		compute the fcomplex exponential of a fcomplex number
clog		compute the fcomplex logarithm of a fcomplex number
************************************************************************
Function Prototypes:
fcomplex cipow(fcomplex a, int p);
fcomplex crpow(fcomplex a, float p);
fcomplex rcpow(float a, fcomplex p);
fcomplex ccpow (fcomplex a, fcomplex p)
fcomplex ccos(fcomplex a);
fcomplex csin(fcomplex a);
fcomplex ccosh(fcomplex a);
fcomplex csinh(fcomplex a);
fcomplex cexp1(fcomplex a);
fcomplex clog(fcomplex a);
************************************************************************
Credits:
	Dave Hale, original version in C++
	Gabriel Alvarez, translation to C
***********************************************************************/
/**************** end self doc ********************************/

#include "cwp.h"

fcomplex cipow(fcomplex a, int p)
{
	fcomplex res;
	fcomplex b;

	if (p==0) {
		return cmplx(1.0,0.0);
	} else if (a.r==0.0 && a.i==0.0) {
		return cmplx(0.0,0.0);
	} else {
		res=cmplx(1.0,0.0);
		b=a;
		if (p<0) {
			p = -p;
			b = cinv(b);
		}
		for(;;) {
			if (p&1) res = cmul(res,b);
			if ((p>>=1)==0)
				return res;
			else
				b = cmul(b,b);
		}
	}
}			

fcomplex crpow(fcomplex a, float p)
{
	float ar,ai,amp,phs;

	if (p==0.0) return cmplx(1.0,0.0);
	if (a.r==0.0 && a.i==0.0) return cmplx(0.0,0.0);

	ar = a.r; ai = a.i;
	amp = (float)(exp(0.5*p*log(ar*ar+ai*ai)));
	phs = (float)(p*atan2(ai,ar));

	return cmplx((float)(amp*cos(phs)),(float)(amp*sin(phs)));	
}

fcomplex rcpow(float a, fcomplex p)
{
	float pr,pi,loga,amp,phs;

	if (p.r==0.0 && p.i==0.0) return cmplx(1.0,0.0);
	if (a==0.0) return cmplx(0.0,0.0);
	pr = p.r; pi = p.i;
	loga = (float)(0.5*log(a*a));
	amp = (float)exp(pr*loga);
	phs = pi*loga;

	return cmplx((float)(amp*cos(phs)),(float)(amp*sin(phs)));	
}

fcomplex ccpow (fcomplex a, fcomplex p)
{
	float ar,ai,pr,pi,loga,arga,amp,phs;

	if (p.r==0.0 && p.i==0.0) return cmplx(1.0,0.0);
	if (a.r==0.0 && a.i==0.0) return cmplx(0.0,0.0);

	pr = p.r; pi = p.i; ar = a.r; ai = a.i;
	loga = (float)(0.5*log(ar*ar+ai*ai));
	arga = (float)(atan2(ai,ar));
	amp = (float)(exp(pr*loga-pi*arga));
	phs = pr*arga+pi*loga;

	return cmplx((float)(amp*cos(phs)),(float)(amp*sin(phs)));
}

fcomplex ccos(fcomplex a)
{
	return cmplx((float)(cos(a.r)*cosh(a.i)),(float)(-sin(a.r)*sinh(a.i)));
}

fcomplex csin(fcomplex a)
{
	return cmplx((float)(sin(a.r)*cosh(a.i)),(float)(cos(a.r)*sinh(a.i)));
}

fcomplex ccosh(fcomplex a)
{
	return cmplx((float)(cos(a.i)*cosh(a.r)),(float)(sin(a.i)*sinh(a.r)));
}

fcomplex csinh(fcomplex a)
{
	return cmplx((float)(cos(a.i)*sinh(a.r)),(float)(sin(a.i)*cosh(a.r)));
}

fcomplex cexp1(fcomplex a)
{
	float r=(float)(exp(a.r));
	return cmplx((float)(r*cos(a.i)),(float)(r*sin(a.i)));
}

fcomplex clog(fcomplex a)
{
	float ar=a.r,ai=a.i,h=(float)(sqrt(ar*ar+ai*ai));
	return cmplx((float)log(h),(float)atan2(ai,ar));
}