/*
*
* ELLIPTIC.C: Functions for computing complete and incomplete elliptic
* integrals of the first and second kind, and an example
* of usage for computing magnetic fields due to a circular
* current filament.
*
* DISCLAIMER
*
* Although every effort has been expended to ensure correctness, this
* software is not guaranteed to be correct.
*
* You may do anything you want with this software, as long as you maintain
* the existing comments intact, including the names of the original
* authors. If you modify this software, please indicate that fact in
* your comments.
*
* E. Dennison (Nov 17, 1998)
*
*/ |
#include "float.h"
#include "math.h"
#include "stdio.h"
#define MAX(x,y) (x>y?x:y)
#define MAX3(x,y,z) MAX(MAX(x,y),z)
#define MIN(x,y) (x>y?y:x)
#define MIN3(x,y,z) MIN(MIN(x,y),z)
/*
* NOTE: constants declared in float.h for the x86 CPU are:
*
* DBL_EPSILON 2.2204460492503131e-016,smallest such that 1.0+DBL_EPSILON!=1.0
* DBL_MAX 1.7976931348623158e+308 maximum possible value for type double
* DBL_MIN 2.2250738585072014e-308 minimum possible positive value for double
*
* If you are compiling this code for a non-Intel CPU, your values for these
* constants may be different.
*/
/*
* F(k): complete elliptic integral of the first kind
*/
#define F(k,ierr) drf(0.0,1.0-pow(k,2),1.0,&ierr)
/*
* drf.c: Compute the complete or incomplete elliptic integral of the
* first kind.
*
* Description:
*
* For x,y and z non-negative and at most one of them zero, drf(x,y,z)
* = integral from zero to infinity of
*
*
* -1/2 -1/2 -1/2
* (1/2)(t+x) (t+y) (t+z) dt.
*
* If x, y or z is zero, the integral is complete.
*
* *piErr returns non-zero if any arguments are invalid.
*
* Credits:
*
* This function is adapted by E. Dennison from FORTRAN code written by:
*
* Carlson, B.C.
* Notis, E.M
* Pexton, R.L.
*
*/
|
double drf(double x, double y, double z, int* piErr)
{
int iErr=0;
double mu,
xn,yn,zn,
xndev,yndev,zndev,
xnroot,ynroot,znroot,
lambda,
epsilon,
e2,e3,
result,
s;
const double c1 = 1.0/24.0;
const double c2 = 3.0/44.0;
const double c3 = 1.0/14.0;
const double errtol = pow(DBL_EPSILON*4.0,1.0/6.0);
const double lolim = 5.0*DBL_MIN;
const double hilim = DBL_MAX/5.0;
if (piErr)
{
if (MIN3(x,y,z)<0.0)
{
iErr = 1;
}
else if (MIN3(x+y,x+z,y+z)<lolim)
{
iErr = 2;
}
else if (MAX3(x,y,z)>hilim)
{
iErr = 3;
}
}
if (iErr)
{
if (piErr)
{
*piErr = iErr;
}
result = 0.0;
} |
else
{
xn = x;
yn = y;
zn = z;
while (1)
{
mu = (xn+yn+zn)/3.0;
xndev = 2.0-(mu+xn)/mu;
yndev = 2.0-(mu+yn)/mu;
zndev = 2.0-(mu+zn)/mu;
epsilon = MAX3(fabs(xndev),fabs(yndev),fabs(zndev));
if (epsilon < errtol) break;
xnroot = sqrt(xn);
ynroot = sqrt(yn);
znroot = sqrt(zn);
lambda = xnroot*(ynroot+znroot) +ynroot*znroot;
xn = (xn+lambda)*0.25;
yn = (yn+lambda)*0.25;
zn = (zn+lambda)*0.25;
}
e2 = xndev*yndev - pow(zndev,2);
e3 = xndev*yndev*zndev;
s = 1.0 + (c1*e2-0.1-c2*e3)*e2 + c3*e3;
if (piErr)
{
*piErr = 0;
}
result = s/sqrt(mu);
}
return result;
}
/* END drf() */
|
/*
* E(k): complete elliptic integral of the second kind
*/
#define E(k,ierr) (drf(0.0,1.0-pow(k,2),1.0,&ierr)\
-(pow(k,2)/3.0)*drd(0.0,1.0-pow(k,2),1.0,&ierr))
/*
* FastE(k): fast, complete elliptic integral of the second kind.
* Use this macro if the complete elliptic integral of
* the first kind was previously computed for the same
* value of k.
*/
#define FastE(F,k,ierr) ((F)-(pow(k,2)/3.0)*drd(0.0,1.0-pow(k,2),1.0,&ierr))
/*
* drd.c: Compute the complete or incomplete elliptic integral of the
* second kind.
*
* Description:
*
* For x and y non-negative, x+y and z positive, drf(x,y,z) = integral
* from zero to infinity of
*
*
* -1/2 -1/2 -3/2
* (3/2)(t+x) (t+y) (t+z) dt.
*
* If x or y is zero, the integral is complete.
*
* *piErr returns non-zero if any arguments are invalid.
*
* Credits:
*
* This function is adapted by E. Dennison from FORTRAN code written by:
*
* Carlson, B.C.
* Notis, E.M
* Pexton, R.L.
*
*/
|
double drd(double x, double y, double z, int* piErr)
{
int iErr=0;
double mu,
xn,yn,zn,
xndev,yndev,zndev,
xnroot,ynroot,znroot,
lambda,
epsilon,
ea,eb,ec,ed,ef,
sigma,
power4,
result,
s1,s2;
const double c1 = 3.0/14.0;
const double c2 = 1.0/6.0;
const double c3 = 9.0/22.0;
const double c4 = 3.0/26.0;
const double errtol = pow(DBL_EPSILON/3.0,1.0/6.0);
double uplim;
const double lolim = 2.0/pow(DBL_MAX,2.0/3.0);
double tuplim = pow(DBL_MIN,1.0/3.0);
tuplim = pow(0.1*errtol,1.0/3.0)/tuplim;
uplim = pow(tuplim,2.0);
if (piErr)
{
if (MIN(x,y)<0.0)
{
iErr = 1;
}
else if (MAX3(x,y,z)>uplim)
{
iErr = 2;
}
else if (MIN(x+y,z)<lolim)
{
iErr = 3;
}
}
if (iErr)
{
if (piErr)
{
*piErr = iErr;
}
result = 0.0;
} |
else
{
xn = x;
yn = y;
zn = z;
sigma = 0.0;
power4 = 1.0;
while (1)
{
mu = (xn+yn+3.0*zn)*0.2;
xndev = (mu-xn)/mu;
yndev = (mu-yn)/mu;
zndev = (mu-zn)/mu;
epsilon = MAX3(fabs(xndev),fabs(yndev),fabs(zndev));
if (epsilon < errtol) break;
xnroot = sqrt(xn);
ynroot = sqrt(yn);
znroot = sqrt(zn);
lambda = xnroot*(ynroot+znroot) +ynroot*znroot;
sigma = sigma+power4/(znroot*(zn+lambda));
power4 = power4*0.25;
xn = (xn+lambda)*0.25;
yn = (yn+lambda)*0.25;
zn = (zn+lambda)*0.25;
}
ea = xndev*yndev;
eb = zndev*zndev;
ec = ea-eb;
ed = ea-6.0*eb;
ef = ed+ec+ec;
s1 = ed*(-c1+0.25*c3*ed-1.5*c4*zndev*ef);
s2 = zndev*(c2*ef+zndev*(-c3*ec+zndev*c4*ea));
if (piErr)
{
*piErr = 0;
}
result = 3.0*sigma+power4*(1.0+s1+s2)/(mu*sqrt(mu));
}
return result;
}
/* END drd() */
|
/*
* MAIN: Sample program using drd and drf
*
* Description:
*
* Demonstrate use of drd and drf in computing the magnetic
* field at any point in space due to a circular current
* filament (one turn coil).
*
* Author:
* E Dennison
*
* Disclaimer:
*
*/
int main()
{
#define PI (3.141592654)
#define MU0 (PI*4.0E-7)
/* input parameters */
float a; /* loop radius */
float r; /* measurement point radius */
float x; /* measurement point axial position */
float i; /* loop current */
/* output parameters */
double Hx,Hr; /* axial, radial field components, A/m */
/* working vars */
double k,q,rq,fk,ek,al,be,ga,al2,be2,alt4,Ht;
int ierr;
/* gather input parameters */
printf("Loop radius (meters): ");
scanf("%f",&a);
printf("Radius at measurement point (meters): ");
scanf("%f",&r);
printf("Axial position at measurement point (meters): ");
scanf("%f",&x);
printf("Loop current (amperes): ");
scanf("%f",&i);
|
/* begin computation here */
al = r/a;
alt4 = al*4.0;
be = x/a;
ga = x/r;
al2 = al*al;
be2 = be*be;
q = pow(1+al,2)+be2;
rq = sqrt(q);
k = sqrt(alt4/q);
fk = F(k,ierr);
ek = FastE(fk,k,ierr);
Ht = i/(2.0*a*PI*rq);
Hx = Ht*(ek*(1-al2-be2)/(q-alt4)+fk);
Hr = (r==0.0)?(0.0):(Ht*ga*(ek*(1+al2+be2)/(q-alt4)-fk));
/* display the output */
printf("Axial field, Bx: %e (Tesla)\n",(float)Hx*MU0);
printf("Radial field, Br: %e (Tesla)\n",(float)Hr*MU0);
return 0;
}
|