/* t_update.f -- translated by f2c (version 20000817).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Subroutine */ int t_update__(p_a__, p_ih__, c_ih__, lpa, dm, da, isa,
                                helm,
                                f1h, f1v, f2h, f2v, f, ei, aki, ef, akf,
                                qhkl, en, hx, hy, hz, if1, if2, qm, ier)
real *p_a__, *p_ih__, *c_ih__;
logical *lpa;
real *dm, *da;
integer *isa;
real *helm, *f1h, *f1v, *f2h, *f2v, *f, *ei, *aki, *ef, *akf, *qhkl, *en,
    *hx, *hy, *hz;
integer *if1, *if2;
real *qm;
integer *ier;
{
  static doublereal dakf, daki, dphi;
  static integer ieri, imod, ieru;
  extern /* Subroutine */ int ex_up__();
  static doublereal df;
  static integer id;
  static real qs;
  static doublereal dbqhkl[3];
  extern /* Subroutine */ int sam_up__();
  static doublereal da2, da3, da4, da6;
  extern /* Subroutine */ int erreso_();
  static doublereal dda, def, dei, ddm, dqm, dhx, dhy, dhz, dqs;
  extern /* Subroutine */ int helm_up__(), flip_up__();

/* 	=================== */

/* dec$	  Ident	    'V01C' */
/* ------------------------------------------------------------------ */
/* 		Note: */
/* 		  IF1,If2 changed to Int*4. They were Real*4 in ILL version. */
/* ------------------------------------------------------------------ */
/* Updates: */
/*  V01A  7-May-1996 DM.	  Put error output to IOLUN, use IMPLICIT NONE and */
/* 			    get the code indented so that it is readable! */
/*  V01C 12-Oct-1998 DM.	  Put in Mag Field calculation for SINQ Helmolz coils. */
/* ----------------------------------------------------------------------- */
/*   Entry points in this file: */
/* 	T_UPDATE :  USE THE MOTOR ANGLES TO CALCULATE VALUES OF */
/* 			EI,KI,EF,KF,QH,QK,QL,EN AND TO DETERMINE */
/* 			WHETHER FLIPPERS F1,F2 ARE 'ON' OR 'OFF'. */
/* 	EX_UP    :  CALCULATE EX AND AKX FORM AX2 VALUES. */
/* 	SAM_UP	 :  CALCULATE THINGS IN RECIPROCICAL LATTICE. */
/* 	HELM_UP	 :  CALCULATE HELMHOLTZ FIELDS. */
/* 	FLIP_UP	 :  CHECK IF FLIPPERS ON OR OFF. */
/* ----------------------------------------------------------------------- */
/*  T_UPDATE: */
/* 	ROUTINE TO USE THE MOTOR ANGLES TO */
/* 	CALCULATE VALUES OF EI,KI,EF,KF,QH,QK,QL,EN */
/* 	USE CURRENT-SUPPLY VALUES TO COMPUTE HELMHOLTZ FIELDS HX,HY,HZ */
/* 	AND TO DETERMINE WHETHER FLIPPERS F1,F2 ARE 'ON' OR 'OFF'. */
/* 	FLIPPERS ARE 'ON' ONLY IF CURRENTS ARE THOSE GIVEN BY */
/* 	F1V,F1H,F2V,F2H */
/* ----------------------------------------------------------------------- */
/* File [MAD.SRC]T_UPDATE.FOR */
/* cc        IMPLICIT DOUBLE PRECISION (A-H,O-Z) */
/* 	include 'iolsddef.inc' */
/* ----------------------------------------------------------------------- */
/* 				    Define the dummy arguments */
/*  Input: */
/* Positions of angles A1-A6. */
/* 	.. helmotz (8 currents). */
/* Position of currents for flippers and .. */
/* Configuration of the machine ... */
/* Conversion factors for Helmotz (4 currents). */
/* True if machine in polarization mode. */
/* Monochr. d-spacing */
/* Analyser d-spacing */
/* 	.. +ve to the left. */
/* Scattering sense at analyser (probably) .. */
/* 	.. ki (probably). */
/* Angle between axis of Helmolz pair one and .. */
/*  .. current is possibly ki*f1h ??) */
/* Flipper 1 hor and vert currents (hor .. */
/*  .. current is possibly kf*f2h ??) */
/* Flipper 2 hor and vert currents (hor .. */
/*  Output: */
/* Energy unit */
/* Incident neutron energy */
/* Incident neutron wave vector */
/* Final neutron energy */
/* Final neutron wave vector */
/* Components of q in reciprocal space */
/* Energy transfer */
/* Components of Helmolz field on sample */
/* Status of Flippers 1 and 2 */
/* Length of q */
/* ----------------------------------------------------------------------- */
/* Error status */

/* ----------------------------------------------------------------------- */
/* 	SET UP */

  /* Parameter adjustments */
  --qhkl;
  --c_ih__;
  --p_ih__;
  --p_a__;

  /* Function Body */
  ddm = *dm;
  dda = *da;
  df = *f;
  da2 = p_a__[2];
  da3 = p_a__[3];
  da4 = p_a__[4];
  da6 = p_a__[6];
/* ----------------------------------------------------------------------- */

  ieri = 0;
  ieru = 9;
  ex_up__(&ddm, &dei, &daki, &da2, &df, &ieri);
  if (ieri == 0) {
    *ei = dei;
    *aki = daki;
    ieru = 0;
  } else {
    imod = 3;
    erreso_(&imod, &ieri);
    ieru = ieri + 8;
  }
  ieri = 0;
  ieru = 9;
  ex_up__(&dda, &def, &dakf, &da6, &df, &ieri);
  if (ieri == 0) {
    *ef = def;
    *akf = dakf;
    ieru = 0;
  } else {
    if (*isa == 0) {
      *ef = dei;
      *akf = daki;
      def = dei;
      dakf = daki;
      ieru = 0;
    } else {
      imod = 3;
      erreso_(&imod, &ieri);
      ieru = ieri + 8;
    }
  }
  if (*isa == 0) {
    *ef = dei;
    *akf = daki;
    def = dei;
    dakf = daki;
    ieru = 0;
  }
  if (ieru == 0) {
    *en = dei - def;
  }
  ieri = 0;
  ieru = 5;
  sam_up__(dbqhkl, &dqm, &dqs, &dphi, &daki, &dakf, &da3, &da4, &ieri);
  if (ieri == 0) {
    for (id = 1; id <= 3; ++id) {
      qhkl[id] = dbqhkl[id - 1];
    }
    *qm = dqm;
    qs = dqs;
    ieru = 0;
  } else {
    imod = 2;
    erreso_(&imod, &ieri);
    ieru = ieri + 4;
  }

  ieri = 0;
  if (*lpa) {
    ieru = 1;
    helm_up__(&dphi, helm, &p_ih__[1], &c_ih__[1], &dhx, &dhy, &dhz,
              &ieri);
    if (ieri != 0) {
/* 	    WRITE(6,*) 'ERROR IN HELM_UP CHECK COEFF' */
    } else {
      *hx = dhx;
      *hy = dhy;
      *hz = dhz;
      ieru = 0;
    }
    flip_up__(if1, if2, &p_ih__[1], f1v, f1h, f2v, f2h, aki, akf, &ieri);
  }
/* ----------------------------------------------------------------------- */
  *ier = ieru;
  return 0;
}                               /* t_update__ */


