/*	nmmm.c --- recording model to determine 'a' parameter.
 *
 *	Implementation of "Modelling Digital Recording in Thin Film
 *	Media", by R.H. Noyau, B.K. Middleton, J.J. Miles, and
 *	N.D. Mackintosh, IEEE Transactions on Magnetics, Vol. 24,
 *	No. 6, Nov. 1988, pp. 2494-2496.
 *
 *	This model, in contrast to the Williams & Comstock model,
 *	does not assume that the operating write current has been
 *	optimized for maximum head field gradient in the media.
 *	Instead, it allows for operation over a range of currents
 *	and their corresponding head field values by approximating
 *	the B-H loop of the media. In addition, the initial
 *	magnetization curve is approximated in a way that provides
 *	information on minor loop operating points.
 */

#include <stdio.h>
#include <float.h>
#include <math.h>

#define pi	M_PI

/* Input parameters are:
 *	Hc	- coercivity of media,
 *	Mr	- media remanence,
 *	Hg	- deep gap field,
 *	d	- total head to media spacing,
 *	s	- squareness (Mr / Ms),
 *	D	- media thickness.
 *
 * Output values are:
 *	errcode - function return value,
 *	a1	- transition length due to head field gradient,
 *	a2	- transition length after demagnetization,
 *	a3	- transition length after remagnetization.
 */
int nmmm(double Hc,double Mr,double Hg,double d,double s,double D,
	double *a1,double *a2,double *a3)
{
	double CHI,Hc1,Hd,Ht,kinv,K,K1,K2,Mr1,Mh,Ms;
	double a,A1,A2,A3,f,r,x1,y,tmp,ps;
	

/* First initialize returned parameter values. */
	*a1 = 0.0;
	*a2 = 0.0;
	*a3 = 0.0;
	
/* Now check for bad input parameter. */
	if (Mr == 0.0) return 1;	/* Avoid possible divide by zero. */
	if (s == 0.0) return 1;

/* Compute auxiliary quantities. */
	ps = pi * s;
	Ms = Mr / s;
	CHI = (Mr / Ms) * sin(ps) / ps;
	y = d + D / 2.0;	

/* Compute the value of K=K1 to establish the hysteresis loop from
 * Eq(3) in the article. (Note negative Hg as in section 2.2 of article.) */
	tmp = atan((-Hg / Hc + 1.0) * tan(ps / 2.0)) 
		- atan((-Hg / Hc - 1.0) * tan(ps / 2.0));
	kinv = 1.0 + (1.0 / (ps)) * tmp;
	K1 = 1.0 / kinv;
	
/* Compute Mr1 from Eq(2b) */
	Mr1 = Mr * (2.0 * K1 - 1.0);
	
/* Compute Hc1 from Eq(2a) using M(h)=Mh=0 */
	Mh = 0.0;
	Hc1 = Hc * (1.0 + tan(((K1 - 1.0) * Mr + Mh) / 
		((2.0 * K1 / pi) * Ms)) / tan(ps / 2.0));

/* Compute Ht from Eq(2a) using M(h)=Mh=Mr1 / 2 */
	Mh = Mr1 / 2.0;
	Ht = Hc * (1.0 + tan(((K1 - 1.0) * Mr + Mh) / 
		((2.0 * K1 / pi) * Ms)) / tan(ps / 2.0));

/* Compute x1 from Eq(7). */
	x1 = y * tan((pi / 2.0) * (2.0 * Hc1 - Hg) / Hg);

/* Compute a1 using wide gap approximation --- Eq(9b) */
	a = .4;		/* starting value for iteration */
	r = a;
	while (r > 1e-8) {
		r /= 2.0;
		f = (Hg / pi) * (pi / 2.0 + atan((x1 + a) / y))
			- 2.0 * Mr1 * D / a - Ht;
		if (f < 0.0) a += r;
		else a -= r;
	}
	A1 = a;
/* Compute K2 to account for relaxation of the transition under the
 * influence of its own field, as the transition moves away from
 * the head. Uses Eq(12) */
 	K2 = ((Mr1 / 2.0) + Mr) / ((2.0 * Ms / pi) * atan((Ht / Hc + 1.0) 
 		* tan(ps / 2.0))+Mr);

/* Now compute a2 using Eq(15) */
 	a = 2.0*A1;		/* Initial estimate of a2 */		
	r = a;
	while (r > 1e-8) {
		r /= 2.0;
		
		Hd = -4.0 * Mr1 * D * A1 / (a * a + A1 * A1);
		f = (2.0 * Ms / pi) * K2 * atan((Hd / Hc + 1.0) *
			 tan(ps / 2.0)) - (1.0 - K2) * Mr;
		f -= (2.0 * Mr1 / pi) * atan(A1 / a);
		if (f < 0.0) a += r;
		else a -= r;
	}
	A2 = a;
	
/* Finally compute a3 using Eq(21). Note that this is the only use of CHI. */
	a = A2;
	r = a;
	while (r > 1e-8) {
		r /= 2.0;
		tmp = A1 / ((a + 2.0 * y) * (a + 2.0 * y) + A1 * A1) +
			A1 / (A2 * A2 + A1 * A1);
		f = CHI * 4.0 * Mr1 * D * (-A1 / (a * a + A1 * A1) + tmp);
		f -= (2.0 * Mr1 / pi) * (atan(A1 / a) - atan(A1 / A2));
		if (f < 0.0) a += r;
		else a -= r;
	}
	*a3 = a;
	*a2 = A2;
	*a1 = A1;
	return 0;
}