#include <math.h>
/* szczech.c --- Implementation of the Szczech head field function.
 *                 (C)1994, C. Bond, all rights reserved.
 *
 *	 Citations:
 *		 1) "An Approach for Deriving Field Equations for Magnetic
 *		  Heads of Different Geometrical Configurations",
 *		  Theodore J. Szczech and Paul R. Iverson, IEEE Trans.
 *		  on Mag., Vol. MAG-22, No. 5, September 1986.
 *
 *		 2) "Improvement of the Coefficients in Field Equations
 *		  for Thin-Film Recording Heads", T. J. Szczech and
 *		  P. R. Iverson, IEEE Trans. on Mag., Vol. MAG-23, No. 5,
 *		  pp. 3866-7, September 1987, Letters.
 *
 *		  This version uses equations from the papers cited above.
 *		  The values for the Hx and Hy fields are normalized to Hg
 *		  and the calling routine must multiply these by Hg to
 *		  recover the correct numeric values.
 *
 *	Undershoot scaling and reduction are added through parameters
 *	usp and usq. These parameters vary from 0.0 to 1.0 and control
 *	the relative amplitude of the corresponding undershoots. That is,
 *	usp = 0.0 removes the undershoot normally associated with 'p'.
 *	usp = 1.0 causes the head field to be calculated with full
 *	undershoot amplitude for pole 'p'. Values between 0 and 1 will
 *	scale the undershoots accordingly. Similarly for 'q'.	
 */
#define pi      3.14159265358979323846
#define pi_o_2	1.57079632679489661923
#define inv_pi	0.31830988618379067154
#define HX	1
#define HY	2
#define HXY	3

/* Szczech coefficients. Values are those recommended in 1986 paper. */
#define AA	0.835
#define BB	0.0433
#define CC	0.512

/* These (inferior) values are used in 1986 paper.
 * #define DD 	(0.008 * p / g + 0.116)
 * #define EE	(0.01 * p / g + 0.89)
 * #define FF	(0.008 * q / g + 0.116)
 * #define GG	(0.01 * q / g + 0.89)
 */

/* These (corrected) values are used in the 1987 paper. */ 
#define DD	(0.004 * p / g + 0.086)
#define EE	0.9
#define FF	(0.004 * q / g + 0.086)
#define GG	0.9
#define II	(0.004 * p / g + 0.086)
#define JJ	0.9
#define KK	(0.004 * q / g + 0.086)
#define LL	0.9

int szczech(double p,double q,double g,double x,double y,
		double usp,double usq,double *Hx,double *Hy,int CMD)
