frenel.c

Numerical Recipes in C的源代码

#include <math.h>
#include "complex.h"
#define EPS 6.0e-8
#define MAXIT 100
#define FPMIN 1.0e-30
#define XMIN 1.5
#define PI 3.1415927
#define PIBY2 (PI/2.0)
#define TRUE 1
#define ONE Complex(1.0,0.0)

void frenel(float x, float *s, float *c)
{
	void nrerror(char error_text[]);
	int k,n,odd;
	float a,ax,fact,pix2,sign,sum,sumc,sums,term,test;
	fcomplex b,cc,d,h,del,cs;

	ax=fabs(x);
	if (ax < sqrt(FPMIN)) {
		*s=0.0;
		*c=ax;
	} else if (ax <= XMIN) {
		sum=sums=0.0;
		sumc=ax;
		sign=1.0;
		fact=PIBY2*ax*ax;
		odd=TRUE;
		term=ax;
		n=3;
		for (k=1;k<=MAXIT;k++) {
			term *= fact/k;
			sum += sign*term/n;
			test=fabs(sum)*EPS;
			if (odd) {
				sign = -sign;
				sums=sum;
				sum=sumc;
			} else {
				sumc=sum;
				sum=sums;
			}
			if (term < test) break;
			odd=!odd;
			n += 2;
		}
		if (k > MAXIT) nrerror("series failed in frenel");
		*s=sums;
		*c=sumc;
	} else {
		pix2=PI*ax*ax;
		b=Complex(1.0,-pix2);
		cc=Complex(1.0/FPMIN,0.0);
		d=h=Cdiv(ONE,b);
		n = -1;
		for (k=2;k<=MAXIT;k++) {
			n += 2;
			a = -n*(n+1);
			b=Cadd(b,Complex(4.0,0.0));
			d=Cdiv(ONE,Cadd(RCmul(a,d),b));
			cc=Cadd(b,Cdiv(Complex(a,0.0),cc));
			del=Cmul(cc,d);
			h=Cmul(h,del);
			if (fabs(del.r-1.0)+fabs(del.i) < EPS) break;
		}
		if (k > MAXIT) nrerror("cf failed in frenel");
		h=Cmul(Complex(ax,-ax),h);
		cs=Cmul(Complex(0.5,0.5),
			Csub(ONE,Cmul(Complex(cos(0.5*pix2),sin(0.5*pix2)),h)));
		*c=cs.r;
		*s=cs.i;
	}
	if (x < 0.0) {
		*c = -(*c);
		*s = -(*s);
	}
}
#undef EPS
#undef MAXIT
#undef FPMIN
#undef XMIN
#undef PI
#undef PIBY2
#undef TRUE
#undef ONE