/* lindholm.c -- Implementation of Lindholm's field equations.
 *
 *	"Magnetic Fields of Finite Track Width Heads",
 *	Dennis A. Lindholm, IEEE Trans. on Mag., Vol. MAG-13,
 *	No. 5, September 1977.
 *
 *	These equations use superposition to provide approximate
 *	solutions for head fields of finite track width heads.
 *	Exact solutions for `wedge' sections, previously derived
 *	by van Herk and Tjaden for the zero gap case, are used
 *	with the Karlquist approximation for the gap potential
 *	to derive equations for the X, Y, and Z fields.
 *
 *	This approach is not accurate for near field calculations
 *	because the surface charge distribution under material
 *	superposition is not preserved. Volume charge, however,
 *	is maintained so far field calculations should be
 *	reasonably accurate. The Karlquist gap potential is also
 *	not accurate in the near field situation.
 *
 *	In spite of the cautions above, these equations are
 *	easy to implement and provide reasonably close results
 *	for many situations.
 */
#include <float.h>
#include <math.h>

#define pi	3.14159265358979323846

double acoth(double x)
{
    return 0.5*log((x+1.0)/(x-1.0));
}

double asinh(double x)
{
    return log(x+sqrt(x*x+1.0));
}

int lindholm(double g, double w, double x, double y, double z,
double *Hx, double *Hy, double *Hz)
{
	double up,um,uupp,uupm,uumm,uump,rhom,rhop,vm,vp,zm,zp;
	double H2x,H2y,H2z,f3xp,f3xm,f3yp,f3ym,f3zp,f3zm;
	double H3xm,H3xp,H3ym,H3yp,H3zm,H3zp;
	double ot,g2,tmp;
	
/* Check for valid inputs */
	if ((g < 0.0) || (y < 0.0) || (w < 0.0)) return 1;
	if (g == 0.0) g = 1e-6;
	if (y == 0.0) y = 1e-6;
	if (w == 0.0) w = 1e-6;
    if (z == 0.0) z = 1e-6;
/*  exception occurs if z approx w */
    if (fabs(z)-fabs(w) < 0.01) z += 0.01;
/* Initialize return values */
	*Hx = 0.0;
	*Hy = 0.0;
	*Hz = 0.0;

/* Initialize auxiliary variables */
	ot = 1.0 / 3.0;
	g2 = g / 2.0;	
	up = (x + g2) / y;
	um = (x - g2) / y;
	 
	zp = z + w / 2.0;
	zm = z - w / 2.0;
	rhop = sqrt(y * y + zp * zp);
	rhom = sqrt(y * y + zm * zm);
	vm = (zm < 0.0) ? (-atan(y / zm)) : (pi - atan(y / zm));
	vp = (zp < 0.0) ? (pi + atan(y / zp)) : (atan(y / zp));
	vm *= ot;
	vp *= ot;
    uupp = ot * asinh((x + g2) / rhop);
    uupm = ot * asinh((x + g2) / rhom);
    uump = ot * asinh((x - g2) / rhop);
    uumm = ot * asinh((x - g2) / rhom);

/* Compute field values */

	H2x = atan(up) - atan(um);
	H2y = -0.5 * (log(up * up + 1.0) - log(um * um + 1.0));
	H2z = 0.0;

	tmp = sin(2.0 * vm);
    f3xp = atan(sinh(2.0 * uupm) / tmp);
    f3xm = atan(sinh(2.0 * uumm) / tmp);
	H3xm = f3xp - f3xm;


	tmp = sin(2.0 * vp);
    f3xp = atan(sinh(2.0 * uupp) / tmp);
    f3xm = atan(sinh(2.0 * uump) / tmp);
	H3xp = f3xp - f3xm;

	tmp = cos(vm);
    f3yp = -4.0 * cosh(uupm) * tmp + acoth(cosh(uupm) / tmp);
    f3ym = -4.0 * cosh(uumm) * tmp + acoth(cosh(uumm) / tmp);
	H3ym = f3yp - f3ym;

	tmp = cos(vp);
    f3yp = -4.0 * cosh(uupp) * tmp + acoth(cosh(uupp) / tmp);
    f3ym = -4.0 * cosh(uump) * tmp + acoth(cosh(uump) / tmp);
	H3yp = f3yp - f3ym;

	tmp = sin(vm);
    f3zp = -4.0 * cosh(uupm) * tmp + acoth(cosh(uupm) / tmp);
    f3zm = -4.0 * cosh(uumm) * tmp + acoth(cosh(uumm) / tmp);
	H3zm = f3zp - f3zm;

	tmp = sin(vp);
    f3zp = -4.0 * cosh(uupp) * tmp + acoth(cosh(uupp) / tmp);
    f3zm = -4.0 * cosh(uump) * tmp + acoth(cosh(uump) / tmp);
	H3zp = f3zp - f3zm;

	*Hx = H3xm + H3xp - H2x;
	*Hy = H3ym + H3yp - H2y;
/*	*Hz = H3zm + H3zp - H2z;	Line below is corrected. */
	*Hz = -H3zm + H3zp - H2z;	
	return 0;
}