/* Subroutine */ int ex_up__(dx, ex, akx, ax2, f, ier)
doublereal *dx, *ex, *akx, *ax2, *f;
integer *ier;
{
  /* System generated locals */
  doublereal d__1;

  /* Builtin functions */
  double sin();

  /* Local variables */
  static doublereal arg;

/* 	================ */

/* ----------------------------------------------------------------------- */
/* 	 CALCULATE EX AND AKX FORM AX2 VALUES */
/* ----------------------------------------------------------------------- */
/* Part of [MAD.SRC]T_UPDATE.FOR */

/* ----------------------------------------------------------------------- */
/* 				    Define the dummy arguments */
/*  Input: */
/* D-SPACING OF CRISTAL */
/* TAKE OFF ANGLE */
/*  Output: */
/* ENERGY UNIT */
/* ENERGY */
/* MOMENTUM */
/* ----------------------------------------------------------------------- */
/* Error - IER=1 if DX OR AX TOO SMALL */
/* ----------------------------------------------------------------------- */
/* !!!!!!!!!! This has to be fixed manually after conversion by f2c. */
/* !!!!!!!!!! The reason is a different definition of the abs function. */
/* !!!!!!!!!! MK, May 2001 */
  arg = (d__1 = *dx * sin(*ax2 / 114.59155902616465), abs(d__1));
  if (arg <= 1e-4) {
    *ier = 1;
  } else {
    *ier = 0;
    *akx = 3.1415926535897932384626433832795 / arg;
/* Computing 2nd power */
    d__1 = *akx;
    *ex = *f * (d__1 * d__1);
  }
/* ----------------------------------------------------------------------- */
  return 0;
}                               /* ex_up__ */