/* The parameter CMD determines which of the field values to calculate.
 * If the calling routine sets CMD to HX, only the HX field value
 * is calculated. If CMD is set to HY, only the Hy field value is
 * calculated. CMD may be set to HX|HY to calculate both field values.
 */
{
	double p1,q1,QQ,x1,y1,x2,y2;
	double a1,a2,a3,ax,t1,t2,t3,x1h1,x1h2;
	double lgc,p1h,q1h,p1hx,q1hx,atn;

/* Check for valid input options. */
	if ((y < 0.0) || (g <= 0.0)) return 1;
	if (y == 0.0) y = g * 1e-6;	/* make y 'effectively' zero */
	if (usp < 0.0) usp = 0.0;
	if (usq < 0.0) usq = 0.0;
	if (usp > 1.0) usp = 1.0;
	if (usq > 1.0) usq = 1.0;	
/* Compute x1,y1 using Eq.(18). */
	x1 = x / g;
	y1 = y / g;

/* Compute auxiliary variables. */
	y2 = y1 * y1;
	x2 = x1 * x1;
	x1h1 = 0.5 + x1;
	x1h2 = 0.5 - x1;
	atn = atan(x1h1 / y1) + atan(x1h2 / y1);
	p1 = 0.99 * p / g;
	q1 = 0.99 * q / g;
	p1h = p1 + 0.5;
	q1h = q1 + 0.5;	
	p1hx = p1h - x1;
	q1hx = q1h + x1;
	lgc = log((CC + 0.5) / (CC - 0.5)) / CC;
	
/* Compute QQ using Eq.(17). */
	QQ = 1.0 / ((y2 + (x1 + CC) * (x1 + CC)) * 
		(y2 + (x1 - CC) * (x1 - CC)));

	if (CMD & HX) {
		*Hx = 0.0;
/* Compute Hx field with Eq.(16). First part (t1) is modified Karlquist. */
		a1 = (x2 + y2 + CC * CC) * lgc;
		a2 = x1 * log((y2 + x1h2 * x1h2) / 
			(y2 + x1h1 * x1h1));
		a3 = QQ * BB * (x2 - y2 - CC * CC) - AA;
		t1 = QQ * y1 * BB * (a1 + a2) - a3 * atn; 

		ax = x1 - EE * p1h;
		a1 = DD / (ax * ax + y2);
		a2 = y1 * 0.5 * log(((p1h - EE * p1h) * (p1h - EE * p1h)) /
			(p1hx * p1hx + y2));
		a3 = ax * (pi_o_2 - atan(p1hx / y1));
		t2 = a1 * (a2 - a3);
	
		ax = x1 + GG * q1h;
		a1 = FF / (ax * ax + y2);
		a2 = y1 * 0.5 * log(((q1h - GG * q1h) * (q1h - GG * q1h)) /
			(q1hx * q1hx + y2));
		a3 = ax * (pi_o_2 - atan(q1hx / y1));
		t3 = a1 * (a2 + a3);
		*Hx = inv_pi * (t1 + usp * t2 + usq * t3);
	}
	if (CMD & HY) {
		*Hy = 0.0;
/* Compute Hy field with Eq.(19). */
		a1 = (x2 + y2 - CC * CC) * lgc;
		a2 = 2.0 * y1 * atn;
		a3 = log((y2 + x1h2 * x1h2) / (y2 + x1h1 * x1h1));
		a3 *= (0.5 * (QQ * BB * ( x2 - y2 - CC * CC) - AA));
		t1 = QQ * x1 * BB * (a1 + a2) + a3;
		
		a1 = y1 * (pi_o_2 - atan(p1hx / y1));
		a2 = log(((p1h - JJ * p1h) * (p1h - JJ * p1h)) /
			(p1hx * p1hx + y2));
		ax = x1 - JJ * p1h;
		a2 *= (0.5 * ax);
		a3 = II / (ax * ax + y2);
		t2 = a3 * (a1 + a2);

		a1 = y1 * (pi_o_2 - atan(q1hx / y1));
		a2 = log(((q1h - LL * q1h) * ( q1h - LL * q1h)) /
			(q1hx * q1hx  + y2));
		ax = x1 + LL * q1h;
		a2 *= (0.5 * ax);
		a3 = KK / (ax * ax + y2);
		t3 = a3 * (a1 - a2);
		*Hy = inv_pi * (-t1 - t2 + t3);
	}
	return 0;
}

int szpulse(double *array,int array_size,double sig_start,double sig_stop,
	double p,double q,double g,double usp,double usq,double y,int CMD)
/* Computes a vector representing an isolated pulse using Szczech head
 * field model.
 */
{
	double step_size;
	double  x,Hx,Hy;
	int i,errcode;
	
	step_size = (sig_stop - sig_start) / (double)array_size;
	x = sig_start;

	if (CMD & 2) {
		for (i=0;i<=array_size;i++) {
			if ((errcode = szczech(p,q,g,x,y,usp,usq,&Hx,&array[i],2)) != 0)
				return errcode;
			x += step_size;
		}
	return 0;
	}
	for (i=0;i<=array_size;i++) {
		if ((errcode = szczech(p,q,g,x,y,usp,usq,&array[i],&Hy,1)) != 0)
			return errcode;
/* the next 2 lines window the head field function */
			if (x < 0) array[i] *= (0.5*(cos(x*pi/sig_start)+1.0));
			else array[i] *= (0.5*(cos(x*pi/sig_stop)+1.0));
		x += step_size;
	}	
	return 0;  
}