// yhtfh.cpp -- Implementation of:
//
//  "Calculation of Exact Head and Image Fields of Recording
//   Heads by Conformal Mapping", by Jyh Shinn Yang and Huei
//   Li Huang, IEEE Trans. on Mag., Vol. 25, No. 3, May 1989.
//
//  The following code implements the thin film head field equations.
//
//  NOTE: The static global 't' can be used to carry convergence information
//  from one invocation to the next. Set 'flag' to zero at the 1st entry
//  to establish the initial solution. After that, set 'flag' to one so
//  that 't' can prime the iteration loop at the nearby last solution.
//
#include <float.h>
#include <math.h>
#include <complex.h>
#include <stdio.h>

#define max(a,b)    (((a) > (b)) ? (a) : (b))
#define pi M_PI
const double twopi = 2.0 * M_PI;
static complex <double> t;

int yhtfh(double g,double p,double x,double y,double *Hx,double *Hy,int FLAG)
{
    complex <double> a,A0,f,fp,H,sq,t0,tsq,tmpc,z,z0,r,theta,I;
	double c,lambda,lambda2,S,tmpr,z0r,z0i;
	int i;

// Check for valid input parameters.
    if ((g < 1.0e-4) || (p < 1.0e-4)) return 1; // Invalid gap or pole
    if (y < 0.0) return 1;                      // Invalid 'y'
    y = max(y,1e-6);

    S = (x < 0.0) ? -1.0 : 1.0;
    x = max(fabs(x),1e-8);      // Use positive x for iteration and
                                // Avoid zero slope at x = 0
// Initialize 
    *Hy = 0.0;                  // Clear result
	*Hx = 0.0;
    z = complex<double>(x,y);
    I = complex<double>(0.0,1.0);

    c = 1.0 + 2.0 * p / g;      // Always positive, real and > 1.0
    tmpr = sqrt(c * c - 1.0);   // Always positive real
    lambda = c + tmpr;          // Always positive, real and > 1.0
    A0 = complex<double>(0.0, -g / (twopi * lambda));
	z0r = p + g;
    z0i = - g * (c * log(c * c - 1.0 + c * tmpr) + log(tmpr)) / twopi;
    z0 = complex<double>(z0r,z0i); 

    tmpc = 2.0 * z / g;         // Initial estimate of t
	t0 = tmpc + sqrt(tmpc * tmpc - 1.0);
	tmpr = (1.0 + lambda * lambda) / 2.0;   
	if (FLAG == 0) t = t0;
	for (i = 0; i < 100; i++) { // Outer loop limits iteration count
        t = complex<double>(fabs(real(t)),imag(t));
		tsq = t * t;
        sq = I * sqrt((tsq - 1.0) * (lambda *lambda - tsq));
		f = sq - tmpr * log(sq + tsq-tmpr);
		f -= (lambda * log((sq + lambda) / tsq  - tmpr / lambda));
		f += ((z0 - z) / A0);
		fp = 2.0 * sq / t;
		a = f /  fp;
		t = t - a;
		if (sqrt(real(f)*real(f) + imag(f)*imag(f)) < 1e-4) goto _10;
	}
	fprintf(stderr,"Convergence failure!\n");
	return 2;
_10:
    if (real(t) < 0.0) t = complex<double>(fabs(real(t)),imag(t));
    tsq = t * t;
    sq = I * sqrt((tsq-1.0)*(lambda*lambda-tsq));

// *** Test for convergence to correct 'z'-value
//    f = sq-((1.0+lambda*lambda)/2.0)*log(sq+tsq-(1.0+lambda*lambda)/2.0);
//    f = f - lambda*log((sq+lambda)/tsq-(1.0+lambda*lambda)/(2.0*lambda));
//    z = A0*f+z0;
//    cout << z - A0*f+z0 << endl;
	
    H = -lambda / sq;
    *Hx = real(H);
    *Hy = -S * imag(H);
	return 0;
}