/* Subroutine */ int sam_up__(qhkl, qm, qs, phi, aki, akf, a3, a4, ier)
doublereal *qhkl, *qm, *qs, *phi, *aki, *akf, *a3, *a4;
integer *ier;
{
  /* Builtin functions */
  double cos(), sin(), atan2();

  /* Local variables */
  static doublereal qpar, qperp;
  extern /* Subroutine */ int sp2rlv_();
  static doublereal qt[3];

/* 	================= */
/* ----------------------------------------------------------------------- */
/* Part of [MAD.SRC]T_UPDATE.FOR */

/* ----------------------------------------------------------------------- */
/* 				    Define the dummy arguments */
/*  Input: */
/* Actual value of KI */
/* Actual value of KF */
/* Actual value of A3 */
/*  Output: */
/* Actual value of A4 */
/* CALCULATED POSITIONS IN RECIP. LATTICE */
/* Q MODULUS */
/* Q MODULUS SQUARED */
/* ANGLE BETWEEN KI AND Q */
/* 				NOT CALCULATED) */
/* ----------------------------------------------------------------------- */
/* 				LOCAL PARAMETERS */
/* Error status. IER=1 if QM TOO SMALL (PHI */
/* ----------------------------------------------------------------------- */
/* 	QPAR AND QPERP ARE COMPONENTS OF Q PARALLEL AND PERP TO KI */
/* 	TRANSFORM FROM SCATTERING PLANE INTO RECIPROCAL LATTICE */

  /* Parameter adjustments */
  --qhkl;

  /* Function Body */
  qpar = *aki - *akf * cos(*a4 / 57.29577951308232087679815481410517);
  qperp = -(*akf) * sin(*a4 / 57.29577951308232087679815481410517);
  qt[0] = qpar * cos(*a3 / 57.29577951308232087679815481410517) + qperp *
      sin(*a3 / 57.29577951308232087679815481410517);
  qt[1] = -qpar * sin(*a3 / 57.29577951308232087679815481410517) + qperp *
      cos(*a3 / 57.29577951308232087679815481410517);
  qt[2] = 0.;
  sp2rlv_(&qhkl[1], qt, qm, qs, ier);
  *ier = 3;
  if (*qm > .001) {
    *ier = 0;
  }
  *phi = 0.;
  if (abs(qperp) > .001 && abs(qpar) > .001) {
    *phi = atan2(qperp, qpar);
  } else if (abs(qpar) < .001) {
    if (*a4 > (float) 0.) {
      *phi = -1.5707963267948966;
    } else {
      *phi = 1.5707963267948966;
    }
  }
/* ----------------------------------------------------------------------- */
  return 0;
}                               /* sam_up__ */


/* Subroutine */ int helm_up__(phi, helm, p_ih__, c_ih__, dhx, dhy, dhz,
                               ier)
