498 lines
15 KiB
C
498 lines
15 KiB
C
/* 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__ */
|