1 #ifdef HAVE_STDLIB_H
2 #include <stdlib.h>
3 #endif
4 #include <math.h>
5
6 #include "yagi.h"
7
8 #define EPS 6.0e-8
9 #define EULER 0.57721566
10 #define MAXIT 100
11 #define PIBY2 1.5707963
12 #define FPMIN 1.0e-30
13 #define TMIN 2.0
14 #define TRUE 1
15 #define ONE Complex(1.0,0.0)
16
cisi(double x,double * ci,double * si)17 void cisi(double x, double *ci, double *si)
18 {
19 int i,k,odd;
20 double a,err,fact,sign,sum,sumc,sums,t,term;
21 fcomplex h,b,c,d,del;
22
23 t=fabs(x);
24 if (t == 0.0) {
25 *si=0.0;
26 *ci = -1.0/FPMIN;
27 return;
28 }
29 if (t > TMIN) {
30 b=Complex(1.0,t);
31 c=Complex(1.0/FPMIN,0.0);
32 d=h=Cdiv(ONE,b);
33 for (i=2;i<=MAXIT;i++) {
34 a = -(i-1)*(i-1);
35 b=Cadd(b,Complex(2.0,0.0));
36 d=Cdiv(ONE,Cadd(RCmul(a,d),b));
37 c=Cadd(b,Cdiv(Complex(a,0.0),c));
38 del=Cmul(c,d);
39 h=Cmul(h,del);
40 if (fabs(del.r-1.0)+fabs(del.i) < EPS) break;
41 }
42 if (i > MAXIT) nrerror("cf failed in cisi_hack.c");
43 h=Cmul(Complex(cos(t),-sin(t)),h);
44 *ci = -h.r;
45 *si=PIBY2+h.i;
46 } else {
47 if (t < sqrt(FPMIN)) {
48 sumc=0.0;
49 sums=t;
50 } else {
51 sum=sums=sumc=0.0;
52 sign=fact=1.0;
53 odd=TRUE;
54 for (k=1;k<=MAXIT;k++) {
55 fact *= t/k;
56 term=fact/k;
57 sum += sign*term;
58 err=term/fabs(sum);
59 if (odd) {
60 sign = -sign;
61 sums=sum;
62 sum=sumc;
63 } else {
64 sumc=sum;
65 sum=sums;
66 }
67 if (err < EPS) break;
68 odd=!odd;
69 }
70 if (k > MAXIT) nrerror("maxits exceeded in cisi");
71 }
72 *si=sums;
73 *ci=sumc+log(t)+EULER;
74 }
75 if (x < 0.0) *si = -(*si);
76 }
77 #undef EPS
78 #undef EULER
79 #undef MAXIT
80 #undef PIBY2
81 #undef FPMIN
82 #undef TMIN
83 #undef TRUE
84 #undef ONE
85