#include <math.h>

#define mu0 1.25663706144e-6

/* Head field function based on Karlquist model */
double Hx0(double Hg,double g,double x,double y)
{
    double g2,dtmp1;

    g2 = 0.5 * g;
    dtmp1 = atan((g2+x)/y)+atan((g2-x)/y);
    return Hg*dtmp1/M_PI;
}