#include <math.h>
/* jones.c --- Implementation of the Jones 'transmission line' model.
 *		by C. Bond, 1994
 *
 *	"Analysis of the Efficiency and Inductance of Multiturn
 *	 Thin Film Magnetic Recording Heads", R. E. Jones, Jr.,
 *	IEEE Trans. on Magnetics, Vol. MAG-14, No. 5, Sept. 1978.
 *
 *	The following code uses multiple turns. Results from the
 *	efficiency calculations will generally not differ much from
 *	those using a single turn model (jones1).
 *
 *	It is assumed that all turns have the same width (= delta).
 *	This assumption leads to some simplifications in the loops
 *	used to compute efficiency and inductance. It also leads
 *	to a single value for 'alpha', further simplifying the loops.	
 *
 *	The determination of the coil geometry is handled externally.
 *	Spacing between coils may be estimated based on some fraction 
 *	of the trace (stripe) width, or on a fixed process limit. In
 *	the latter case, 2.5 microns is a good figure for current
 *	(1993) technology.
 */

/*
#define La	LYK	// length of yoke body
#define Lb	th	// throat height
#define ga	GA	// spacing between yoke top and bottom
#define gb	g	// gap length
#define t1a	YK1	// top yoke thickness
#define t2a	YK2	// bottom yoke thickness
#define t1b	p1t	// p1 thickness
#define t2b	p2t	// p2 thickness
// mu is the permeability
// W  is the (average) yoke width
// N  is the number of turns
// D1 is the distance from the back gap to the start of the coils
// D2 is the length of the coil region
*/

#define mu0	1.25663706144e-6

void jones(double t1b,double t2b,double gb,double Lb,int N,
	double delta,double D1, double D2,double ga,double t1a,
	double t2a,double La,double mu,double W,
	double *eff,double *induct)
{
	double alpha,Claa,Clbb,d1a,d2a,K,lambda_a,lambda_b;
	double Laa,Lbb,lsq,M,SCla,Sd1a,Sd2a,Slaa,Slbb,SSla;
	double T1,T2,tmp,tmp2,Si,Sj;
	int i,j;

/* Initialize returned parameters */
	*eff = 0.0;
	*induct = 0.0;
	
/* Compute alpha using Eqs.(8) and (9) */
	lsq = ga / (1.0 / (mu * t1a) + 1.0 / (mu * t2a));
	alpha = lsq * mu0 * W / delta / ga;

/* Compute lambda_a using Eq.(8) */
	lambda_a = sqrt(lsq);

/* Compute lambda_b using Eq.(14) */
	lambda_b = sqrt(gb / (1.0 / (mu * t1b) + 1.0 / (mu * t2b)));

/*****************************
 *	Auxiliary variables
 *****************************/
	Laa = La / lambda_a;
	Lbb = Lb / lambda_b;
	d1a = delta / lambda_a;
	d2a = d1a * 0.5;
	K = ga * lambda_b / ( gb * lambda_a);

	Slaa = sinh(Laa);
	Claa = cosh(Laa);
	Slbb = sinh(Lbb);
	Clbb = cosh(Lbb);
	Sd1a = sinh(d1a);
	Sd2a = sinh(d2a);

/* Compute sums */
	SCla = 0.0;
	SSla = 0.0;
	for (i = 0; i < N; i++) {
		Si = D1 + (double)i * D2 / (double)N;
		SCla += cosh(Si / lambda_a);
		SSla += sinh(Si / lambda_a);
	}
	M = 0.0;
	for (i = 0;i < N; i++) {
		tmp = 0.0;
		tmp2 = 0.0;
		for (j = i; j < N; j++) {
			Sj = D1 + (double)j * D2 / (double)N;
			tmp  += cosh(Sj / lambda_a);
			tmp2 += sinh(Sj / lambda_a);
		}
		Si = D1 + (double)i * D2 / (double)N;
		M += (sinh(Si / lambda_a) * tmp 
			- cosh(Si / lambda_a) * tmp2);
	}
	
/* Compute T1 and T2 using Eqs.(22) and (23) */
	T1 = Slaa * Clbb + K * Claa * Slbb;
	T2 = Claa * Clbb + K * Slaa * Slbb;		

/* Compute efficiency using Eqs.(28) and (29) */		
	*eff = Sd2a * SCla / (d2a * N * T2);

/* Compute inductance using Eqs.(31) */
 	tmp = Sd2a * Sd2a *(M + SCla*(SCla * (T1 / T2) - SSla));
	tmp *= (4.0 * lambda_a);
	tmp += (N*(delta - lambda_a * Sd1a));
	*induct = alpha * tmp / delta;		
}