#include <stdio.h>
#include <math.h>
#include "maxdefc.h"

extern double bx1d[NZMAX+2], by1d[NZMAX+2], bz1d;
extern double bx2d[NNMAX], by2d[NNMAX];    /*bx2d(NX,NY), by2d(NX,NY)*/
extern int NX, NY, NZ;
extern double lenx, leny, lenz;
extern double delx, dely, delz;
extern double b0, dbob;
extern int nkmin1d, nkmax1d;
extern double alpha1, lambda_parx, lambda_pary;
extern double eslab, fracx;
extern int nxhalf, nyhalf;
extern double kmin, kmax, xi;
extern double alpha2, lambda_per;
extern double e2d;
extern int iseed_1, iseed_2, iseed_3;
extern int intinput[100];
extern double dbleinput[100];
extern int printb;
extern double alpha3;
extern double rr0, th0, ph0;

#define PI (4.0*atan(1.0))

void emfldspar3(t,x,y,z,bx,by,bz,nStts)
int *nStts;
double *t, *x, *y, *z, *bx, *by, *bz;
{
  double t1;
  double xx, yy, zz;
  double bxslab, byslab;
/*  double exslab,eyslab,ezslab; */
  double bx2d2d, by2d2d;
  void mhdfldsslab();
  void mhdflds2d();
  void lint();
  void lintarea();

  printb  = intinput[1];
  iseed_1 = intinput[2];
  iseed_2 = intinput[3];
  iseed_3 = intinput[4];
  nkmin1d = intinput[5];
  nkmax1d = intinput[6];
  nxhalf  = intinput[7];
  nyhalf  = intinput[8];

  eslab       = dbleinput[1];
  fracx       = dbleinput[2];
  lambda_parx = dbleinput[3];
  lambda_pary = dbleinput[4];
  alpha1      = dbleinput[5];

  e2d         = dbleinput[6];
  lambda_per  = dbleinput[7];
  alpha2      = dbleinput[8];

  xi          = dbleinput[9];
  kmin       = dbleinput[10];
  kmax       = dbleinput[11];

  alpha3      = dbleinput[12];
  th0         = dbleinput[13];
  ph0         = dbleinput[14];
  rr0         = dbleinput[15];

  printf("iseed1,iseed2,iseed3,nkmin,nkmax = %d, %d, %d, %lf, %lf\n",iseed_1,iseed_2,iseed_3,kmin,kmax);
  printf("eslab=%lf, fracx=%lf\n",eslab,fracx);

  if (*x < 0.0) *nStts = 0;
  if (*y < 0.0) *nStts = 0;
  if (*z < 0.0) *nStts = 0;

  t1 = 0.0;

  mhdfldsslab(&t1);
  mhdflds2d(&t1);

  zz = *z/delz; 
  xx = *x/delx;
  yy = *y/dely;

  if(zz > NZ) nStts = 0; 
  if(zz < 0) nStts = 0;

  lint(bx1d,NZ,zz,bxslab);
  lint(by1d,NZ,zz,byslab);

  lintarea(bx2d,NX,NY,xx,yy,bx2d2d);
  lintarea(by2d,NX,NY,xx,yy,by2d2d);

  if (eslab == 0.0) {
    bxslab = 0.0;
    byslab = 0.0;
  }
  if (e2d == 0.0) {
    bx2d2d = 0.0;
    by2d2d = 0.0;
  }

  *bx = (bxslab+bx2d2d)*rr0*rr0/((*z)*(*z));
  *by = (byslab+by2d2d)*rr0*rr0/((*z)*(*z));
  *bz = b0*rr0*rr0/((*z)*(*z));

  return;
}
