/* bertero.c --- Implementation of Bertero, et al. head field model.
 *
 *		"Fields and Transforms for Thin Film Heads",
 *       Gerardo A. Bertero, H. Neal Bertram, and
 *		 David M. Barnett, IEEE Trans. on Mag.,
 *		 Vol. 29, No. 1, January 1993.
 *
 */
#include <math.h>
#include <stdio.h>

#define pi	3.14159265358979323846
#define pi_o_2	1.57079632679489661923
#define two_pi	6.2831853071795864769
#define sqrt2	1.4142135623730950488
#define HX	1
#define HY	2
#define HXY	3

int bertero(double g,double p,double q,double x,double y,
			double *Hx,double *Hy,int CMD)
{
	double g2,g22,a1,a2,a3,a4,a5,t1,t2,t3,t4,x2,y2;
	double L1,L2,pt1,pt2,pw1,pw2,sq,xyg1,xyg2;
	
/* Check for valid input options. */
	if (y < 0.0) return 1;
	if ((CMD == 0) || (CMD > 3)) return 1;
	if ((Hx == NULL) && (CMD & HX)) return 1;
	if ((Hy == NULL) && (CMD & HY)) return 1;

/* Compute auxiliary variables. */	
	L1 = 2.0 * p * g;
	L2 = 2.0 * q * g;
	g2 = g / 2.0;
	g22 = g2 * g2;
	x2 = x * x;
	y2 = y * y;
	xyg1 = x2 - y2 - g22;
	xyg2 = x2 + y2 - g22;	
	sq = sqrt(xyg1 * xyg1 + 4.0 * x2 * y2);
	pt1 = (pi_o_2 - atan((x + L1 * 0.5) / y)) / 3.0;
	pt2 = (pi_o_2 + atan((x - L2 * 0.5) / y)) / 3.0;
	pw1 = pow(L1 * L1 * 12.0 * pi,1.0 / 3.0); 
	pw1 *= pow((x + L1 * 0.5) * (x + L1 * 0.5) + y2,1.0 / 6.0);
	pw2 = pow(L2 * L2 * 12.0 * pi,1.0 / 3.0); 
	pw2 *= pow((x - L2 * 0.5) * (x - L2 * 0.5) + y2,1.0 / 6.0);
	

	if (CMD & HX) {
		*Hx = 0.0;
/* Compute Hx using Eq.(8) */
		a1 = atan((g2 + x) / y) + atan((g2 - x) / y);	
		t1 = a1 / two_pi;	

		a1 = sqrt((sq - xyg1) /(xyg2 * xyg2 + 4.0 * y2 * g22));
		t2 = a1 / (two_pi * sqrt2);

		a1 = sin(pt1);
		t3 = a1 / pw1;

		a1 = sin((pi_o_2 + atan((x - L2 / 2.0)/y)) / 3.0);
		t4 = a1 / pw2;

		*Hx = (t1 + g * (t2 - t3 - t4));
	}
	if (CMD & HY) {
		*Hy = 0.0;
/* Compute Hy using Eq.(9) */
		a1 = log(((x + g2) * (x + g2) + y2) /
			((x - g2) * (x - g2) + y2));
		t1 = 0.25 * a1 / pi;
		
		a1 = sqrt((sq + xyg1) /(xyg2 * xyg2 + 4.0 * y2 * g22));
		if (x >= 0.0) t2 = a1 / (two_pi * sqrt2);
		else t2 = -a1 / (two_pi * sqrt2);
		
		a1 = cos(pt1);
		t3 = a1 / pw1;

		a1 = cos(pt2);
		t4 = a2 / pw2;
		
		*Hy = (-t1 + g * (-t2 + t3 - t4)); 	
	}
	return 0;
}