Listing 2

/* programs for computations of kelvin functions and their
   derivatives

DEPENDENCIES:
   header file complex.h required

*/

#include "complex.h"
#include <stdio.h>
#define max(a,b) (((a)>(b))? (a): (b))
#define abs(x) ((x)? (x):-(x))

/* test driver*/
void main (argc,argv) int argc; char **argv;
{
int i, maxit=50;
double x,ber,bei,beip,berp,kei,ker,keip,kerp;
FILE *fileid;
fileid=fopen("PLOT.KE","w");
fprintf(fileid," 4 \n");
for(i=1;i<=maxit;i++)
   {x=i*.25;
   ke(x,&ker,&kei,&kerp,&keip);
   fprintf(fileid," %f %e %e %e %e\n",x,ker,kei,kerp,keip);
   };
fclose(fileid);
fileid=fopen("PLOT.BE","w");
fprintf(fileid," 4 \n");
for(i=1;i<=maxit;i++)
   {x=i*.25;
   be(x,&ber,&bei,&berp,&beip);
   fprintf(fileid," %f %e %e %e %e \n",x,ber,bei,berp,beip);
   };
fclose(fileid);
exit(0);
}

/* print a complex number*/
printc(x) struct complex *x;
{printf(" %e %e ",(x->x),(x->y));return;
}

/* complex exponential*/
cexp( x,ans) struct complex *x,*ans;
{
double y,exp(),sin(),cos();
struct complex c1,c2;
y = exp (x->x);
c2.x= cos (x->y);
c2.y= sin (x->y);
CTREAL(c1,c2,y);
CLET(*ans,c1);
return;
}

/* Kelvin functions and their derivatives
 Rational approximations from Abramowitz & Stegun
 section 9.9 pp.384-5*/
ke(x,ker,kei,kerp,keip)double x,*ker,*kei,*kerp,*keip;
{
double y,z,al,sqrt(),log(),ber,bei,berp,beip;
struct complex CC,cex,cdum,f,phi,theta;
if (x<=8.)
 {
 y=x*x/64.;
 z=y*y;
 be(x,&ber,&bei,&berp,&beip);
 al=-log(.5*x);

 *ker=al* ber+.7853981634* bei
 +(((((((-.00002458*z+.00309699)*z-.19636347)*z+5.65539121)*z
 -60.60977451)*z+171.36272133)*z-59.O5819744)*z-.57721566);
 *kei=al* bei-.7853981634* ber
 +(((((((.00029532*z-.02695875)*z+1.17509064)*z-21.30060904)*z
 +124.2356965)*z-142.91827687)*z+6.76454936)*y);

 *kerp=al* berp- ber/x+beip*.7853981634
 +x*((((((-.00001075*z+.00116137)*z-.06136358)*z+1.4138478)*z
 -11.36433272)*z+21.42034017)*z-3.69113734)*y;

 *keip=al* beip- bei/x-.7853981634* berp
 +x*((((((.00011997*z-.00926707)*z+.33049424)*z-4.65950823)*z
 +19.41182758)*z-13.39858846)*z+.21139217);
 return;
   }
/*else*/

CMPLX(CC,-.7071067812,-.7071067812);
CTREAL(CC,CC,x);
CTHET(-x,&y,&z);
CMPLX(theta,y,z);
CADD(CC,theta,CC);
cexp(&CC,&cex);
y=1.253314137/sqrt(x);
CTREAL(f,cex,y);
*ker-f.x;
*kei-f.y;
CPHI(-x,&y,&z);
CMPLX(phi,y,z);/*cex now phi of AS*/
CMULT(cdum,f,phi);
*kerp=-cdum.x;
*keip=-cdum.y;
return;
}

be(X,BER,BEI,BERP,BEIP)double X,*BER,*BEI,*BERP,*BEIP;
{
struct complex CC,cex,cdum,g,theta,phi;
double Y,Z,sqrt(),FKER,FKEI,FKERP,FKEIP;
static double PII=.3183098862;
if(X<=8.)
 {
 Y=(X/8.);Y=Y*Y;
 Z=Y*Y ;
 *BER=((((((-.00000901*Z+.00122552)*Z-.08349609)*Z
 +2.64191397)*Z-32.36345652)*Z+113.77777774)*Z-64.)*Z+1.;
 *BEI=Y*((((((.00011346*Z-.01103667)*Z+.52185615)*Z
 -10.56765779)*Z+72.81777742)*Z-113.77777774)*Z+16.);
 *BERP=X*((((((-.00000394*Z+.00045957)*Z-.02609253)*Z
 +.66047849)*Z-6.06814810)*Z+14.22222222)*Z-4.)*Y;
 *BEIP-X*((((((.00004609*Z-.00379386)*Z+.14677204)*Z
 -2.31167514)*Z+11.37777772)*Z-10.66666666)*Z+.5);
 return;
 }
/*else*/
   CMPLX(CC,.7071067812,.7071067812);
   CTREAL(CC,CC,X);
   CTHET(X,&Y,&Z);
   CMPLX(theta,Y,Z);
   CADD(CC,theta,CC);
   cexp(&CC,&cex);
   Y=.3989422804/sqrt(X);
   CTREAL(g,cex,Y);
   ke(X,&FKER,&FKEI,&FKERP,&FKEIP);
   *BER=g.x-PII*FKEI;
   *BEI=g.y+PII*FKER;
   CPHI(X,&Y,&Z);
   CMPLX(phi ,Y,Z);
   CMULT(cdum,phi,g);
   *BERP= cdum.x-PII*FKEIP;
   *BEIP= cdum.y+PII*FKERP;
   return;
}

CTHET(X,PARTR,PARTI)double X,*PARTI,*PARTR;
{
double Y;
 Y=8./X;
 *PARTI=((((.0000019*Y+.0000051)*Y*Y-.0000901)*Y-.0009765)
 *Y-.0110485)*Y-.3926991;
 *PARTR=((((.0000006*Y-.0000034)*Y-.0000252)*Y-.0000906)*Y*Y
 +.0110486)*Y;
return;
}
CPHI(X,PARTR,PARTI)double X,*PARTR,*PARTI;
{double Y;
  Y=8./X;
  *PARTI=(((((-.0000032*Y-.0000024)*Y+.0000338)*Y+.0002452)*Y
  +.0013811)*Y-.0000001)*Y+.7071068;
  *PARTR=((((((.0000016*Y+.0000117)*Y+.0000346)*Y+.0000005)*Y
  -.0013813)*Y -.0625001)*Y+.7071068);
  return;
}