/* ciureanu.c --- Efficiency and inductance of multi-turn head. Uses a
 *			three region model.
 *
 *	"Efficiency and Inductance Analysis of a Three-Region
 *	 Thin-Film Head", Petru C. Ciureanu, IEEE Trans. on Mag.,
 *	 Vol. MAG-22, No. 1, January 1988.
 *
 *	This model is an extension of the Jones model from two regions
 *	to three regions. The third region is a 'back gap' which will
 *	usually be smaller than the main part of the yoke. For the
 *	case in which this region is of zero length, the results are
 *	the same as Jones'.
 */
#include <stdio.h>
#include <math.h>
#include <malloc.h>

struct chead {
    double th;
    double lyk;
    double bg;
    double g;
    double ga;
    double t;
    double mur;
    double l;
    double *sp;
    double w;
    int N;
};

#define La	hd.th
#define Lb	hd.lyk
#define Lc	hd.bg
#define Ga	hd.g
#define Gb	hd.ga
#define T	hd.t
#define Mur	hd.mur
#define L	hd.l
#define S	hd.sp
#define W	hd.w
#define N	hd.N
#define mu0 1.2566370614e-6
void ciur2(struct chead hd,double *eff,double *induct)
{
	double lambda_a,lambda_b;
	double alpha,L1;
	double Laa,Lba,Lca,Lbb,L1b,L2b,Slaa,Claa,Slca,Clca,Slbb,Clbb;
	double Sl2b,tmp,K,K0,K1,K2,M1,M2;
	double SS1,SC1,SS2,SC2,SS3,SC3;
	int i,j;

/* Initialize returned parameters. */
	*eff = 0.0;
	*induct = 0.0;
	
/* Compute lambda_a, lambda_b and alpha using Eq.(12). */
/* (Note the assumption of uniform thickness and permeability.) */
	lambda_a = sqrt(Mur * T * Ga / 2.0);
	lambda_b = sqrt(Mur * T * Gb / 2.0);
	alpha = mu0 * Mur * T * W / ( 2.0 * L);

/* Compute auxiliary variables. */
	Laa = La / lambda_a;
	Lca = Lc / lambda_a;
	Lbb = Lb / lambda_b;
	L1b = L / lambda_b;
	L2b = L / (2.0 * lambda_b);
    Sl2b = sinh(L2b);
    Claa = cosh(Laa);
    Slaa = sinh(Laa);
    Clca = cosh(Lca);
    Slca = sinh(Lca);
    Clbb = cosh(Lbb);
    Slbb = sinh(Lbb);
	K = (Ga * lambda_b) / (Gb * lambda_a); 
	SS1 = 0.0;
	SC1 = 0.0;
	SS2 = 0.0;
	SC2 = 0.0;
	SS3 = 0.0;
	SC3 = 0.0;
	for (i = 0; i < N; i++) {
        SS1 += sinh(S[i] / lambda_b);
        SC1 += cosh(S[i] / lambda_b);
        SS2 += sinh((S[i] - La) / lambda_b); 
        SC2 += cosh((S[i] - La) / lambda_b); 
        SS3 += sinh((La + Lb - S[i]) / lambda_b);
        SC3 += cosh((La + Lb - S[i]) / lambda_b);
	}
/* Compute efficiency using Eqs.(15) and (16). */
	tmp = Claa*Clbb+ Slaa * Slbb / K;
	tmp *= (L2b * (Clca + K * Slca));
	K0 = Sl2b / tmp; 	

	*eff = K0 * (Clca * SC3 + K * Slca * SS3) / N;

/* Compute inductance using Eq.(23) and (24). */
	tmp = Clca + K * Slca;
	tmp *= (Claa * Clbb + Slaa * Slbb / K);
	L1 = 4.0 * alpha * Sl2b * Sl2b / (K * L1b * tmp);  

	M1 = 0.0;
	for (i = 0; i < N; i++) {
		tmp = 0.0;
		for (j = i; j < N; j++)
            tmp += sinh(S[j] / lambda_b);
        M1 += cosh(S[i] / lambda_b) * tmp;
	}
	M2 = 0.0; 		
	for (i = 0; i < N; i++) {
		tmp = 0.0;
		for (j = i; j < N; j++)
            tmp += cosh(S[j] / lambda_b);
        M2 += sinh(S[i] / lambda_b) * tmp;
	}
	tmp = Clca * Slaa * SC2 * SC3;
	tmp += (K * Clca * Claa * SS2 * SC3);
	tmp += (K * Slca * Slaa * SC2 * SS3);
	tmp += (K * K * Slca * Claa * SS2 * SS3);
	L = L1 * tmp;
	L += 4.0 * alpha * Sl2b * Sl2b * (M2 - M1) / L1b;

/* NOTE: In the following equation, the original paper shows '6' in
 *	 place of 'N'. This causes the computed inductance to differ
 *	 from that of Jones in cases where  lc = 0.0.
 */
    L += N * alpha * (1.0 - sinh(L1b) / L1b); 
	*induct = L;
}