doublereal *phi;
real *helm, *p_ih__, *c_ih__;
doublereal *dhx, *dhy, *dhz;
integer *ier;
{
  /* System generated locals */
  real r__1;
  doublereal d__1, d__2;

  /* Builtin functions */
  double sqrt(), atan2(), cos(), sin();

  /* Local variables */
  static doublereal hdir, hmod, hpar, hdir2, h__[4], hperp;
  static integer ic;
  static logical at_sinq__;

/* 	================== */

/* ---------------------------------------------------------------- */
/* 	    HELMHOLTZ COILS UPDATE  (ONLY IF LPA IS TRUE) */
/* ---------------------------------------------------------------- */
/* Part of [MAD.SRC]T_UPDATE.FOR */

/* 	include 'common_sinq.inc' */

/* ----------------------------------------------------------------------- */
/*  At ILL: */
/*  	There are 3 coils for Hx/Hy at 120 degrees to each other. */

/* 	There is a 4th coil for Hz. */

/*  At SINQ: */
/*  	There is an Hx coil and an Hy coil (actually each is 4 coils powered */
/*  	in series). They are mounted on a ring (SRO). The value of HELM is */
/* 	the angle between the Hx coil and ki. */

/* 	There is a 3rd coil for Hz. */
/* ----------------------------------------------------------------------- */
/* 				    Define the dummy arguments */
/*  Input: */
/* Angle between KI and Q in scatt'g plane (rad) */
/* Angle between first coil and KI (degrees) */
/* 	    with Helmolz coils */
/* Current values of 4 currents associated .. */
/*  Output: */
/* Conversion factor between Amperes and Gauss */
/* \ */
/*  > The calculated fields */
/* / */
/* ---------------------------------------------------------------- */
/* 				    LOCAL VARIABLES */
/* Error status. IER=1 if C_I HAS A ZERO COEF */
/* ---------------------------------------------------------------- */
/* 	H 	!FIELD OF 4 COILS AROUND SAMPLE */
/* 	HPAR 	!FIELD PAR COIL 1 */
/* 	HPERP	!FIELD PERP. TO COIL 1 */
/* 	HDIR2	!ANGLE BETWEEN FIELD AND COIL 1 */
/* 	HDIR	!ANGLE BETWEEN FIELD AND Q */
/* ---------------------------------------------------------------- */
  /* Parameter adjustments */
  --c_ih__;
  --p_ih__;

  /* Function Body */
  *ier = 0;
  at_sinq__ = TRUE_;
  if (!at_sinq__) {
    for (ic = 1; ic <= 4; ++ic) {
      h__[ic - 1] = 0.;
      if ((r__1 = c_ih__[ic], dabs(r__1)) > .001) {
        h__[ic - 1] = p_ih__[ic + 4] * c_ih__[ic];
      } else {
        *ier = 1;
      }
    }
    hpar = h__[0] - (h__[1] + h__[2]) * .5;
    hperp = sqrt(3.) * .5 * (h__[1] - h__[2]);
    *dhz = h__[3];
  } else {
    for (ic = 1; ic <= 3; ++ic) {
      h__[ic - 1] = 0.;
      if ((r__1 = c_ih__[ic], dabs(r__1)) > .001) {
        h__[ic - 1] = p_ih__[ic + 4] * c_ih__[ic];
      } else {
        *ier = 1;
      }
    }
    hpar = h__[0];
    hperp = h__[1];
    *dhz = h__[2];
  }

/* Computing 2nd power */
  d__1 = hpar;
/* Computing 2nd power */
  d__2 = hperp;
  hmod = sqrt(d__1 * d__1 + d__2 * d__2);
  if (abs(hpar) > .001 && abs(hperp) > .001) {
    hdir2 = atan2((abs(hperp)), (abs(hpar)));
    if (hpar > 0. && hperp < 0.) {
      hdir2 = -hdir2;
    } else if (hpar < 0. && hperp > 0.) {
      hdir2 = 3.1415926535897932384626433832795 - hdir2;
    } else if (hpar < 0. && hperp < 0.) {
      hdir2 += -3.1415926535897932384626433832795;
    }
  } else if (abs(hpar) > .001) {
    if (hpar >= (float) 0.) {
      hdir2 = (float) 0.;
    }
    if (hpar < (float) 0.) {
      hdir2 = 3.1415926535897932384626433832795;
    }
  } else if (abs(hperp) > .001) {
    if (hperp >= (float) 0.) {
      hdir2 = 1.5707963267948966;
    }
    if (hperp < (float) 0.) {
      hdir2 = -1.5707963267948966;
    }
  } else {
    hdir2 = (float) 0.;
  }
  hdir = hdir2 + *helm / 57.29577951308232087679815481410517 - *phi;
  *dhx = hmod * cos(hdir);
  *dhy = hmod * sin(hdir);

  return 0;
}                               /* helm_up__ */


/* Subroutine */ int flip_up__(if1, if2, p_ih__, f1v, f1h, f2v, f2h, aki,
                               akf,
                               ier)
integer *if1, *if2;
real *p_ih__, *f1v, *f1h, *f2v, *f2h, *aki, *akf;
integer *ier;
{
  /* System generated locals */
  real r__1, r__2;

/* 	================== */
/* ---------------------------------------------------------------- */
/* 	 FLIPPERS; ONLY IF LPA */
/* ---------------------------------------------------------------- */
/* Part of [MAD.SRC]T_UPDATE.FOR */
/* ----------------------------------------------------------------------- */
/* 				    Define the dummy arguments */
/*  Input: */
/* 	    with flipper coils */
/* Current values of 4 currents associated .. */
/*  .. current is possibly aki*f1h ??) */
/* Flipper 1 vert and hor currents (hor .. */
/*  .. current is possibly akf*f2h ??) */
/* Flipper 2 vert and hor currents (hor .. */
/* Incident neutron wave vector */
/*  Output: */
/* Final neutron wave vector */
/* Status of Flipper 1 (=1 if on and set OK) */
/* Status of Flipper 2 (=1 if on and set OK) */
/* ---------------------------------------------------------------- */
/* Error status. 0 if no error. */
  /* Parameter adjustments */
  --p_ih__;

  /* Function Body */
  *ier = 0;
  *if1 = 0;
  *if2 = 0;
  if ((r__1 = p_ih__[1] - *f1v, dabs(r__1)) < (float) .05
      && (r__2 = p_ih__[2] - *aki * *f1h, dabs(r__2)) < (float) .05) {
    *if1 = 1;
  }
  if ((r__1 = p_ih__[3] - *f2v, dabs(r__1)) < (float) .05
      && (r__2 = p_ih__[4] - *akf * *f2h, dabs(r__2)) < (float) .05) {
    *if2 = 1;
  }
/* ----------------------------------------------------------------------- */
  return 0;
}                               /* flip_up__ */