44 lines
1.0 KiB
C
44 lines
1.0 KiB
C
|
|
#include <math.h>
|
|
#define PI 3.1415927
|
|
#define THIRD (1.0/3.0)
|
|
#define TWOTHR (2.0*THIRD)
|
|
#define ONOVRT 0.57735027
|
|
|
|
void airy(float x, float *ai, float *bi, float *aip, float *bip)
|
|
{
|
|
void bessik(float x, float xnu, float *ri, float *rk, float *rip,
|
|
float *rkp);
|
|
void bessjy(float x, float xnu, float *rj, float *ry, float *rjp,
|
|
float *ryp);
|
|
float absx,ri,rip,rj,rjp,rk,rkp,rootx,ry,ryp,z;
|
|
|
|
absx=fabs(x);
|
|
rootx=sqrt(absx);
|
|
z=TWOTHR*absx*rootx;
|
|
if (x > 0.0) {
|
|
bessik(z,THIRD,&ri,&rk,&rip,&rkp);
|
|
*ai=rootx*ONOVRT*rk/PI;
|
|
*bi=rootx*(rk/PI+2.0*ONOVRT*ri);
|
|
bessik(z,TWOTHR,&ri,&rk,&rip,&rkp);
|
|
*aip = -x*ONOVRT*rk/PI;
|
|
*bip=x*(rk/PI+2.0*ONOVRT*ri);
|
|
} else if (x < 0.0) {
|
|
bessjy(z,THIRD,&rj,&ry,&rjp,&ryp);
|
|
*ai=0.5*rootx*(rj-ONOVRT*ry);
|
|
*bi = -0.5*rootx*(ry+ONOVRT*rj);
|
|
bessjy(z,TWOTHR,&rj,&ry,&rjp,&ryp);
|
|
*aip=0.5*absx*(ONOVRT*ry+rj);
|
|
*bip=0.5*absx*(ONOVRT*rj-ry);
|
|
} else {
|
|
*ai=0.35502805;
|
|
*bi=(*ai)/ONOVRT;
|
|
*aip = -0.25881940;
|
|
*bip = -(*aip)/ONOVRT;
|
|
}
|
|
}
|
|
#undef PI
|
|
#undef THIRD
|
|
#undef TWOTHR
|
|
#undef ONOVRT
|