#include <math.h>
/*	wandc.c --- Williams & Comstock write model.
 *			by C. Bond, 1994
 *
 *	An implementation of "An Analytical Model of the Write
 *	Process in Digital Recording Media", by M.L. Williams
 *	and R.L. Comstock, 17th Annual AIP Conference
 *	Proceedings, Part 1, No. 5, 1971, pp. 738-742.
 *
 *	This model is used to calculated the so-called 'a'
 *	parameter which is a measure of transition length
 *	in the media. The effective transition length is
 *	influenced by the write field gradient, the media
 *	properties and image fields. It is assumed that
 *	the write current has been adjusted for the given
 *	fly height, etc., so that the head field gradient
 *	is maximum. An arctangent transition shape is
 *	also assumed.
 *
 *	An input option, SHUNT/NOSHUNT, is provided to
 *	select transition lengths with or without the
 *	effects of shunting by the head.
 *
 *	The central equation governing the model is:
 *
 *              dM(I)           dM(I)   dH(h)   dH(dI)
 *              ----- (x=0) =   ----- ( ----- + ------ )
 *               dx             dH(I)    dx       dx
 *
 *	where the terms and relations are explained in the
 *	above cited paper.
 */
 
#define NOSHUNT	0
#define SHUNT	1

/* Input parameters are:
 *	Hc	- coercivity of the media,
 *	Mr	- media remanence,
 *  d   - total head to media spacing (bottom of head to top of media)
 *	delta	- media thickness,
 *  SS  - s-star, W&C media squareness parameter,
 *  g   - head gap length,
 *	option	- selection command for SHUNT or NOSHUNT
 *
 * Output values are:
 *  errcode - the function return value: 0 = success 
 *	a1	- transition length for head field (intermediate),
 *	a2	- transition length due to demagnetization fields,
 *	ad	- media dominant transition length (for reference).
 */
 
int wandc(double Hc,double Mr,double d,double delta,
	double SS,double g,double *a1,double *a2,double *ad,int option) 
{
	double CHI,r,y,Q,K1,a,A1,A2,A3,tmp;
	double f,V,B,rr;

/* First initialize return values. */
	*a1 = 0.0;
	*a2 = 0.0;
	
/* Now test for invalid input parameters. */
	if ((Hc == 0) || (g == 0)) return 1;

/* Compute auxiliary variables. */
	CHI = Mr / (4.0 * Hc);
	r = (3 + SS) / 4.0;
	y = sqrt(d * (d + delta));

/* Compute media dominant transition length for reference. W&C
 * recommend using this for the transition length if it is
 * larger than a2. */
	*ad = 2 * Mr * delta * (3.3 - 2.3 * SS) / Hc;
	
/* Compute Q using simplified quadratic. NOTE: the value of Q used in
 * the following equations depends on the assumption of maximum head
 * field gradient in the media. The relations are determined from the
 * Karlquist head field model and are open to challenge. A better
 * model would follow from replacing the Q parameter with a value
 * determined from a calculation of the actual head field gradient.
 * The following equation is good for the range :    0 << y/g <= 2 
 */

 	Q = 0.65 + (y / g) * ( 0.2 - 0.05 * (y / g));

/* Other equations for Q are: (see exact solution elsewhere)
 *	1) Q = 0.7115 + (y / g) * (0.3265 - 0.114 * (y / g))
 *		(best fit to published curve)
 *	2) Q = 0.65 + 0.1 * y / g
 *		(straight line)
 *	3) Q = 0.65 + (y/g) * (0.375 - 0.175 * (y/g))
 *		(2nd order fit)
 */
 
/* Compute simplified intermediate transition parameter, a1, using Eq(2). */
 	K1 = y * (1.0 - SS) / (pi * Q);
 	A1 = r * K1 + r * sqrt(K1 * K1 + (4.0 * Mr * delta * y) / (Hc * Q * r));

	if (option == NOSHUNT) {
/* Compute simplified 'final' transition, a2, without including the
 * shunting effect of the head on the demagnetization field, using Eq(4). */
		tmp = sqrt((A1 / (2.0 * r)) * (A1 / (2.0 * r)) + 
			2.0 * pi * CHI * delta * (A1 / r)); 	 	
		A2 = A1 / (2.0 * r) + tmp;
		*a1 = A1;
		*a2 = A2;
		return 0;
	}

/* Compute modified a1 with shunting (imaging) effects using the simplified
 * value of a1 as an estimate and bisection on Eq(5). */
	a = A1;
	rr = a;
	while (rr > 1e-8) {
		rr /= 2.0;
		V = 4.0 * y * (a + y) / ((a + 2.0 * y) * (a + 2.0 * y));
		B = (1.0 - SS) * y / (pi * Q);
		f = B + sqrt(B * B + 4.0 * Mr * delta * y * V / 
			(Hc * Q * r )) - a / r;
		if (f > 0.0) a += rr;
		else a -= rr;
	}
	A1 = a;
	
/* Compute the final transition, a2, from the modified a1. */
	tmp = sqrt((A1 / (2.0 * r)) * (A1 / (2.0 * r)) + 
		2.0 * pi * CHI * delta * (A1 / r)); 	 	
	A2 = A1 / (2.0 * r) + tmp;
	*a1 = A1;
	*a2 = A2;	
	return 0;
}